'************************************************************** '* Program to test Cholesky Method for symmetric * '* and positive definite matrices * '* * '* Basic version by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ---------------------------------------------------------- * '* Reference: * '* * '* "Numerical Algorithms with C, By Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* ---------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Input file: matsym.dat * '* * '* 4 * '* 5.0 -1.0 -1.0 -1.0 1.0 * '* 5.0 -1.0 -1.0 1.0 * '* 5.0 -1.0 1.0 * '* 5.0 1.0 * '* * '* Output file: matsym.lst * '* * '* ---------------------------------------------------------- * '* Solving a linear system by Cholesky Method * '* ---------------------------------------------------------- * '* The matrix of left hand coefficients must be a symmetric * '* positive definite matrix (else error code = 2) * '* * '* Size of system: 4 * '* * '* Input Linear System: * '* * '* 5.000000 -1.000000 -1.000000 -1.000000 1.000000 * '* * '* -1.000000 5.000000 -1.000000 -1.000000 1.000000 * '* * '* -1.000000 -1.000000 5.000000 -1.000000 1.000000 * '* * '* -1.000000 -1.000000 -1.000000 5.000000 1.000000 * '* * '* Solution: * '* * '* 0.500000 0.500000 0.500000 0.500000 * '* ---------------------------------------------------------- * '* * '************************************************************** defint i-n defdbl a-h,o-z ' Only diagonal and upper-diagonal elements are read ' main variables ' a, org : matrix (n,n) ' b, x : vector (n) ' n, i, j, mode, irc : INTEGER ' finput\$, foutput\$, s\$ : STRINGS 'begin main program EPSQUAD = 1e-12# ZERO = 0# cls print input " Input data file name (without *.dat): ", s\$ finput\$ = s\$ + ".dat" foutput\$ = s\$ + ".lst" ' open input, output files open finput\$ for input as #1 open foutput\$ for output as #2 ' print header to output file print #2,"----------------------------------------------------------" print #2," Solving a linear system by Cholesky Method " print #2,"----------------------------------------------------------" print #2," The matrix of left hand coefficients must be a symmetric " print #2," positive definite matrix (else error code = 2) " print #2, ' read size n of system from input file input #1, n if n < 1 then print " n must be > 0 !" irc = 1 end end if dim a(n,n), org(n,n), b(n), x(n) for i = 0 to n-1 'read left hand coefficients of line i for j = i to n-1 input #1, tmp org(i,j) = tmp org(j,i) = tmp 'symmetric term next j 'read right hand coefficient of line i input #1, b(i) next i close #1 'print data to output file print #2," Size of system: "; n print #2, print #2," Input Linear System:" print #2, f\$="#####.######" for i = 0 to n-1 for j = 0 to n-1 a(i,j) = org(i,j) print #2,using f\$; a(i,j); next j print #2,using f\$; b(i) print #2, next i mode = 0 'Factor matrix a and solve linear system 'call Cholesky procedure gosub 1000 'print results to output file if irc = 0 then print #2," Solution:" print #2, for i=0 to n-1 print #2,using f\$; x(i); next i else print #2," ERROR choly: code = "; irc end if print #2, print #2,"----------------------------------------------------------" close #2 print print " Results in file "; foutput\$ END 'of main program 1000 'Subroutine Choly '====================================================================== '* * '* The cholesky subroutine solves linear systems : a * x = b * '* for positive definite symmetric n x n matrices mat using the * '* Cholesky method. * '* * '* a must be symmetric and positive definite, or the method will * '* fail. i.e. for all nonzero vectors y we must have y' * a * y > 0 * '* * '* cholesky uses only the lower triangle of matrix a. * '* * '*====================================================================* '* * '* Application: * '* ============ * '* * '* Solve linear systems with positive definite system matrices * '* efficeiently. * '* * '*====================================================================* '* * '* Control parameter: * '* ================== * '* mode integer * '* = 0 Factor matrix a and solve linear system * '* = 1 Compute factorization only * '* = 2 Solve system only; for this call the factorization * '* must be stored in matrix a from chodec. Used to * '* solve for many right hand sides. * '* * '* Input parameters: * '* ================ * '* n integer (n > 0) * '* Dimension of a, size of b and x * '* a REAL matrix (n,n) * '* mode = 0, 1: Matrix of the linear system * '* mode = 2 : Cholesky factor * '* b REAL vector (n) (for mode = 0, 2) * '* Right hand side of system of equations * '* * '* Output parameters: * '* ================== * '* a REAL matrix (n,n) (for mode = 0, 1) * '* Cholesky decomposition of input matrix mat * '* x REAL vector (n) (for mode = 0, 2) * '* solution vector * '* * '* Return value rc: * '* =============== * '* = 0 all ok * '* = 1 n < 1 or false input parameter * '* = 2 Matrix is not positive definite * '* = 3 wrong value for mode * '* * '*====================================================================* '* * '* Subroutines in use: * '* ================== * '* 2000: chodec: Factorization * '* 3000: chosol: Solving system * '* * '====================================================================== irc=3 if n < 1 then irc=1 return end if if mode=0 then 'Factor matrix and solve system 'chodec (n, a, irc) gosub 2000 if irc = 0 then 'chosol(n, a, b, x, irc) gosub 3000 end if end if if mode=1 then 'Factor only 'chodec (n, a, irc) gosub 2000 end if if mode=2 then 'Ssolve only 'chosol (n, a, b, x, irc) gosub 3000 end if Return 2000 'Procedure chodec '====================================================================== '* * '* chodec decomposes the symmetric positive definite matrix mat. * '* * '*====================================================================* '* * '* Input parameters: * '* ================ * '* n integer (n > 0) * '* Dimension of matrix a * '* a REAL matrix (n,n) * '* Matrix of left hand coefficients * '* * '* Output parameter: * '* ================ * '* a REAL matrix (n,n) * '* Cholesky decomposition in the lower triangle * '* * '* Return value rc: * '* =============== * '* = 0 all ok * '* = 1 n < 1 * '* = 2 Matrix not positive definite * '* * '*====================================================================* '* * '* Functions in use: SQR (square root in Basic) * '* ================ * '* * '*====================================================================* '* * '* Constants in use: EPSQUAD (small number) * '* ================ * '* * '====================================================================== if n < 1 then irc=1 ' n < 1 error ! return end if if a(0,0) < EPSQUAD then ' matrix a not positive definite irc=2 return end if a(0,0) = SQR(a(0,0)) for j = 1 to n-1 a(j,0) = a(j,0) / a(0,0) next j for i = 1 to n-1 sum = a(i,i) for j = 0 to i-1 sum = sum - a(i,j)^2 next j if sum < EPSQUAD then ' matrix a not positive definite irc=2 return end if a(i,i) = SQR(sum) for j = i + 1 to n-1 sum = a(j,i) for k = 0 to i-1 sum = sum - a(i,k) * a(j,k) next k a(j,i) = sum / a(i,i) next j next i irc = 0 return 3000 'Procedure chosol '====================================================================== '* * '* chosol finds the solution x of the linear system B' * B * x = b * '* for a lower triangular nonsingular matrix B as supplied in chodec.* '* * '*====================================================================* '* * '* Input parameters: * '* ================ * '* n integer (n > 0) * '* Dimension of matrix a, size of vectors b and x * '* a REAL matrix (n,n) * '* lower triangular matrix as supplied by chodec * '* b REAL vector b(n) * '* Right hand side * '* * '* Output parameter: * '* ================= * '* x REAL vector x(n) * '* solution vector * '* * '* Return value rc: * '* =============== * '* = 0 all ok * '* = 1 improper lwer triangular matrix or n < 1 * '* * '====================================================================== if n < 1 then ' n < 1 error ! irc=1 return end if if a(0,0) = ZERO then ' improper factor matrix irc=1 return end if x(0) = b(0) / a(0,0) ' update right hand side for k = 1 to n-1 sum = ZERO for j = 0 to k-1 sum = sum + a(k,j) * x(j) next j if a(k,k) = ZERO then irc=1 return end if x(k) = (b(k) - sum) / a(k,k) next k x(n-1) = x(n-1) / a(n-1,n-1) ' back substitution for k=n-2 to 0 step -1 sum = ZERO for j = k + 1 to n-1 sum = sum + a(j,k) * x(j) next j x(k) = (x(k) - sum) / a(k,k) next k irc = 0 Return 'end of file tcholy1.bas