C CGDILU SOURCE CHAT 05/01/12 21:53:22 5004 SUBROUTINE CGDILU(N,NNZ,ROWPTR,COLIND,VAL, $ B,X, $ DINV, $ LRES, $ R,P,Q,Z, $ ITER,RESID,IMPR,IRET) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C NOM : CGDILU C DESCRIPTION : C Résolution d'un système linéaire Ax=b C par la méthode du gradient conjugué C préconditionné par la factorisation D-ILU de A. 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 (functions) C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur C PSDILU : Preconditionner solve C (Montée-Descente adaptée à D-ILU) C*********************************************************************** C ENTREES : N, NNZ, C ROWPTR, COLIND, VAL, C B, DINV, IMPR C ENTREES/SORTIES : X, R, P, Q, Z, ITER, RESID C SORTIES : LRES, IRET C CODE RETOUR (IRET) : 0 si ok C <0 si problème (non utilisé !) 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 DINV : vecteur contenant les inverses des pivots C de la factorisation D-ILU de A. 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,P,Q,Z: vecteurs de travail initialisés et détruits C dans prcg (on peut en utiliser moins ?). 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, 02/04/98, version initiale C HISTORIQUE : v1, 02/04/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,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 DINV(N) REAL*8 LRES(ITER) REAL*8 R(N),P(N),Q(N),Z(N) * * * .. Variables locales * .. Parameters REAL*8 ZERO , ONE PARAMETER (ZERO=0.0D0, ONE=1.0D0) * .. INTEGER MAXIT REAL*8 ALPHA,BETA,RHO,RHO1,TOL,BNRM2 REAL*8 DDOT,DNRM2 * .. * .. External subroutines and functions.. EXTERNAL DAXPY,DCOPY,DDOT,DNRM2 * EXTERNAL STOPTEST2 * .. * .. Executable Statements .. * IRET = 0 MAXIT = ITER TOL = RESID IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgdilu.eso' * * Calcul du résidu initial * C r(0)=b CALL DCOPY(N,B,1,R,1) IF(IMPR.GT.6) THEN WRITE(IOIMP,*) 'r(0)=b' ENDIF IF (DNRM2(N,X,1).NE.ZERO) THEN C r(0)=b-Ax CALL DMOMV ('No transpose',N,-ONE, $ NNZ,ROWPTR,COLIND,VAL, $ X,ONE,R) IF(IMPR.GT.6) THEN WRITE(IOIMP,*) 'r(0)=b-Ax' ENDIF * Je ne comprends pas cette ligne de code qui * vient de CGREVCOM.f * IF (DNRM2(N,R,1).LT.TOL) GOTO 30 ENDIF * ITER = 0 * On ne peut pas normer par B à cause des 1.E30 BNRM2 = DNRM2(N,B,1) IF(IMPR.GT.5) 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.5) THEN WRITE(IOIMP,*) 'Résidu initial=',RESID ENDIF IF (RESID.LE.TOL) THEN IRET = 0 GOTO 30 ENDIF * 10 CONTINUE * * Perform (Preconditioned) Conjugate Gradient Iteration. * ITER = ITER +1 IF(IMPR.GT.6) THEN WRITE(IOIMP,*) 'ITER=',ITER ENDIF * * Preconditioner Solve (D-ILU) * C Mz(i-1)=r(i-1) CALL PSDILU(N,NNZ,ROWPTR,COLIND,VAL, $ DINV, $ Z,R) * * C rho(i-1)=r(i-1)'z(i-1) RHO = DDOT(N,R,1,Z,1) * * Compute direction vector P * IF (ITER.GT.1) THEN C beta(i-1)=rho(i-1)/rho(i-2) BETA = RHO / RHO1 C p(i)=z(i-1)+beta(i-1)p(i-1) CALL DAXPY(N,BETA,P,1,Z,1) CALL DCOPY(N,Z,1,P,1) ELSE C p(1)=z(0) CALL DCOPY(N,Z,1,P,1) ENDIF * * Compute scalar ALPHA (save A*p to q) * C q(i)=Ap(i) CALL DMOMV('No transpose',N,ONE, $ NNZ,ROWPTR,COLIND,VAL, $ P,ZERO,Q) C alpha(i)=rho(i-1)/(p(i)'q(i)) ALPHA = RHO / DDOT(N,P,1,Q,1) * * Compute current solution vector x * C x(i)=x(i-1)+alpha(i)p(i) CALL DAXPY(N,ALPHA,P,1,X,1) IF(IMPR.GT.7) THEN WRITE(IOIMP,*) 'X=',(X(III),III=1,N) ENDIF * * Compute residual vector R, find norm, * then check for tolerance * C r(i)=r(i-1)-alpha(i)r(i) CALL DAXPY(N,-ALPHA,Q,1,R,1) * RESID = DNRM2(N,R,1) / BNRM2 * RESID = DNRM2(N,R,1) IF(IMPR.GT.6) THEN WRITE(IOIMP,*) 'Résidu=',RESID ENDIF LRES(ITER)=RESID * * Iteration successful * IF (RESID.LE.TOL) THEN IRET = 0 * * Iteration fails * ELSE IF (ITER.EQ.MAXIT) THEN IRET = 1 * * Do some more * ELSE RHO1 = RHO GOTO 10 ENDIF * * Normal termination * 30 CONTINUE RETURN * * End of CGDILU. * END