/************************************************************* * Program to test Cholesky Method for symmetric * * and positive definite matrices * * * * C++ 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 * * ---------------------------------------------------------- * * * *************************************************************/ // Program to link with basis_r.cpp,vmblock.cpp, fcholy.cpp #include #include #include "Fcholy.h" // Only diagonal and upper-diagonal elements are read int main () { REAL **a, *b, *x, **org; REAL tmp; int n, i, j, cas, rc; void *vmblock; FILE *fp1, *fp2; char input[20], output[20], s[20]; printf("\n Input data file name (without *.dat): "); scanf("%s",s); strcpy(input,s); strcat(input,".dat"); strcpy(output,s); strcat(output,".lst"); fp1=fopen(input,"r"); fp2=fopen(output,"w"); WriteHead(fp2," Solving a linear system by Cholesky Method"); fprintf(fp2," The matrix of left hand coefficients must be a symmetric\n"); fprintf(fp2," positive definite matrix (else error code = 2)\n\n"); if (fscanf (fp1,"%d", &n) <= 0) { LogError ("Input Stream", 0, __FILE__, __LINE__); return 1; } if (n < 1) { LogError ("n must be > 0", 0, __FILE__, __LINE__); return 1; } //dynamic allocations vmblock = vminit(); a = (REAL **)vmalloc(vmblock, MATRIX, n, n); org = (REAL **)vmalloc(vmblock, MATRIX, n, n); x = (REAL *) vmalloc(vmblock, VEKTOR, n, 0); b = (REAL *) vmalloc(vmblock, VEKTOR, n, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 1; } for (i = 0; i < n; i++) { //read left hand coefficients of line i for (j = i; j < n; j++) { if (fscanf (fp1,FORMAT_IN, &tmp) <= 0) { LogError ("Input Stream", 0, __FILE__, __LINE__); return 1; } org[i][j] = tmp; org[j][i] = tmp; } //read right hand coefficient of line i fscanf(fp1,FORMAT_IN,&tmp); b[i]=tmp; } fclose(fp1); fprintf(fp2," Size of system: %d\n\n",n); fprintf(fp2," Input Linear System:\n\n"); for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { a[i][j] = org[i][j]; fprintf(fp2,FORMAT_126LF, a[i][j]); } fprintf(fp2,FORMAT_126LF, b[i]); fprintf(fp2,"\n"); } cas = 0; //call Cholesky subroutine rc = choly ( cas, n, a, b, x ); if (rc == 0) { fprintf(fp2,"\n Solution:\n\n"); WriteVec(fp2,n,x); } else { fprintf(fp2,"\n ERROR choly: code=%d\n",rc); } WriteEnd(fp2); fclose(fp2); printf("\n Results in %s.\n\n",output); return (0); } // end of file tcholy.cpp