'********************************************************************* '* Test program for Gauss Seidel method * '* ----------------------------------------------------------------- * '* SAMPLE RUN: * '* (Solve a linear system by Gauss Seidel iterative method). * '* * '* Input file tseidel2.dat contains: * '* * '* 16 * '* 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 * '* -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 * '* 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 * '* 0 0 -1 4 0 0 0 -1 0 0 0 0 0 0 0 0 * '* -1 0 0 0 4 -1 0 0 -1 0 0 0 0 0 0 0 * '* 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 * '* 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 * '* 0 0 0 -1 0 0 -1 4 0 0 0 -1 0 0 0 0 * '* 0 0 0 0 -1 0 0 0 4 -1 0 0 -1 0 0 0 * '* 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 * '* 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 * '* 0 0 0 0 0 0 0 -1 0 0 -1 4 0 0 0 -1 * '* 0 0 0 0 0 0 0 0 -1 0 0 0 4 -1 0 0 * '* 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 * '* 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 * '* 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 * '* 2 * '* 1 * '* 1 * '* 2 * '* 1 * '* 0 * '* 0 * '* 1 * '* 1 * '* 0 * '* 0 * '* 1 * '* 2 * '* 1 * '* 1 * '* 2 * '* * '* See full results in output file tseidel2.lst. * '* * '* Number of iterations: 30 * '* * '* Solution Vector: * '* 1.000000001332439 * '* .9999999987492376 * '* .9999999990155366 * '* .9999999993882775 * '* .9999999987492376 * '* .9999999968816945 * '* .9999999987457741 * '* .9999999988433159 * '* .9999999990155365 * '* .9999999987457738 * '* .9999999994186047 * '* .9999999996013789 * '* .9999999993882775 * '* .9999999988433159 * '* .9999999996013789 * '* .9999999993760317 * '* * '* ----------------------------------------------------------------- * '* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and * '* Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************************* 'PROGRAM TEST_SEIDEL DEFDBL A-H, O-Z DEFINT I-N DIM rc AS INTEGER 'return code from Gauss Seidel (must be zero) ONE = 1# TWO = 2# ZERO = 0# omega = 1.5 'relaxation factor krit = 0 'no special criterion 'open input file OPEN "tseidel2.dat" FOR INPUT AS #1 INPUT #1, n IF n < 1 THEN PRINT #1, PRINT #1, " Dimension must be > 0" END END IF DIM xmat(n, n), b(n), x(n), residu(n) FOR i = 1 TO n x(i) = ZERO NEXT i FOR i = 1 TO n FOR j = 1 TO n INPUT #1, xmat(i, j) NEXT j NEXT i FOR i = 1 TO n INPUT #1, b(i) NEXT i CLOSE #1 'open output file OPEN "tseidel2.lst" FOR OUTPUT AS #2 PRINT #2, PRINT #2, "-----------------------------------------------------------------------" PRINT #2, " Gauss Seidel method" PRINT #2, "-----------------------------------------------------------------------" PRINT #2, " Dimension of the input xmatrix = "; n PRINT #2, PRINT #2, " Input xmatrix:" FOR i = 1 TO n FOR j = 1 TO n - 1 PRINT #2, xmat(i, j); NEXT j PRINT #2, xmat(i, n) NEXT i PRINT #2, PRINT #2, " Second member:" FOR i = 1 TO n - 1 PRINT #2, b(i); NEXT i PRINT #2, b(n) PRINT #2, PRINT #2, " Solution vector:" GOSUB 1000 'call seidel(krit, n, xmat, b, omega, x, residu, iter, rc) IF rc <> 0 THEN PRINT #2, " Error in seidel: rc = "; rc END END IF FOR i = 1 TO n PRINT #2, x(i) NEXT i PRINT #2, PRINT #2, " Residual vector (must be near zero):" FOR i = 1 TO n PRINT #2, residu(i) NEXT i PRINT #2, PRINT #2, " Number of iterations: "; iter PRINT #2, "-----------------------------------------------------------------------" CLS PRINT PRINT " Results in tseidel2.lst..." PRINT INPUT "", R$ CLOSE #2 END '*====================================================================* '* * '* seidel solves the linear system xmat * x = b iteratively. * '* Here xmat is a nonsingular n x n xmatrix, 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 xmat, * '* convergence of the scheme is guaranteed. * '* * '*====================================================================* '* * '* Input parameters: * '* ================ * '* krit integer * '* select criterion * '* =1 : row sum criterion * '* =2 : column sum criterion * '* =3 : criterion of Schmidt-v.Mises * '* other : no check * '* n integer ( n > 0 ) * '* size of xmat, b and x * '* xmat Double xmat(n,n) * '* Matrix of the liear system * '* b Double b(n) * '* Right hand side * '* omega Double omega ( 0.0 < omega < 2.0 ) * '* Relaxation coefficient. * '* x Double x(n) * '* Starting vector for iteration * '* * '* Output parameters: * '* ================== * '* x Double x(n) * '* solution vector * '* residu Double residu(n) * '* residual vector b - xmat * x; close to zero vector * '* iter integer * '* Number of iterations performed * '* * '* Return code rc: * '* ============== * '* = 0 solution has been found * '* = 1 n < 1 or omega <= 0 or omega >= 2 * '* = 2 improper xmat or b or x (not used here) * '* = 3 one diagonal element of xmat vanishes * '* = 4 Iteration number exceeded * '* = 11 column sum criterion violated * '* = 12 row sum criterion violated * '* = 13 Schmidt-v.Mises criterion violated * '* * '*====================================================================* 1000 'Subroutine Seidel 'Labels; 10, 20 ITERMAX = 300 'Maximum number of iterations rc = 0 'Initialize return code iter = 0 'Initialize iteration counter IF n < 1 OR omega <= ZERO OR omega >= TWO THEN 'Check omega rc = 1 RETURN END IF eps = 1E-08 'required precision FOR i = 1 TO n 'transform xmat so that all IF xmat(i, i) = ZERO THEN 'diagonals equal 1. rc = 3 RETURN END IF tmp = ONE / xmat(i, i) FOR j = 1 TO n xmat(i, j) = xmat(i, j) * tmp NEXT j b(i) = b(i) * tmp 'adjust right hand side b NEXT i 'check convergence criteria IF krit = 1 THEN FOR i = 1 TO n 'row sum criterion tmp = ZERO FOR j = 1 TO n tmp = tmp + ABS(xmat(i, j)) IF tmp >= TWO THEN rc = 11 RETURN END IF NEXT j NEXT i ELSEIF krit = 2 THEN FOR j = 1 TO n 'column sum criterion tmp = ZERO FOR i = 1 TO n tmp = tmp + ABS(xmat(i, j)) IF tmp >= TWO THEN rc = 12 RETURN END IF NEXT i NEXT j ELSEIF krit = 3 THEN tmp = ZERO FOR i = 1 TO n FOR j = 1 TO n 'criterion of Schmidt, tmp = tmp + xmat(i, j) ^ 2 'Von Mises. NEXT j tmp = SQR(tmp - ONE) IF tmp >= ONE THEN rc = 13 RETURN END IF NEXT i END IF 'end check criteria FOR i = 1 TO n residu(i) = x(i) 'store x in residu NEXT i WHILE iter <= ITERMAX 'Begin iteration iter = iter + 1 FOR i = 1 TO n tmp = b(i) FOR j = 1 TO n tmp = tmp - xmat(i, j) * residu(j) NEXT j residu(i) = residu(i) + omega * tmp NEXT i FOR i = 1 TO n 'check break-off criterion tmp = x(i) - residu(i) IF ABS(tmp) <= eps THEN x(i) = residu(i) 'if rc = 0 at end of loop rc = 0 ' -> stop iteration. ELSE FOR j = 1 TO n x(j) = residu(j) NEXT j rc = 4 GOTO 10 'continue iteration END IF NEXT i IF rc = 0 THEN GOTO 20 'solution found 'End iteration 10 WEND 20 FOR i = 1 TO n 'find residual vector tmp = b(i) FOR j = 1 TO n tmp = tmp - xmat(i, j) * x(j) NEXT j residu(i) = tmp NEXT i RETURN ' J-P Moreau April, 2005