1 : C CGSN SOURCE GOUNAND 11/10/07 21:15:11 7149
2 : SUBROUTINE CGSN(KMORS,KISA,B,X,
3 : $ KPREC,INVDIA,ILUM,ILUI,
4 : $ LRES,LNMV,
5 : $ R,RTLD,P,PHAT,
6 : $ Q,QHAT,UHAT,
7 : $ XP,BP,
8 : $ ITER,INMV,BRTOL,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
9 : IMPLICIT REAL*8 (A-H,O-Z)
10 : IMPLICIT INTEGER(I-N)
11 : C***********************************************************************
12 : C NOM : CGSN
13 : C DESCRIPTION :
14 : C Résolution d'un système linéaire Ax=b
15 : C par une méthode de gradient conjugué au carré (CGS) avec la
16 : C stratégie de stabilisation de Neumaier préconditionnée ou non.
17 : C
18 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
19 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
20 : C mél : gounand@semt2.smts.cea.fr
21 : C REFERENCE (bibtex-like) :
22 : C @BOOK{templates,
23 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
24 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
25 : C H. Van der Vorst},
26 : C TITLE={Templates for the Solution of Linear Systems :
27 : C Building Blocks for Iterative Methods},
28 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
29 : C -> URL : http://www.netlib.org/templates/Templates.html
30 : C***********************************************************************
31 : C APPELES (BLAS 1) : GCOPY, GAXPY (subroutines)
32 : C GNRM2, GDOT (functions)
33 : C APPELES (Calcul) : GMOMV : produit Matrice Morse - Vecteur
34 : C PSOLVE : Preconditionner solve
35 : C***********************************************************************
36 : C ENTREES :
37 : C ENTREES/SORTIES :
38 : C SORTIES :
39 : C CODE RETOUR (IRET) : 0 si ok
40 : C <0 si problème (non utilisé !)
41 : C >0 si l'algorithme n'a pas convergé
42 : C***********************************************************************
43 : C VERSION : v1, 10/02/06
44 : C HISTORIQUE : v1, 10/02/06, création
45 : C HISTORIQUE :
46 : C HISTORIQUE :
47 : C***********************************************************************
48 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
49 : C en cas de modification de ce sous-programme afin de faciliter
50 : C la maintenance !
51 : C***********************************************************************
52 : *
53 : * .. Includes .. a commenter si mise au point ok
54 : -INC CCREEL
55 : -INC CCOPTIO
56 : -INC SMLREEL
57 : POINTEUR LRES.MLREEL
58 : -INC SMLENTI
59 : POINTEUR LNMV.MLENTI
60 : -INC SMMATRIK
61 : POINTEUR KMORS.PMORS
62 : POINTEUR KISA.IZA
63 : POINTEUR INVDIA.IZA
64 : POINTEUR ILUM.PMORS
65 : POINTEUR ILUI.IZA
66 : POINTEUR B.IZA
67 : POINTEUR X.IZA
68 : * .. Work Vectors for CGSN
69 : POINTEUR R.IZA,RTLD.IZA,P.IZA,PHAT.IZA
70 : POINTEUR Q.IZA,QHAT.IZA,UHAT.IZA
71 : POINTEUR U.IZA,VHAT.IZA
72 : * Note QHAT and U, UHAT,VHAT share the same memory space
73 : POINTEUR XP.IZA,BP.IZA
74 : *
75 : INTEGER IMPR
76 : * .. Scalar Arguments ..
77 : INTEGER N,NNZ
78 : INTEGER ITER,IRET
79 : REAL*8 BRTOL,RESID
80 : * ..
81 : *
82 : *
83 : * .. Variables locales
84 : * .. Parameters
85 : REAL*8 ZERO,ONE
86 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
87 : * ..
88 : INTEGER MAXIT,III
89 : REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
90 : REAL*8 RHOTOL
91 : REAL*8 GDOT,GNRM2,RNRM2
92 : REAL*8 MU,MUP
93 : * ..
94 : * .. External subroutines and functions..
95 : EXTERNAL GAXPY,GCOPY,GDOT,GNRM2,GSCAL
96 : * ..
97 : * .. Executable Statements ..
98 : *
99 : *
100 : IRET = 0
101 : INMV = 0
102 : IWARN = 0
103 : MAXIT = ITER
104 : TOL = RESID
105 : *
106 : * Same memory space
107 : U=QHAT
108 : VHAT=UHAT
109 : *
110 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgsn.eso'
111 : *
112 : * Set breakdown parameter tolerances.
113 : *
114 : RHOTOL = BRTOL
115 : TAUTOL = BRTOL
116 : *
117 : * Set initial residual.
118 : *
119 : C r(0)=b
119 : CALL GCOPY(B,R)
120 : C r(0)=b-Ax
121 : * IF (GNRM2(X).NE.ZERO) THEN
122 : C r(0)=b-Ax
122 : * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
123 : * ENDIF
123 : CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
124 : INMV=INMV+1
125 : *
126 : ITER = 0
127 : BNRM2 = GNRM2(B)
128 : IF(IMPR.GT.6) THEN
129 : WRITE(IOIMP,*) '||B||=',BNRM2
130 : ENDIF
131 : RNRM2 = GNRM2(R)
132 : IF (ICALRS.EQ.1) BNRM2=RNRM2
133 : *
134 : IF (BNRM2.LT.XPETIT) BNRM2 = ONE
135 : RESID = RNRM2 / BNRM2
136 : IF(IMPR.GT.6) THEN
137 : WRITE(IOIMP,*) 'Résidu initial=',RESID
138 : ENDIF
139 : LRES.PROG(ITER+1)=RESID
140 : LNMV.LECT(ITER+1)=INMV
141 : IF (RESID.LE.TOL) GOTO 30
142 : C r[tilde]=r(0)
142 : CALL GCOPY(R,RTLD)
143 : C Neumaier's Strategy
144 : C x=x(0) x'=0
145 : C Normalement, c'est déjà le cas
146 : C b'=r(0) mu' = ||b'||
146 : CALL GCOPY(R,BP)
147 : MUP=RNRM2
148 : *! WRITE(IOIMP,*) 'MUP0=',MUP
149 : *! WRITE(IOIMP,*) '||XP||',GNRM2(XP)
150 : *
151 : 10 CONTINUE
152 : *
153 : * Perform Preconditionned Conjugate Gradient Squared iteration.
154 : *
155 : ITER = ITER + 1
156 : IF(IMPR.GT.7) THEN
157 : WRITE(IOIMP,*) 'ITER=',ITER
158 : ENDIF
159 : *
160 : C rho(i-1)=r[tilde]'r(i-1)
161 : RHO = GDOT(RTLD,R,IDDOT)
162 : * IF ((ITER.LE.10).OR.((ITER/100)*100.EQ.ITER)) THEN
163 : * WRITE(IOIMP,*) 'ITER=',ITER,' RHO=',RHO
164 : * ENDIF
165 : C if rho(i-1)=0 method fails
166 : IF (ABS(RHO).LT.RHOTOL) THEN
167 : RHO=SIGN(XZPREC,RHO)
168 : IWARN=IWARN+1
169 : ENDIF
170 : *
171 : * Compute vector P.
172 : *
173 : IF (ITER.GT.1) THEN
174 : C beta(i-1)=(rho(i-1)/rho(i-2))
175 : BETA = RHO / RHO1
176 : C u(i)=r(i-1)+beta(i-1)*q(i-1)
176 : CALL GCOPY(R,U)
176 : CALL GAXPY(BETA,Q,U)
177 : C p(i)=u(i)+beta(i-1)*(q(i-1)+beta(i-1)*p(i-1))
177 : CALL GSCAL(BETA**2,P)
177 : CALL GAXPY(BETA,Q,P)
177 : CALL GAXPY(ONE,U,P)
178 : ELSE
179 : C u(1)=r(0)
180 : C p(1)=u(1)
180 : CALL GCOPY(R,U)
180 : CALL GCOPY(U,P)
181 : ENDIF
182 : *
183 : * Compute direction adjusting vector PHAT and scalar ALPHA.
184 : *
185 : C
186 : C Preconditionner solve
187 : C Mp[hat]=p(i)
187 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
188 : $ ,PHAT,P)
189 : C vhat=Ap[hat]
189 : CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,PHAT,ZERO,VHAT)
190 : INMV=INMV+1
191 : C alpha(i)=rho(i-1)/(r[tilde]'vhat)
192 : TAU = GDOT(RTLD,VHAT,IDDOT)
193 : IF (ABS(TAU).LT.TAUTOL) THEN
194 : TAU=SIGN(XZPREC,TAU)
195 : IWARN=IWARN+1
196 : ENDIF
197 : ALPHA = RHO / TAU
198 : C q(i)=u(i)-alpha(i)vhat
198 : CALL GCOPY(U,Q)
198 : CALL GAXPY(-ALPHA,VHAT,Q)
199 : *
200 : * Compute direction adjusting vector UHAT.
201 : * PHAT is being used as temporary storage here.
202 : *
203 : C M uhat=u(i)+q(i)
203 : CALL GCOPY(Q,PHAT)
203 : CALL GAXPY(ONE,U,PHAT)
203 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
204 : $ ,UHAT,PHAT)
205 : *
206 : * Compute new solution approximation vector X.
207 : *
208 : * x'(i)=x'(i-1)+alpha(i)uhat
208 : CALL GAXPY(ALPHA,UHAT,XP)
209 : *
210 : * Compute residual R and check for tolerance
211 : * using Neumaier's strategy
212 : *
213 : C r(i)=b'-Ax'(i)
213 : CALL GCOPY(BP,R)
213 : CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,XP,ONE,R)
214 : INMV=INMV+1
215 : MU=GNRM2(R)
216 : RESID = MU / BNRM2
217 : CC qhat= A uhat
217 : C CALL GMOMV('N',ONE,KMORS,KISA,UHAT,ZERO,QHAT)
218 : C INMV=INMV+1
219 : CC r(i)=r(i-1)-alpha(i) qhat
219 : C CALL GAXPY(-ALPHA,QHAT,R)
220 : C RESID = GNRM2(R) / BNRM2
221 :
222 :
223 : IF(IMPR.GT.6) THEN
224 : WRITE(IOIMP,*) 'Résidu =',RESID
225 : ENDIF
226 : LRES.PROG(ITER+1)=RESID
227 : LNMV.LECT(ITER+1)=INMV
228 : IF (RESID.LE.TOL) THEN
229 : * x(i)=x(i)+x'(i)
229 : CALL GAXPY(ONE,XP,X)
230 : GOTO 30
231 : ENDIF
232 : * Neumaier's strategy
233 : IF (MU.LT.MUP) THEN
234 : * x(i)=x(i)+x'(i)
234 : CALL GAXPY(ONE,XP,X)
235 : * x'(i)=0
235 : CALL GSCAL(ZERO,XP)
236 : * b'=r(i)
236 : CALL GCOPY(R,BP)
237 : * mu'=mu
238 : MUP = MU
239 : * WRITE(IOIMP,*) 'ITER=',ITER
240 : * WRITE(IOIMP,*) 'MUP=',MUP
241 : ENDIF
242 : *
243 : * Iteration fails.
244 : *
245 : * IF (ITER.EQ.MAXIT) THEN
246 : IF (INMV.GE.MAXIT) THEN
247 : IRET = 1
248 : IF (IMPR.GT.0) THEN
249 : C WRITE(IOIMP,*)
250 : C $ 'cgsn.eso : Convergence to tolerance not achieved'
251 : C WRITE(IOIMP,*) 'INMV=',INMV,
252 : C $ ' TOL=',TOL,
253 : C $ ' RESID=',RESID
254 : IF (IWARN.GT.0) THEN
255 : WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN
256 : $ ,' times.'
257 : ENDIF
258 : ENDIF
259 : RETURN
260 : ENDIF
261 : *
262 : * Do some more
263 : *
264 : RHO1 = RHO
265 : GOTO 10
266 : *
267 : * Iteration successful
268 : *
269 : 30 CONTINUE
270 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
271 : WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN,' times.'
272 : ENDIF
273 : IRET = 0
274 : RETURN
275 : *
276 : * End of CGSN.
277 : *
278 : END
279 :
280 :
281 :
282 :
283 :
284 :
285 :
286 :
287 :
288 :
289 :
© Cast3M 2003 - All rights reserved.
Disclaimer