'****************************************************************** '* Inversion of a symmetric matrix by Choleski decomposition. * '* The matrix must be positive definite. * '* -------------------------------------------------------------- * '* REFERENCE: * '* From a Java Library Created by Vadim Kutsyy, * '* "http://www.kutsyy.com". * '* -------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Inversion of a square real symetric matrix by Cholevsky method * '* (The matrix must positive definite). * '* * '* Size = 4 * '* * '* Matrix A: * '* 5.000000 -1.000000 -1.000000 -1.000000 * '* -1.000000 5.000000 -1.000000 -1.000000 * '* -1.000000 -1.000000 5.000000 -1.000000 * '* -1.000000 -1.000000 -1.000000 5.000000 * '* * '* Determinant = 432.000000 * '* * '* Matrix Inv(A): * '* 0.250000 0.083333 0.083333 0.083333 * '* 0.083333 0.250000 0.083333 0.083333 * '* 0.083333 0.083333 0.250000 0.083333 * '* 0.083333 0.083333 0.083333 0.250000 * '* * '* Basic Release By Jean-Pierre Moreau, Paris. * '* (www.jpmoreau.fr) * '* -------------------------------------------------------------- * '* Release 1.1 : added verification Inv(A) * A = I. * '****************************************************************** ' List of subroutines: ' ------------------- ' 500: Print matrix A(n,n) with caption ' 501: Print matrix aa(n,n) with caption ' 502: Print matrix C(n,n) with caption ' 550: Copy matrix A in matrix A1 ' 600: Calculate determinant ' 700: Check if matrix is positive definite ' 800: Matrix inverse using Cholesky decomposition ' 900: Matrix multiplication C=A1*aa ' 1000: Cholesky decomposition ' 1200: Inverse of Cholesky decomposition ' 1500: Main method for Cholesky decomposition ' ------------------------------------------------- 'Program Cholesky DEFDBL A-H, O-Z DEFINT I-N CLS 'clear screen n = 4 DIM A(n, n), A1(n, n), aa(n, n), C(n, n), p(n) PRINT " Inversion of a square real symetric matrix by Cholesky method" PRINT " (The matrix must positive def.)." PRINT PRINT " Size = "; n ' define lower half of matrix A(0, 0) = 5 A(1, 0) = -1: A(1, 1) = 5 A(2, 0) = -1: A(2, 1) = -1: A(2, 2) = 5 A(3, 0) = -1: A(3, 1) = -1: A(3, 2) = -1: A(3, 3) = 5 ' define upper half by symmetry FOR i = 0 TO n - 1 FOR J = i + 1 TO n - 1 A(i, J) = A(J, i) NEXT J NEXT i s$ = "Matrix A:" GOSUB 500 'call MatPrint(s$,n,A) PRINT GOSUB 700 'check matrix A IF icheck = 1 THEN GOSUB 550 'A1=A GOSUB 600 'd=choldet(n,A) PRINT USING " Determinant = ####.######"; d GOSUB 800 'CALL cholsl(n, A, aa) s$ = "Matrix Inv(A):" GOSUB 501 'call MatPrint(s$,n,aa) ELSE PRINT " This matrix is not positive definite !" END IF PRINT INPUT " Do you want a verification (y/n) ?: ", answer$ IF answer$ = "y" THEN GOSUB 900 'call MatMult(n,A1,aa,C); s$ = "Verification A * Inv(A) = I:" GOSUB 502 'call MatPrint(s$,n,C) END IF END 'of main program ' ----------------------------------------------------- ' print a square real matrix A of size n with caption s ' (n items per line). ' input s caption string ' input n size of matrix ' input A matrix(n,n) ' ----------------------------------------------------- 500 : PRINT PRINT " "; s$ FOR i = 0 TO n - 1 FOR J = 0 TO n - 1 PRINT USING " ###.######"; A(i, J); NEXT J PRINT NEXT i RETURN ' ------------------------------------------------------ ' print a square real matrix aa of size n with caption s ' (n items per line). ' input s caption string ' input n size of matrix ' input aa matrix(n,n) ' ------------------------------------------------------ 501 : PRINT PRINT " "; s$ FOR i = 0 TO n - 1 FOR J = 0 TO n - 1 PRINT USING " ###.######"; aa(i, J); NEXT J PRINT NEXT i RETURN ' ----------------------------------------------------- ' print a square real matrix C of size n with caption s ' (n items per line). ' input s caption string ' input n size of matrix ' input C matrix(n,n) ' ----------------------------------------------------- 502 : PRINT PRINT " "; s$ FOR i = 0 TO n - 1 FOR J = 0 TO n - 1 PRINT USING " ###.######"; C(i, J); NEXT J PRINT NEXT i RETURN 'copy matrix A in matrix A1 550 FOR i = 0 TO n - 1 FOR J = 0 TO n - 1 A1(i, J) = A(i, J) NEXT J NEXT i RETURN ' ---------------------------------------------------------------------- ' Computation of Determinant of the matrix using Cholesky decomposition ' input n size of matrix ' input A Symmetric positive def. matrix ' output d det(A) ' uses choldc(int,MAT,MAT) ' ---------------------------------------------------------------------- 'choldet 600 : d = 1# GOSUB 1000 'call choldc(n,A,aa) FOR i = 0 TO n - 1 d = d * aa(i, i) NEXT i d = d * d RETURN 'check if matrix A is positive definite (icheck=1) 700 : icheck = 1 FOR i = 0 TO n - 1 FOR J = i TO n - 1 sum = A(i, J) FOR k = i - 1 TO 0 STEP -1 sum = sum - A(i, k) * A(J, k) NEXT k IF i = J THEN IF sum <= 0 THEN icheck = 0 END IF END IF NEXT J NEXT i RETURN ' --------------------------------------------------- ' Matrix inverse using Cholevsky decomposition ' input n size of matrix ' input A Symmetric psitive def. matrix ' output aa inverse of A ' uses choldc1(int,MAT,VEC) ' --------------------------------------------------- 'cholsl 800 : GOSUB 1200 'call choldcsl(n,A,aa) FOR i = 0 TO n - 1 FOR J = i + 1 TO n - 1 aa(i, J) = 0# NEXT J NEXT i FOR i = 0 TO n - 1 aa(i, i) = aa(i, i) * aa(i, i) FOR k = i + 1 TO n - 1 aa(i, i) = aa(i, i) + aa(k, i) * aa(k, i) NEXT k FOR J = i + 1 TO n - 1 FOR k = J TO n - 1 aa(i, J) = aa(i, J) + aa(k, i) * aa(k, J) NEXT k NEXT J NEXT i FOR i = 0 TO n - 1 FOR J = 0 TO i - 1 aa(i, J) = aa(J, i) NEXT J NEXT i RETURN '******************************************* '* MULTIPLICATION OF TWO SQUARE REAL * '* MATRICES * '* --------------------------------------- * '* INPUTS: A1 MATRIX N*N * '* aa MATRIX N*N * '* n INTEGER * '* --------------------------------------- * '* OUTPUTS: C MATRIX N*N PRODUCT A1*aa * '******************************************* 900 : FOR i = 0 TO n - 1 FOR J = 0 TO n - 1 sum = 0# FOR k = 0 TO n - 1 sum = sum + A1(i, k) * aa(k, J) NEXT k C(i, J) = sum NEXT J NEXT i RETURN ' ----------------------------------------------- ' Cholesky decomposition. ' input n size of matrix ' input A Symmetric positive def. matrix ' output aa lower deomposed matrix ' uses choldc1(int,pMAT,pVEC) ' ----------------------------------------------- ' choldc 1000 : FOR i = 0 TO n - 1 FOR J = 0 TO n - 1 aa(i, J) = A(i, J) NEXT J NEXT i GOSUB 1500 'call choldc1(n, aa, p) FOR i = 0 TO n - 1 aa(i, i) = p(i) FOR J = i + 1 TO n - 1 aa(i, J) = 0# NEXT J NEXT i RETURN ' --------------------------------------------------- ' Inverse of Cholesky decomposition. ' input n size of matrix ' input A Symmetric positive def. matrix ' output aa inverse of lower decomposed matrix ' uses choldc1(int,MAT,VEC) ' --------------------------------------------------- 'choldcsl 1200 : FOR i = 0 TO n - 1 FOR J = 0 TO n - 1 aa(i, J) = A(i, J) NEXT J NEXT i GOSUB 1500 'call choldc1(n, aa, p) FOR i = 0 TO n - 1 aa(i, i) = 1# / p(i) FOR J = i + 1 TO n - 1 sum = 0# FOR k = i TO J - 1 sum = sum - aa(J, k) * aa(k, i) NEXT k aa(J, i) = sum / p(J) NEXT J NEXT i RETURN ' ---------------------------------------------------- ' main method for Cholesky decomposition. ' input n size of matrix ' input/output A Symmetric positive def. matrix ' output p vector of resulting diag of a ' author: ' ----------------------------------------------------- 'choldc1 1500 : FOR i = 0 TO n - 1 FOR J = i TO n - 1 sum = A(i, J) FOR k = i - 1 TO 0 STEP -1 sum = sum - A(i, k) * A(J, k) NEXT k IF i = J THEN IF sum <= 0 THEN PRINT " the matrix is not positive definite!" END IF p(i) = SQR(sum) ELSE A(J, i) = sum / p(i) END IF NEXT J NEXT i RETURN 'end of file cholev.bas