/* ---------------------- MODULE fseidel.cpp ------------------------ */ #include "basis.h" #define ITERMAX 300 /* Maximal number of iterations */ int seidel /* Gauss Seidel Method with relaxation .......*/ ( int crit, /* crit = 0, 1, 2, 3 ...............*/ int n, /* size of matrix ..................*/ REAL * mat[], /* matrix ..........................*/ REAL b[], /* Right hand side .................*/ REAL omega, /* Relaxaktion coefficient .........*/ REAL x[], /* solution ........................*/ REAL residu[], /* Residuum vector .................*/ int * iter /* # of iterations .................*/ ) /*====================================================================* * * * seidel solves the linear system mat * x = b iteratively. * * Here mat is a nonsingular n x n matrix, b is the right hand * * side for the linear system and x is the solution. * * * * seidel uses the Gauss Seidel Method with relaxation for a given * * relaxation coefficient 0 < omega < 2. * * If omega = 1, the standard Gauss Seidel method (without * * relaxation) is performed. * * * *====================================================================* * * * Applications: * * ============= * * Solve linear systems with nonsingular system matrices that * * satisfy one of the following criteria: row sum criterion, * * column sum criterion or the criterion of Schmidt and v. Mises.* * Only if at least one of these criteria is satisfied for mat, * * convergence of the scheme is guaranteed. * * * *====================================================================* * * * Input parameters: * * ================ * * crit int crit; * * select criterion * * =1 : row sum criterion * * =2 : column sum criterion * * =3 : criterion of Schmidt-v.Mises * * other : no check * * n int n; ( n > 0 ) * * size of mat, b and x * * mat REAL *mat[n]; * * Matrix of the liear system * * b REAL b[n]; * * Right hand side * * omega REAL omega; ( 0.0 < omega < 2.0 ) * * Relaxation coefficient. * * x REAL x[n]; * * Starting vector for iteration * * * * Output parameters: * * ================== * * x REAL x[n]; * * solution vector * * residu REAL residu[n]; * * residual vector b - mat * x; close to zero vector * * iter int *iter; * * Number of iterations performed * * * * Return value : * * ============= * * = 0 solution has been found * * = 1 n < 1 or omega <= 0 or omega >= 2 * * = 2 improper mat or b or x * * = 3 one diagonal element of mat vanishes * * = 4 Iteration number exceeded * * = 11 column sum criterion violated * * = 12 row sum criterion violated * * = 13 Schmidt-v.Mises criterion violated * * * *====================================================================* * * * Constants used : NULL, MACH_EPS * * ================ * * * * Macros : ABS, SQR, SQRT * * * * Reference: [BIBLI 11] * *====================================================================*/ { int i, j, rc = 0; REAL tmp, eps; *iter = 0; /* Initialize iteration counter */ if ( n < 1 || /* Check omega */ omega <= ZERO || omega >= TWO ) return (1); if (mat == NULL) return (2); for (i = 0; i < n; i++) if (mat[i] == NULL) return (2); if (x == NULL || b == NULL || residu == NULL) return (2); eps = (REAL) (MACH_EPS * (REAL) 128.0); for (i = 0; i < n; i++) /* transform mat so that all*/ { /* diagonals equal 1 */ if (mat[i][i] == ZERO) return (3); tmp = ONE / mat[i][i]; for (j = 0; j < n; j++) mat[i][j] *= tmp; b[i] *= tmp; /* adjust right hand side b */ } switch (crit) /* check convergence criteria */ { case 1: for (i = 0; i < n; i++) /* row sum criterion */ { for (tmp = ZERO, j = 0; j < n; j++) tmp += ABS (mat[i][j]); if (tmp >= TWO) return (11); } break; case 2: for (j = 0; j < n; j++) /* column sum criterion */ { for (tmp = ZERO, i = 0; i < n; i++) tmp += ABS (mat[i][j]); if (tmp >= TWO) return (12); } break; case 3: for (tmp = ZERO, i = 0; i < n; i++) for (j = 0; j < n; j++) /* criterion of Schmidt, */ tmp += SQR (mat[i][j]); /* v. Mises */ tmp = SQRT (tmp - ONE); if (tmp >= ONE) return (13); break; default: break; /* do not check */ } for (i = 0; i < n; i++) residu[i] = x[i]; /* store x in residu */ while (*iter <= ITERMAX) /* Begin iteration */ { (*iter)++; for (i = 0; i < n; i++) { for (tmp = b[i], j = 0; j < n; j++) tmp -= mat[i][j] * residu[j]; residu[i] += omega * tmp; } for (i = 0; i < n; i++) /* check break-off criterion */ { tmp = x[i] - residu[i]; if (ABS (tmp) <= eps) { x[i] = residu[i]; /* If rc = 0 at end of loop */ rc = 0; /* -> stop iteration */ } else { for (j = 0; j < n; j++) x[j] = residu[j]; rc = 4; break; } } if (rc == 0) break; /* solution found */ } /* End iteration */ for (i = 0; i < n; i++) /* find residual vector */ { for (tmp = b[i], j = 0; j < n; j++) tmp -= mat[i][j] * x[j]; residu[i] = tmp; } return (rc); } /* ----------------------- END fseidel.cpp -------------------------- */