/*********************************************************************** * This is a test program for the function cg_method() from the module * * cg.cpp to solve a linear system * * A * X = Y * * for a symmetric positive definite matrix A using the conjugate * * gradient method. * * -------------------------------------------------------------------- * * Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and * * Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * * -------------------------------------------------------------------- * * * SAMPLE RUN: * * * * Solve the following linear system: * * * * 4 size N of matrix * * 5.0 -1.0 -1.0 -1.0 the upper triangle (with diagonal) * * 5.0 -1.0 -1.0 of the positive definite matrix A * * 5.0 -1.0 * * 5.0 * * 1.0 right hand side of the linear system * * 1.0 * * 1.0 * * 1.0 * * * * ----------------------------------------------------- * * Symmetric Linear System (Conjugate Gradient Method) * * ----------------------------------------------------- * * Test data (Matrix and right hand side): * * 5.00000 -1.00000 -1.00000 -1.00000 1.00000 * * -1.00000 5.00000 -1.00000 -1.00000 1.00000 * * -1.00000 -1.00000 5.00000 -1.00000 1.00000 * * -1.00000 -1.00000 -1.00000 5.00000 1.00000 * * * * Solution for the linear system: * * 0.50000 0.50000 0.50000 0.50000 * * ----------------------------------------------------- * * * * C++ version without dynamic allocations by J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***********************************************************************/ #include #include #include "cg.cpp" /* for cg_method() (with header) */ int main() { MAT a; /* the upper triangle of a positive definite real */ /* matrix. */ VEC y, /* right hand side of the linear system */ x; /* solution vector of the system */ int n, /* size of matrix A */ i, /* row index */ j, /* column index */ fehler; /* return value of cg_method() */ /* ----------------------------- print header --------------------- */ printf("-----------------------------------------------------\n"); printf(" Symmetric Linear System (Conjugate Gradient Method) \n"); printf("-----------------------------------------------------\n"); /* ----------------------------- Set inputs ----------------------- */ n=4; /* size of system to solve */ /* define a matrix */ a[0][0]=5.0; a[0][1]=-1.0; a[0][2]=-1.0; a[0][3]=-1.0; a[1][1]= 5.0; a[1][2]=-1.0; a[1][3]=-1.0; a[2][2]= 5.0; a[2][3]=-1.0; a[3][3]= 5.0; /*define b vector */ y[0]=1.0; y[1]=1.0; y[2]=1.0; y[3]=1.0; /* -------------- Print input data as a safeguard ---------------- */ printf(" Test data (Matrix and right hand side):\n"); for (i = 0; i < n; i++) { for (j = 0; j < i; j++) printf(" %9.5f", a[j][i]); for (j = i; j < n; j++) printf(" %9.5f", a[i][j]); printf(" %9.5f\n", y[i]); } /* ---------------- solve linear system --------------------------- */ fehler = cg_method(n, a, y, x); /* perform CG method */ /* -------------------- Print output results ---------------------- */ if (fehler) /* Error in CG method? */ { printf("\n Error in executing CG method.\n\n"); return 4 + fehler; } printf("\n Solution for the linear system:\n"); for (i = 0; i < n; i++) printf(" %9.5f", x[i]); printf("\n"); printf("-----------------------------------------------------\n\n"); return 0; } // end of file cgtst1.cpp