/***************************************************************** * Solving a symmetric linear system by Gauss method * * * * C++ version by J-P Moreau * * (www.jpmoreau.fr) * * -------------------------------------------------------------- * * Reference: * * * * "ALGEBRE Algorithmes et programmes en Pascal * * de Jean-Louis Jardrin - Dunod BO-PRE 1988" [BIBLI 10]. * * -------------------------------------------------------------- * * SAMPLE RUN: * * * * Input file (matsym.dat): * * * * 4 * * 8 2 3 12 25 * * 4 7 0.25 13.25 * * 3 5 18 * * 2 19.25 * * * * Output file (matsym.lst): * * * * -------------------------------------------------------------- * * SYMMETRIC SYSTEM TO BE SOLVED: * * -------------------------------------------------------------- * * * * N=4 * * * * 8.000000 2.000000 3.000000 12.000000 25.000000 * * 2.000000 4.000000 7.000000 0.250000 13.250000 * * 3.000000 7.000000 3.000000 5.000000 18.000000 * * 12.000000 0.250000 5.000000 2.000000 19.250000 * * * * -------------------------------------------------------------- * * System solution: * * -------------------------------------------------------------- * * * * X1 = 1.000000 * * X2 = 1.000000 * * X3 = 1.000000 * * X4 = 1.000000 * * * * -------------------------------------------------------------- * * End of file matsym.lst. * * * *****************************************************************/ #include #include /******************************************************************* * This function solves a symmetric linear system by Gauss method * * ---------------------------------------------------------------- * * The input matrix includes the right hand column, only the upper * * triangle terms are given. * * * * INPUTS: * * eps: desired precision (real) * * n: size of system (integer) * * A: extended matrix of system of type MAT * * OUTPUTS: * * it: error indicator (0=singular system, 1=OK) * * X: system solution (x1,...,xn) of type VEC * * ---------------------------------------------------------------- * * LIMITATION: No term in main diagonal must equal zero. * *******************************************************************/ void RSLSG (REAL eps, int n, REAL *A[], int *it, REAL X[]) { REAL q0, s; int i,j,k; *it=1; k=1; while ((*it != 0) && (k < n)) { if (fabs(A[k][k]) < eps) *it=0; else { for (i=k+1; i=eps)) { X[n] = A[n][n+1]/A[n][n]; for (i=n-1; i>0; i--) { s=0.0; for (j=i+1; j