%macro mat_LowerDiagInv(LowDiagInput=, LowDiagInv=, dim=);
/**********************************************************************/
/* This macro computes the Cholesky decomposition for a square matrix */
/* Subsequently, the Cholesky decomposition is employed to compute */
/* the inverse of the input matrix. */
/* Note that if the original matrix is not full rank but of rank M, */
/* then the inverse matrix which is computed will be the inverse only */
/* through rows 1:M and columns 1:M such that the submatrix */
/* Cholesky[1:M, 1:M] has rank M. */
/* */
/* Macro variables: */
/* LowDiagInput - lower diagonal (square) matrix to be inverted */
/* LowDiagInv - inverse of LowDiagInput matrix */
/* Dim - number of rows in LowDiagInput */
/* */
/* Author: Dale McLerran */
/* Fred Hutchinson Cancer Research Center */
/* */
/**********************************************************************/

%global ScratchMatLDI;
%if &ScratchMatLDI=. %then
%let ScratchMatLDI = 1;
%else %let ScratchMatLDI = %eval(&ScratchMatLDI + 1);

array _scratch_&ScratchMatLDI {&dim,&dim} _temporary_;
do row=1 to &dim;
do col=1 to &dim;
_scratch_&ScratchMatLDI{row,col} = &LowDiagInput{row,col};
if row=col then &LowDiagInv{row,col} = 1;
else &LowDiagInv{row,col}=0;
end;
end;
do pivot=1 to &dim;
do row=pivot to &dim;
scale = _scratch_&ScratchMatLDI[row,pivot];
if scale>1E-10 then do;
do col=1 to &dim;
_scratch_&ScratchMatLDI{row,col} = _scratch_&ScratchMatLDI{row,col} / scale;
&LowDiagInv{row,col} = &LowDiagInv{row,col} / scale;
end;
if row>pivot then do;
do col=1 to &dim;
_scratch_&ScratchMatLDI{row,col} =
_scratch_&ScratchMatLDI{row,col} - _scratch_&ScratchMatLDI{pivot,col};
&LowDiagInv{row,col} = &LowDiagInv{row,col} - &LowDiagInv{pivot,col};
end;
end;
end;
end;
end;

%mend mat_LowerDiagInv;