C[BA*) C[KA{F 5} C[ {Iterative Methods for Linear Systems} C[ {Iterative Methods for Linear Systems}*) C[FE{F 5.4} C[ {The Gau"s-Seidel Iteration} C[ {The Gau"s-Seidel Iteration}*) C[LE*) c SUBROUTINE ADSOR(A,N,IA,B,X,KADAPT,EPS,KMAX,IMETH,ISWITC, C[IX{ADSOR}*) c * OMEGA,WORK,RES,ITNUMB,IERR) SUBROUTINE ADSOR(A,N,IA,B,X,IERR) C C***************************************************************** C * C This program solves an inhomogeneous linear system AX = B of * C equations with a nonsingular system matrix A. The method of * C Jacobi is used jointly with relaxation, where the relaxation * C parameter OMEGA is adjusted during the iteration (adaptive * C SOR method). * C[BE*) C For a suitable choice of parameters (refer to the remark * C below), this program can perform the Gauá-Seidel method or * C a non-adaptive SOR method. * C * C * C INPUT PARAMETERS: * C ================= * C A : 2-dimensional array A(1:IA,1:N), containing the * C system matrix for the linear equations * C N : size of the linear system * C IA : leading dimension of A, as specified in the calling * C program * C B : N-vector B(1:N), the right hand side of the system * C X : N-vector X(1:N) containing the starting value for * C iteration * C KADAPT : Number of iterations, after which the relaxation * C parameter is to be redefined * C EPS : desired accuracy; the iteration is stopped when the * C maximum norm of the relative error does not exceed * C EPS * C KMAX : Maximal number of iterations allowed * C IMETH : parameter that determines the method used: * C = 0, adaptive SOR method * C = 1, SOR method for a given relaxation parameter * C = 2, Gauá-Seidel method * C ISWITC : parameter that determines the convergence criterion * C to be used: * C = 0, none * C = 1, row sum criterion * C = 2, column sum criterion * C = 3, criterion of Schmidt and v. Mises * C OMEGA : in case IMETH=1, the optimal relaxation parameter * C must be part of the input; otherwise only the name * C must be declared in the callimng program. * C * C * C REMARKS: * C ======== * C For the adaptive SOR method (IMETH=0) we recommend to set * C KADAPT=4 or KADAPT=5. * C If the optimal relaxationcoefficient Wopt is known for A, then* C one should set IMETH=1 and OMEGA = Wopt, i.e., the SOR method * C with given optimal relaxation coefficient should be used. * C If IMETH=2, then the program performs the Gauá-Seidel method. * C * C * C AUXILIARY PARAMETERS: * C ===================== * C WORK : 2-dim. array WORK(1:N,1:3) * C * C * C OUTPUT PARAMETERS: * C ================== * C A : 2-dim. array A(1:IA,1:N), the input matrix A is over-* C written by: A(I,J)=A(I,J)/A(I,I) for I,J=1, ..., N * C B : N-vector B(1:N), the right hand side is replaced by * C B(I)=B(I)/A(I,I); I=1,N * C OMEGA : - if IMETH = 0, the program returns the adaptively * C computed relaxations parameter. * C - if IMETH = 1, the optimal relaxation parameter * C is returned as put in externally. * C - if IMETH = 2, then on output OMEGA = 1. * C X : N-vector X(1:N) that contains the solution vector * C RES : N-vector RES(1:N) containing the residuum B - AX; * C the residuum is available even if the desired * C accuracy EPS could not be achieved with the given * C maximum number of iterations. * C ITNUMB : num,bert of iterations actually performed * C IERR : error parameter: * C = 0, the desired convergence criterium has not been * C met * C = 1, the solution X has been found * C = 2, the desired accuracy has not been achieved after* C KMAX iterations * C = 3, input data incorrect * C = 4, system matrix A is numerically singular * C * C----------------------------------------------------------------* C * C Required subroutines: GAUSEI, MNORM, CONV, RESID, MACHPD * C * C***************************************************************** C * C Author : Gisela Engeln-Mllges * C Date : 06.09.1992 * C Source : FORTRAN 77 * C * C[BA*) C***************************************************************** C[BE*) C C Declarations C DOUBLE PRECISION A(1:IA,1:N),B(1:N),X(1:N),WORK(1:N,1:3), * RES(1:N),EPS,OMEGA,FMACHP,HELP,DIFFN,Q, * RELERR,SUM,XN C c The following 5 lines is added by GU EPS=1.0D-06 KADAPT=4 KMAX=2000 IMETH=2 ISWITC=0 OMEGA=1.0d0 c C Checking the inputs EPS, KMAX, IMETH and ISWITC C IF(EPS .LE. 0.0D0 .OR. KMAX .LT. 1 .OR. ISWITC .LT. 0 .OR. * ISWITC .GT. 3 .OR. IMETH .LT. 0 .OR. IMETH .GT. 2) THEN IERR=3 RETURN ENDIF C C Initialize the parameters KADAPT and OMEGA depending on the method C IF(IMETH .EQ. 0) THEN OMEGA=1.0D0 ELSE IF(IMETH .EQ. 1) THEN KADAPT=KMAX ELSE IF(IMETH .EQ. 2) THEN KADAPT=KMAX OMEGA=1.0D0 ENDIF C C Compute the machine constant and initialize the relative error bound C FMACHP=1.0D0 10 FMACHP=0.5D0*FMACHP IF(MACHPD(1.0D0+FMACHP) .EQ. 1) GOTO 10 RELERR=FMACHP*8.0D0 C C Initialize C Q=1.0D0 ITNUMB=0 C C Check whether A is singular; if so, set IERR = 4. C DO 20 I=1,N SUM=DABS(A(I,1)) DO 30 K=2,N SUM=SUM+DABS(A(I,K)) 30 CONTINUE IF(SUM .EQ. 0.0D0) THEN IERR=4 RETURN ELSE IF(DABS(A(I,I))/SUM .LT. RELERR) THEN IERR=4 RETURN ENDIF 20 CONTINUE C C Redefine the entries in A and B: A(I,J) := A(I,J)/A(I,I) C and B(I) := B(I)/A(I,I) . C DO 40 I=1,N HELP=1.0D0/A(I,I) DO 50 J=1,N A(I,J)=A(I,J)*HELP 50 CONTINUE B(I)=B(I)*HELP 40 CONTINUE C C Check for convergence C IF(ISWITC .NE. 0) THEN CALL CONV(ISWITC,A,N,IA,IERR) IF(IERR .EQ. 0) RETURN ENDIF C C The vector RES serves as auxiliary storage for the previous solution C vektor. Initially RES contains the staring vector. C DO 60 I=1,N RES(I)=X(I) 60 CONTINUE C C One iteration with the Gauá-Seidel method gives the first iterate X C CALL GAUSEI(A,N,IA,B,OMEGA,X) C C Up the iteration counter C ITNUMB=ITNUMB+1 C C Compute the difference of the last two iterates C DO 70 I=1,N WORK(I,1)=X(I)-RES(I) 70 CONTINUE C C Iteration loop for the chosen method C DO 80 K=1,KMAX-1 C C Check break-off criterion C CALL MNORM(WORK(1,1),N,DIFFN) CALL MNORM(X,N,XN) IF(DIFFN .LE. EPS*XN) THEN IERR=1 ITNUMB=K CALL RESID(A,N,IA,B,X,RES) RETURN ENDIF IF(K .EQ. KMAX-1) THEN ITNUMB=KMAX IERR=2 CALL RESID(A,N,IA,B,X,RES) RETURN ENDIF C C RES contains the previous iterate C DO 90 I=1,N RES(I)=X(I) 90 CONTINUE C C One iteration step using Gauá-Seidel for a fixed OMEGA C CALL GAUSEI(A,N,IA,B,OMEGA,X) C C Compute the difference of the last two iterates C DO 100 I=1,N WORK(I,2)=X(I)-RES(I) 100 CONTINUE C C If the number of performed iterations K is divisible by KADAPT, C then we compute Q in order to adjust the relaxation parameter; C Q is an estimate of the spectral radius of the iteration matrix. C IF(MOD(K,KADAPT) .EQ. 0) THEN DO 110 I=1,N IF(DABS(WORK(I,1)) .LT. FMACHP) THEN WORK(I,3)=1.0D0 ELSE WORK(I,3)=WORK(I,2)/WORK(I,1) ENDIF 110 CONTINUE CALL MNORM(WORK(1,3),N,Q) C C If Q > 1, then the iteration counter is upped by one and C the next Gauá-Seidel step is executed; otherwise a new C relaxation parameter is calculated. C IF(Q .LE. 1.0D0) THEN Q=MAX(Q,OMEGA-1.0D0) OMEGA=2.0D0/(1.0D0+DSQRT(1.0D0-((Q+OMEGA-1.0D0) * /OMEGA)**2/Q)) ENDIF ENDIF C C The difference vector of the last two iterations is replaced C by the one of the previous two iterations for the approximate solution C DO 120 I=1,N WORK(I,1)=WORK(I,2) 120 CONTINUE 80 CONTINUE END C C C[BA*) C[LE*) SUBROUTINE GAUSEI(A,N,IA,B,OMEGA,X) C[IX{GAUSEI}*) C C***************************************************************** C * C This subroutine performs one iteration with the Gauá-Seidel * C method for a given relaxation parameter. * C[BE*) C * C * C INPUT PARAMETERS: * C ================= * C A : 2-dim. array A(1:IA, 1:N), that contains the * C modified system matrix A : A(I,J)=A(I,J)/A(I,I) for * C I,J=1, ..., N * C N : order of the system * C IA : leading dimension of A, as specified in the calling * C program * C B : N-vector B(1:N) with the modified right hand side: * C B(I)=B(I)/A(I,I); I=1, ..., N * C OMEGA : relaxation parameter * C X : N-vector X(1:N) containing the starting vector for * C the iteration * C * C * C OUTPUT PARAMETERS: * C ================== * C X : N-vector X(1:N) containing the next iteration vector * C * C----------------------------------------------------------------* C * C Required subroutines: none * C * C***************************************************************** C * C Author : Gisela Engeln-Mllges * C Date : 06.09.1992 * C Source : FORTRAN 77 * C * C[BA*) C***************************************************************** C[BE*) C DOUBLE PRECISION A(1:IA,1:N),B(1:N),X(1:N),OMEGA,S C DO 10 I=1,N S=B(I) DO 20 J=1,N S=S-A(I,J)*X(J) 20 CONTINUE X(I)=X(I)+OMEGA*S 10 CONTINUE RETURN END C C C[BA*) C[LE*) SUBROUTINE MNORM(X,N,XNORM) C[IX{MNORM}*) C C***************************************************************** C * C This subroutine calculates the maximum norm XNORM of an * C N-vector X. * C * C----------------------------------------------------------------* C[BE*) C * C Required subroutines: none * C * C***************************************************************** C * C Author : Gisela Engeln-Mllges * C Date : 06.09.1992 * C Source : FORTRAN 77 * C * C***************************************************************** C DOUBLE PRECISION X(1:N),XNORM C XNORM=DABS(X(1)) DO 10 I=2,N XNORM=DMAX1(XNORM,DABS(X(I))) 10 CONTINUE RETURN END C C C[BA*) C[LE*) SUBROUTINE CONV(ISWITC,A,N,IA,IERR) C[IX{CONV}*) C C***************************************************************** C * C This subroutine helps check convergence. * C[BE*) C * C * C INPUT PARAMETERS: * C ================= * C ISWITC : Parameter that determines the convergence criterion * C to be checked: * C = 0, none * C = 1, row sum criterion * C = 2, column sum criterion * C = 3, criterion of Schmidt and v. Mises * C A : 2-dim. array A(1:IA, 1:N), containing the matrix for * C which we want to check convergence of the iterates * C from the various SOR algorithms * C N : order of the matrix A * C IA : leading dimension of A, as prescribed in the calling * C program * C * C * C OUTPUT PARAMETERS: * C ================== * C IERR : error parameter: * C = 0, the desired convergence criterion has not been * C met * C = 1, the desired criterion is satified * C * C----------------------------------------------------------------* C * C Required subroutines: none * C * C***************************************************************** C * C Author : Gisela Engeln-Mllges * C Date : 06.09.1992 * C Source : FORTRAN 77 * C * C[BA*) C***************************************************************** C[BE*) C DOUBLE PRECISION A(1:IA,1:N),SUM C C Row sum criterion C IF(ISWITC .EQ. 1) THEN DO 10 I=1,N SUM=-1.0D0 DO 20 J=1,N SUM=SUM+DABS(A(I,J)) 20 CONTINUE IF(SUM .LT. 1.0D0) THEN IERR=1 ELSE IERR=0 RETURN ENDIF 10 CONTINUE C C Column sum criterion C ELSE IF(ISWITC .EQ. 2) THEN DO 30 J=1,N SUM=-1.0D0 DO 40 I=1,N SUM=SUM+DABS(A(I,J)) 40 CONTINUE IF(SUM .LT. 1.0D0) THEN IERR=1 ELSE IERR=0 RETURN ENDIF 30 CONTINUE C C Criterion of Schmidt and v. Mises C ELSE IF(ISWITC .EQ. 3) THEN SUM=-N DO 50 I=1,N DO 60 J=1,N SUM=SUM+A(I,J)*A(I,J) 60 CONTINUE 50 CONTINUE SUM=DSQRT(SUM) IF(SUM .LT. 1.0D0) THEN IERR=1 ELSE IERR=0 RETURN ENDIF ENDIF END C C C[BA*) C[LE*) SUBROUTINE RESID(A,N,IA,B,X,RES) C[IX{RESID}*) C C***************************************************************** C * C This subroutine computes the residuum RES = B - AX, where * C both A and B are given in modified form. * C * C----------------------------------------------------------------* C[BE*) C * C Required subroutines: none * C * C***************************************************************** C * C Author : Gisela Engeln-Mllges * C Date : 09.06.1992 * C Source : FORTRAN 77 * C * C***************************************************************** C DOUBLE PRECISION A(1:IA,1:N), B(1:N), X(1:N), RES(1:N),DSUM C DO 10 I=1,N DSUM=B(I) DO 20 J=1,N DSUM=DSUM-A(I,J)*X(J) 20 CONTINUE RES(I)=DSUM 10 CONTINUE RETURN END c C[KA{F 0}{Auxiliary Library}{Auxiliary Library}*) INTEGER FUNCTION MACHPD(X) C[IX{MACHPD}*) DOUBLE PRECISION X MACHPD=0 IF (1.0D0 .LT. X) MACHPD=1 RETURN END