I translated into Visual Basic the routine on page 91 of Numerical Recipes in Fortran 77 for computing the inverse of the Cholesky decomposition, and it produces the wrong result. I made a few minor changes because I produce the Cholesky dc as a separate lower triangular matrix C() rather than overwriting A() -- the original matrix, which in this case is a correlation matrix. Since the principal diagonal elements are in C(), I don't need P(). Here is the resulting routine:

ReDim CInv(n, n)

For i = 1 To n
CInv(i, i) = 1# / C(i, i)
For j = i + 1 To n
Sum = 0#
For k = i To j - 1
If k = i Then
Sum = Sum - (C(j, k) * a(k, i))
Else
Sum = Sum - (C(j, k) * C(k, i))
End If
Next
CInv(j, i) = Sum / C(j, j)
Next
Next

In comparing the result to that produced by Matlab, the resulting lower triangular matrix has correct main diagonal elements, and element (2,1) is correct, but all other elements are way off. I am posting here the original correlation matrix, the Cholesky dc, the correct inverse of the Choldc (from Matlab), and the erroneous result produced by the above routine.

I would be most grateful if anyone can help me see what I am doing wrong.

Correlation matrix:
1.0000 0.7149 0.6720 0.5339 0.5340 0.4430 0.4060
0.7149 1.0000 0.6100 0.5099 0.4900 0.4100 0.3250
0.6720 0.6100 1.0000 0.2450 0.4820 0.4000 0.2810
0.5339 0.5099 0.2450 1.0000 0.4199 0.1920 0.2170
0.5340 0.4900 0.4820 0.4199 1.0000 0.2420 0.2320
0.4430 0.4100 0.4000 0.1920 0.2420 1.0000 0.1000
0.4060 0.3250 0.2810 0.2170 0.2320 0.1000 1.0000

Cholesky decomposition:
1.0000 0 0 0 0 0 0
0.7149 0.6992 0 0 0 0 0
0.6720 0.1853 0.7170 0 0 0 0
0.5339 0.1834 -0.2061 0.7993 0 0 0
0.5340 0.1548 0.1318 0.1671 0.8035 0 0
0.4430 0.1334 0.1082 -0.0585 -0.0245 0.8776 0
0.4060 0.0497 -0.0014 -0.0115 0.0120 -0.0988 0.9070

Correct inverse of Cholesky decomposition:
1.0000 0 0 0 0 0 0
-1.0224 1.4302 0 0 0 0 0
-0.6730 -0.3697 1.3947 0 0 0 0
-0.6070 -0.4235 0.3597 1.2511 0 0 0
-0.2310 -0.1268 -0.3036 -0.2603 1.2446 0 0
-0.3133 -0.2036 -0.1565 0.0761 0.0348 1.1395 0
-0.4315 -0.1048 -0.0063 0.0276 -0.0126 0.1241 1.1025

Incorrect inverse produced by above routine:
1.0000 0 0 0 0 0 0
-1.0224 1.4302 0 0 0 0 0
-1.1221 -0.2585 1.3947 0 0 0 0
-0.6587 -0.1816 0.2579 1.2511 0 0 0
-1.0237 -0.2612 -0.1211 -0.2080 1.2446 0 0
-0.6459 -0.1584 -0.1334 0.0713 0.0279 1.1395 0
-0.4377 -0.0397 0.0090 0.0041 -0.0159 0.1089 1.1025