C GMRD SOURCE CHAT 05/01/13 00:18:38 5004 SUBROUTINE GMRD(N,NNZ,ROWPTR,COLIND,VAL, $ B,X, $ MINV, $ LRES, $ R,W,AV, $ Y,GCOS,GSIN,S, $ V,H, $ ITER,RESTRT,RESID,IMPR,IRET) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C NOM : GMRD C DESCRIPTION : C Résolution d'un système linéaire Ax=b C par une méthode GMRES(m) (restarted Generalized Minimal C Residual) préconditionné par la diagonale de A (Jacobi). C C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S) C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF) C mél : gounand@semt2.smts.cea.fr C REFERENCE (bibtex-like) : C @BOOK{templates, C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato, C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine, C H. Van der Vorst}, C TITLE={Templates for the Solution of Linear Systems : C Building Blocks for Iterative Methods}, C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} } C -> URL : http://www.netlib.org/templates/Templates.html C*********************************************************************** C APPELES (BLAS 1) : DCOPY, DAXPY (subroutines) C DNRM2, DDOT, DSCAL (functions) C APPELES (BLAS 2) : DTRSV (subroutines) C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur C*********************************************************************** C ENTREES : N, NNZ, C ROWPTR, COLIND, VAL, C B, MINV, RESTRT, IMPR C ENTREES/SORTIES : X, C R,W,AV, C Y,GCOS,GSIN,S, C V,H, C ITER, RESID C SORTIES : LRES, IRET C CODE RETOUR (IRET) : 0 si ok C <0 si problème (cf. BRTOL) C >0 si l'algorithme n'a pas convergé C N : nombre de degrés de liberté C NNZ : nombre de valeurs non nulles de la matrice Morse C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne C et valeurs de la matrice Morse C B : vecteur second membre C MINV : vecteur contenant les inverses des termes C diagonaux de A C RESTRT : INTEGER C paramètre de redémarrage (restart) (i.e. m) pour C la méthode GMRES(m). (cf. prgmr.eso) C IMPR : niveau d'impression (0..4) C X : vecteur des inconnues. C En entrée, contient l'estimation x(0). C En sortie, la solution trouvée (si la méthode C a convergé). C R,W,AV : vecteurs de travail initialisés et détruits C Y,GCOS,GSIN,S dans prgmr. C V,H : tableaux de travail... C ITER : type INTEGER. C En entrée, contient le nombre maximum C d'itérations à effectuer. C En sortie, contient le nombre d'itérations C réellement effectuées. C RESID : type REAL*8. C En entrée, la valeur maximale du résidu à C convergence de l'algorithme mesurée comme suit : C norme[L2](b - A*x) / norme[L2]( b ) C En sortie, la valeur effective de cette mesure. C LRES : pointeur sur segment MLREEL (include SMLREEL) C contient l'historique de la valeur de RESID en C fonction du nombre d'itérations effectuées. C C*********************************************************************** C VERSION : v1, 16/06/98, version initiale C HISTORIQUE : v1, 16/06/98, création C HISTORIQUE : C HISTORIQUE : C*********************************************************************** C Prière de PRENDRE LE TEMPS de compléter les commentaires C en cas de modification de ce sous-programme afin de faciliter C la maintenance ! C*********************************************************************** * * .. Includes .. a commenter si mise au point ok -INC CCOPTIO INTEGER IMPR * .. Scalar Arguments .. INTEGER N,NNZ INTEGER ITER,RESTRT,IRET REAL*8 RESID * .. * .. Array Arguments .. * .. Matrices stockées en Morse INTEGER ROWPTR(N+1) INTEGER COLIND(NNZ) REAL*8 VAL(NNZ) * .. Vecteurs REAL*8 B(N) REAL*8 X(N) REAL*8 MINV(N) REAL*8 LRES(ITER) * .. Tableaux de travail REAL*8 R(N),W(N),AV(N) REAL*8 Y(RESTRT),GCOS(RESTRT),GSIN(RESTRT) REAL*8 S(RESTRT+1) REAL*8 V(N,RESTRT+1) REAL*8 H(RESTRT+1,RESTRT) * * * .. Variables locales * .. Parameters REAL*8 ZERO,ONE PARAMETER (ZERO=0.0D0, ONE=1.0D0) * .. INTEGER I,MAXIT,J,K,IN REAL*8 BNRM2,TOL,DDOT,DNRM2 REAL*8 H1,H2,RAY * .. * .. External subroutines and functions.. EXTERNAL DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DTRSV EXTERNAL DMOMV * .. * .. Executable Statements .. * * IRET = 0 MAXIT = ITER TOL = RESID IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans gmrd.eso' * * Set initial residual. * C r=b CALL DCOPY(N,B,1,R,1) IF (DNRM2(N,X,1).NE.ZERO) THEN C r=b-Ax CALL DMOMV ('No transpose',N,-ONE, $ NNZ,ROWPTR,COLIND,VAL, $ X,ONE,R) ENDIF * ITER = 0 BNRM2 = DNRM2(N,B,1) IF(IMPR.GT.6) THEN WRITE(IOIMP,*) '||B||=',BNRM2 ENDIF IF (BNRM2.EQ.ZERO) BNRM2 = ONE RESID = DNRM2(N,R,1) / BNRM2 * RESID = DNRM2(N,R,1) IF(IMPR.GT.6) THEN WRITE(IOIMP,*) 'Résidu initial=',RESID ENDIF IF (RESID.LE.TOL) GOTO 30 * 10 CONTINUE * * Construct the first column of V, and initialize S to the * elementary vector E1 scaled by RNORM. * C Preconditionner solve : Mv(1)=r C CALL DCOPY(N,R,1,V(1,1),1) DO 100 IN=1,N V(IN,1)=R(IN)*MINV(IN) 100 CONTINUE C s=||v(1)|| e(1) C v(1)=v(1) / ||v(1)|| S(1) = DNRM2(N,V(1,1),1) CALL DSCAL(N,ONE/S(1),V(1,1),1) * DO 50 I = 1, RESTRT ITER = ITER + 1 C av=A v(i) CALL DMOMV ('No transpose',N,ONE, $ NNZ,ROWPTR,COLIND,VAL, $ V(1,I),ZERO,AV) C Preconditionner solve : Mw=av C CALL DCOPY(N,AV,1,W,1) DO 200 IN=1,N W(IN)=AV(IN)*MINV(IN) 200 CONTINUE * * Construct I-th column of H so that it is orthonormal to * the previous I-1 columns. * C CALL ORTHOH(I,N,RESTRT,WORK2(1,I+H-1),WORK(1,V),WORK(1,W)) DO 60 K = 1, I H(K,I) = DDOT(N,W,1,V(1,K),1) CALL DAXPY(N,-H(K,I),V(1,K),1,W,1) 60 CONTINUE H(I+1,I)=DNRM2(N,W,1) CALL DCOPY(N,W,1,V(1,I+1),1) CALL DSCAL(N,ONE/H(I+1,I),V(1,I+1),1) * C IF (I.GT.1) * * Apply Givens rotations to the I-th column of H. This * effectively reduces the Hessenberg matrix to upper * triangular form during the RESTRT iterations. * C $ CALL APPGIV(I,RESTRT,WORK2(1,I+H-1),WORK2(1,GIV)) DO 70 K = 1, I-1 H1= (GCOS(K)*H(K,I)) + (GSIN(K)*H(K+1,I)) H2=-(GSIN(K)*H(K,I)) + (GCOS(K)*H(K+1,I)) H(K,I) =H1 H(K+1,I)=H2 70 CONTINUE * * Approximate residual norm. Check tolerance. If okay, compute * final approximation vector X and quit. * RAY=SQRT((H(I,I)*H(I,I)) + (H(I+1,I)*H(I+1,I))) GCOS(I) = H(I,I) / RAY GSIN(I) = H(I+1,I) / RAY H(I,I)=RAY H(I+1,I)=ZERO S(I+1)=-GSIN(I)*S(I) S(I) = GCOS(I)*S(I) RESID=ABS(S(I+1)) / BNRM2 LRES(ITER)=RESID IF (RESID.LE.TOL.OR.ITER.EQ.MAXIT) THEN * Solve the system Hy=S * DO 501 J=I,1,-1 * Y(J)=S(J) / H(J,J) * DO 502 K=J-1,1,-1 * S(K)=S(K)-H(K,J)*Y(J) * 502 CONTINUE * 501 CONTINUE CALL DCOPY(I,S,1,Y,1) CALL DTRSV('UPPER','NOTRANS','NONUNIT',I,H,RESTRT+1,Y,1) * Updating solution C CALL UPDATE(I,N,RESTRT,X,WORK2(1,H),WORK(1,Y), C $ WORK(1,S),WORK(1,V)) DO 503 J=I,1,-1 CALL DAXPY(N,Y(J),V(1,J),1,X,1) 503 CONTINUE IF (RESID.LE.TOL) THEN GOTO 30 ENDIF IF (ITER.EQ.MAXIT) THEN GOTO 40 ENDIF ENDIF 50 CONTINUE * * Compute current solution vector X. * * Solve the system Hy=S * DO 301 J=RESTRT,1,-1 * Y(J)=S(J) / H(J,J) * DO 302 K=J-1,1,-1 * S(K)=S(K)-H(K,J)*Y(J) * 302 CONTINUE * 301 CONTINUE CALL DCOPY(RESTRT,S,1,Y,1) CALL DTRSV('UPPER','NOTRANS','NONUNIT', $ RESTRT,H,RESTRT+1,Y,1) * Updating solution C CALL UPDATE(RESTRT,N,RESTRT,X,WORK2(1,H),WORK(1,Y), C $ WORK(1,S),WORK(1,V)) DO 303 J=RESTRT,1,-1 CALL DAXPY(N,Y(J),V(1,J),1,X,1) 303 CONTINUE * * Compute residual vector R, find norm, * then check for tolerance. * CALL DCOPY(N,B,1,R,1) CALL DMOMV ('No transpose',N,-ONE, $ NNZ,ROWPTR,COLIND,VAL, $ X,ONE,R) S(RESTRT+1) = DNRM2(N,R,1) RESID = S(RESTRT+1) / BNRM2 LRES(ITER)=RESID IF (RESID.LE.TOL) GOTO 30 * * Iteration fails. * IF (ITER.EQ.MAXIT) THEN GOTO 40 ENDIF * * Do some more * GOTO 10 * * Iteration successful * 30 CONTINUE IRET = 0 RETURN * * Iteration fails. * 40 CONTINUE IRET = 1 RETURN * * End of GMRD. * END