1 : C CG2 SOURCE GOUNAND 11/10/07 21:15:10 7149
2 : SUBROUTINE CG2(KMORS,KISA,B,X,
3 : $ KPREC,INVDIA,ILUM,ILUI,
4 : $ LRES,LNMV,
5 : $ R,P,Q,Z,
6 : $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
7 : IMPLICIT INTEGER(I-N)
8 : IMPLICIT REAL*8 (A-H,O-Z)
9 : C***********************************************************************
10 : C NOM : CG2
11 : C DESCRIPTION :
12 : C Résolution d'un système linéaire Ax=b
13 : C par la méthode du gradient conjugué préconditionné.
14 : C
15 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
16 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
17 : C mél : gounand@semt2.smts.cea.fr
18 : C REFERENCE (bibtex-like) :
19 : C @BOOK{templates,
20 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
21 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
22 : C H. Van der Vorst},
23 : C TITLE={Templates for the Solution of Linear Systems :
24 : C Building Blocks for Iterative Methods},
25 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
26 : C -> URL : http://www.netlib.org/templates/Templates.html
27 : C***********************************************************************
28 : C APPELES (BLAS 1) : GCOPY, GAXPY (subroutines)
29 : C GNRM2, GDOT (functions)
30 : C APPELES (Calcul) : GMOMV : produit Matrice Morse - Vecteur
31 : C PSOLVE : Preconditionner solve
32 : C***********************************************************************
33 : C ENTREES :
34 : C ENTREES/SORTIES :
35 : C SORTIES :
36 : C CODE RETOUR (IRET) : 0 si ok
37 : C <0 si problème (non utilisé !)
38 : C >0 si l'algorithme n'a pas convergé
39 : C***********************************************************************
40 : C VERSION : v1, 09/02/06
41 : C HISTORIQUE : v1, 09/02/06, création
42 : C HISTORIQUE :
43 : C HISTORIQUE :
44 : C***********************************************************************
45 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
46 : C en cas de modification de ce sous-programme afin de faciliter
47 : C la maintenance !
48 : C***********************************************************************
49 : *
50 : * .. Includes .. a commenter si mise au point ok
51 : -INC CCOPTIO
52 : -INC CCREEL
53 : -INC SMLREEL
54 : POINTEUR LRES.MLREEL
55 : -INC SMLENTI
56 : POINTEUR LNMV.MLENTI
57 : -INC SMMATRIK
58 : POINTEUR KMORS.PMORS
59 : POINTEUR KISA.IZA
60 : POINTEUR INVDIA.IZA
61 : POINTEUR ILUM.PMORS
62 : POINTEUR ILUI.IZA
63 : POINTEUR B.IZA
64 : POINTEUR X.IZA
65 : * .. Work Vectors for CG
66 : POINTEUR R.IZA,P.IZA,Q.IZA,Z.IZA
67 : * .. Scalar Arguments ..
68 : INTEGER ITER, KPREC, IMPR, IRET
69 : REAL*8 RESID
70 : *
71 : *
72 : * .. Variables locales
73 : * .. Parameters
74 : REAL*8 ZERO , ONE
75 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
76 : * ..
77 : INTEGER MAXIT
78 : REAL*8 ALPHA,BETA,RHO,RHO1,TOL,GNRM2,OMEGA
79 : REAL*8 GDOT,BNRM2
80 : * ..
81 : * .. External subroutines and functions..
82 : EXTERNAL GAXPY,GCOPY,GDOT,GNRM2
83 : * ..
84 : * .. Executable Statements ..
85 : *
86 :
87 : IRET = 0
88 : INMV = 0
89 : IWARN = 0
90 : IGRAN = 0
91 : MGRAN = 5
92 : MAXIT = ITER
93 : TOL = RESID
94 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cg2.eso'
95 : *
96 : * Set parameter tolerances.
97 : *
98 : RHOTOL = BRTOL
99 : OMETOL = BRTOL
100 : * RRMAX = 10.D0
101 : RRMAX = 1.D0/(XZPREC**0.25D0)
102 : *
103 : * Calcul du résidu initial
104 : *
105 : C r(0)=b
105 : CALL GCOPY(B,R)
106 :
107 : C r(0)=b-Ax
107 : CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
108 : INMV=INMV+1
109 : *
110 : ITER = 0
111 : BNRM2 = GNRM2(B)
112 : IF(IMPR.GT.5) THEN
113 : WRITE(IOIMP,*) '||B||=',BNRM2
114 : ENDIF
115 : RNRM2 = GNRM2(R)
116 : IF (ICALRS.EQ.1) BNRM2=RNRM2
117 : IF (BNRM2.LT.XPETIT) BNRM2 = ONE
118 : RESID = RNRM2 / BNRM2
119 : IF(IMPR.GT.5) THEN
120 : WRITE(IOIMP,*) 'Initial residual=',RESID
121 : ENDIF
122 : RESMAX = RESID * RRMAX
123 : LRES.PROG(ITER+1)=RESID
124 : LNMV.LECT(ITER+1)=INMV
125 : *
126 : IF (RESID.LE.TOL) THEN
127 : IRET = 0
128 : GOTO 30
129 : ENDIF
130 :
131 : *
132 : 10 CONTINUE
133 : *
134 : * Perform (Preconditioned) Conjugate Gradient Iteration.
135 : *
136 : ITER = ITER +1
137 : IF(IMPR.GT.7) THEN
138 : WRITE(IOIMP,*) 'ITER=',ITER
139 : ENDIF
140 : *
141 : * Preconditioner Solve
142 : C Mz(i-1)=r(i-1)
142 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
143 : $ ,Z,R)
144 : C rho(i-1)=r(i-1)'z(i-1)
145 : RHO = GDOT(R,Z,IDDOT)
146 : * WRITE(IOIMP,*) ' RHO=',RHO
147 : IF (ABS(RHO).LT.RHOTOL) THEN
148 : RHO=SIGN(XZPREC,RHO)
149 : IWARN=IWARN+1
150 : ENDIF
151 : *
152 : * Compute direction vector P
153 : *
154 : IF (ITER.GT.1) THEN
155 : C beta(i-1)=rho(i-1)/rho(i-2)
156 : BETA = RHO / RHO1
157 : C p(i)=z(i-1)+beta(i-1)p(i-1)
157 : CALL GAXPY(BETA,P,Z)
157 : CALL GCOPY(Z,P)
158 : ELSE
159 : C p(1)=z(0)
159 : CALL GCOPY(Z,P)
160 : ENDIF
161 : *
162 : * Compute scalar ALPHA (save A*p to q)
163 : *
164 : C q(i)=Ap(i)
164 : CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,P,ZERO,Q)
165 : INMV=INMV+1
166 : C alpha(i)=rho(i-1)/(p(i)'q(i))
167 : OMEGA = GDOT(P,Q,IDDOT)
168 : * WRITE(IOIMP,*) ' OMEGA=',OMEGA
169 : IF (ABS(OMEGA).LT.OMETOL) THEN
170 : OMEGA=SIGN(XZPREC,OMEGA)
171 : IWARN=IWARN+1
172 : ENDIF
173 : ALPHA = RHO / OMEGA
174 : *
175 : * Compute current solution vector x
176 : *
177 : C x(i)=x(i-1)+alpha(i)p(i)
177 : CALL GAXPY(ALPHA,P,X)
178 : *
179 : * Compute residual vector R, find norm,
180 : * then check for tolerance
181 : *
182 : C r(i)=r(i-1)-alpha(i)r(i)
182 : CALL GAXPY(-ALPHA,Q,R)
183 : *
184 : RESID = GNRM2(R) / BNRM2
185 : IF(IMPR.GT.7) THEN
186 : WRITE(IOIMP,*) 'Résidu=',RESID
187 : ENDIF
188 : LRES.PROG(ITER+1)=RESID
189 : LNMV.LECT(ITER+1)=INMV
190 : IF (RESID.GT.RESMAX) THEN
191 : IGRAN=IGRAN+1
192 : ELSE
193 : IGRAN=0
194 : ENDIF
195 : *
196 : * Iteration successful
197 : *
198 : IF (RESID.LE.TOL) THEN
199 : IRET = 0
200 : *
201 : * Iteration fails
202 : *
203 : * ELSE IF (ITER.EQ.MAXIT.OR.IGRAN.GT.MGRAN) THEN
204 : ELSE IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
205 : IRET = 1
206 : *
207 : * Do some more
208 : *
209 : ELSE
210 : RHO1 = RHO
211 : GOTO 10
212 : ENDIF
213 : *
214 : * Normal termination
215 : *
216 : 30 CONTINUE
217 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
218 : WRITE(IOIMP,*) 'CG nearly failed ',IWARN
219 : $ ,' times.'
220 : ENDIF
221 : RETURN
222 : *
223 : * End of CG2
224 : *
225 : END
226 :
227 :
228 :
229 :
230 :
231 :
232 :
233 :
234 :
235 :
© Cast3M 2003 - All rights reserved.
Disclaimer