/***************************************************************** * Inversion of a symmetric matrix by Cholesky 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 * * * * Determinant = 432.000000 * * * * 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 * * * * 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 * * * * C++ Release By Jean-Pierre Moreau, Paris. * * (www.jpmoreau.fr) * * -------------------------------------------------------------- * * Release 1.1 : added verification Inv(A) * A = I. * *****************************************************************/ #include #include #define SIZE 25 typedef double MAT[SIZE][SIZE], VEC[SIZE]; void choldc1(int,MAT,VEC); /* ----------------------------------------------- Cholesky decomposition. input n size of matrix input A Symmetric positive def. matrix output a lower deomposed matrix uses choldc1(int,MAT,VEC) ----------------------------------------------- */ void choldc(int n,MAT A, MAT a) { int i,j; VEC p; for (i = 0; i < n; i++) for (j = 0; j < n; j++) a[i][j] = A[i][j]; choldc1(n, a, p); for (i = 0; i < n; i++) { a[i][i] = p[i]; for (j = i + 1; j < n; j++) { a[i][j] = 0; } } } /* ----------------------------------------------------- Inverse of Cholesky decomposition. input n size of matrix input A Symmetric positive def. matrix output a inverse of lower deomposed matrix uses choldc1(int,MAT,VEC) ----------------------------------------------------- */ void choldcsl(int n, MAT A, MAT a) { int i,j,k; double sum; VEC p; for (i = 0; i < n; i++) for (j = 0; j < n; j++) a[i][j] = A[i][j]; choldc1(n, a, p); for (i = 0; i < n; i++) { a[i][i] = 1 / p[i]; for (j = i + 1; j < n; j++) { sum = 0; for (k = i; k < j; k++) { sum -= a[j][k] * a[k][i]; } a[j][i] = sum / p[j]; } } } /* ----------------------------------------------------------------------------- Computation of Determinant of the matrix using Cholesky decomposition input n size of matrix input a Symmetric positive def. matrix return det(a) uses choldc(int,MAT,MAT) ------------------------------------------------------------------------------ */ double choldet(int n, MAT a) { MAT c; double d=1; int i; choldc(n,a,c); for (i = 0; i < n; i++) d *= c[i][i]; return d * d; } /* --------------------------------------------------- Matrix inverse using Cholesky decomposition input n size of matrix input A Symmetric positive def. matrix output a inverse of A uses choldc1(MAT, VEC) --------------------------------------------------- */ void cholsl(int n, MAT A, MAT a) { int i,j,k; choldcsl(n,A,a); for (i = 0; i < n; i++) { for (j = i + 1; j < n; j++) { a[i][j] = 0.0; } } for (i = 0; i < n; i++) { a[i][i] *= a[i][i]; for (k = i + 1; k < n; k++) { a[i][i] += a[k][i] * a[k][i]; } for (j = i + 1; j < n; j++) { for (k = j; k < n; k++) { a[i][j] += a[k][i] * a[k][j]; } } } for (i = 0; i < n; i++) { for (j = 0; j < i; j++) { a[i][j] = a[j][i]; } } } /* ---------------------------------------------------- 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: ----------------------------------------------------- */ void choldc1(int n, MAT a, VEC p) { int i,j,k; double sum; for (i = 0; i < n; i++) { for (j = i; j < n; j++) { sum = a[i][j]; for (k = i - 1; k >= 0; k--) { sum -= a[i][k] * a[j][k]; } if (i == j) { if (sum <= 0) { printf(" a is not positive definite!\n"); } p[i] = sqrt(sum); } else { a[j][i] = sum / p[i]; } } } } //print a square real matrix A of size n with caption s //(n items per line). void MatPrint(char *s, int n, MAT A) { int i,j; printf("\n %s\n", s); for (i=0; i=0; k--) sum -= A[i][k] * A[j][k]; if (i == j) if (sum <= 0.0) result=0; } } return result; } /****************************************** * MULTIPLICATION OF TWO SQUARE REAL * * MATRICES * * --------------------------------------- * * INPUTS: A MATRIX N*N * * B MATRIX N*N * * N INTEGER * * --------------------------------------- * * OUTPUTS: C MATRIX N*N PRODUCT A*B * * * ******************************************/ void MatMult(int n, MAT A,MAT B, MAT C) { double SUM; int I,J,K; for (I=0; I