' ******************************************************************
' Cholesky factorization of a positive definite symmetric matrix
' ******************************************************************

#INCLUDE "fbmath.bi"

CONST N = 3

DIM AS DOUBLE A(1 TO N, 1 TO N) = _  ' Original matrix
{{60, 30, 20}, _
 {30, 20, 15}, _
 {20, 15, 12}}

DIM AS DOUBLE  L(1 TO N, 1 TO N)     ' Cholesky factor
DIM AS DOUBLE  B(1 TO N, 1 TO N)     ' Product L * L'
DIM AS INTEGER I, J, K               ' Loop variables

CLS

PRINT "Original matrix:": PRINT

FOR I = 1 TO N
  FOR J = 1 TO N
    PRINT USING "#####.######"; A(I, J);
  NEXT J
  PRINT
NEXT I

' Perform Cholesky factorization
Cholesky A(), L()
IF MathErr <> 0 THEN PRINT "Matrix is not positive definite": STOP

PRINT : PRINT "Cholesky factor (L) :": PRINT

FOR I = 1 TO N
  FOR J = 1 TO I
    PRINT USING "#####.######"; L(I, J);
  NEXT J
  PRINT
NEXT I

' Compute the product L * L' which must
' be equal to the original matrix

FOR I = 1 TO N
  FOR J = 1 TO I
    B(I, J) = 0
    FOR K = 1 TO J
      B(I, J) = B(I, J) + L(I, K) * L(J, K)
    NEXT K
    IF I <> J THEN B(J, I) = B(I, J)
  NEXT J
NEXT I

PRINT : PRINT "Product L * L' :" : PRINT

FOR I = 1 TO N
  FOR J = 1 TO N
    PRINT USING "#####.######"; B(I, J);
  NEXT J
  PRINT
NEXT I

SLEEP
