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