1 : C BICG SOURCE GOUNAND 11/10/07 21:15:09 7149
2 : SUBROUTINE BICG(KMORS,KISA,B,X,
3 : $ KPREC,INVDIA,ILUM,ILUI,
4 : $ LRES,LNMV,
5 : $ R,RTLD,Z,ZTLD,P,PTLD,Q,QTLD,
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 : BICG
11 : C DESCRIPTION :
12 : C Résolution d'un système linéaire Ax=b
13 : C par la méthode du bi-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 : POINTEUR RTLD.IZA,PTLD.IZA,QTLD.IZA,ZTLD.IZA
68 : * .. Scalar Arguments ..
69 : INTEGER ITER, KPREC, IMPR, IRET
70 : REAL*8 RESID
71 : *
72 : *
73 : * .. Variables locales
74 : * .. Parameters
75 : REAL*8 ZERO , ONE
76 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
77 : * ..
78 : INTEGER MAXIT
79 : REAL*8 ALPHA,BETA,RHO,RHO1,TOL,GNRM2,OMEGA
80 : REAL*8 GDOT,BNRM2
81 : * ..
82 : * .. External subroutines and functions..
83 : EXTERNAL GAXPY,GCOPY,GDOT,GNRM2
84 : * ..
85 : * .. Executable Statements ..
86 : *
87 :
88 : IRET = 0
89 : INMV = 0
90 : IWARN = 0
91 : IGRAN = 0
92 : MGRAN = 5
93 : MAXIT = ITER
94 : TOL = RESID
95 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bicg.eso'
96 : *
97 : * Set parameter tolerances.
98 : *
99 : RHOTOL = BRTOL
100 : OMETOL = BRTOL
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
129 : CALL GCOPY(R,RTLD)
130 : *
131 : 10 CONTINUE
132 : *
133 : * Perform (Preconditioned) BiConjugate Gradient Iteration.
134 : *
135 : ITER = ITER +1
136 : IF(IMPR.GT.7) THEN
137 : WRITE(IOIMP,*) 'ITER=',ITER
138 : ENDIF
139 : *
140 : * Preconditioner Solve
141 : C Mz(i-1)=r(i-1)
141 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
141 : *! CALL PSOLVE(KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
142 : $ ,Z,R)
143 : C M'z'(i-1)=r'(i-1)
143 : *! CALL PSOLVE(KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
143 : CALL PSOLVE('T',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
144 : $ ,ZTLD,RTLD)
145 : C rho(i-1)=z(i-1)'r'(i-1)
146 : RHO = GDOT(Z,RTLD,IDDOT)
147 : C WRITE(IOIMP,*) ' RHO=',RHO
148 : C XR1 = GDOT(Z,R)
149 : C WRITE(IOIMP,*) ' XR1=',XR1
150 : C XR2 = GDOT(ZTLD,RTLD)
151 : C WRITE(IOIMP,*) ' XR2=',XR2
152 : IF (ABS(RHO).LT.RHOTOL) THEN
153 : RHO=SIGN(XZPREC,RHO)
154 : IWARN=IWARN+1
155 : ENDIF
156 : *
157 : * Compute direction vector P
158 : *
159 : IF (ITER.GT.1) THEN
160 : C beta(i-1)=rho(i-1)/rho(i-2)
161 : BETA = RHO / RHO1
162 : C p(i)=z(i-1)+beta(i-1)p(i-1)
163 : C p'(i)=z'(i-1)+beta(i-1)p'(i-1)
163 : CALL GAXPY(BETA,P,Z)
163 : CALL GCOPY(Z,P)
163 : CALL GAXPY(BETA,PTLD,ZTLD)
163 : CALL GCOPY(ZTLD,PTLD)
164 : ELSE
165 : C p(1)=z(0)
166 : C p'(1)=z'(0)
166 : CALL GCOPY(Z,P)
166 : CALL GCOPY(ZTLD,PTLD)
167 : ENDIF
168 : *
169 : * Compute scalar ALPHA (save A*p to q)
170 : *
171 : C q(i)=Ap(i)
171 : CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,P,ZERO,Q)
172 : INMV=INMV+1
173 : C q'(i)=A'p'(i)
173 : CALL GMOMV(IMVEC,'T',ONE,KMORS,KISA,PTLD,ZERO,QTLD)
174 : INMV=INMV+1
175 : C alpha(i)=rho(i-1)/(p'(i)'q(i))
176 : OMEGA = GDOT(PTLD,Q,IDDOT)
177 : C WRITE(IOIMP,*) ' OMEGA=',OMEGA
178 : C XOM1 = GDOT(P,Q)
179 : C WRITE(IOIMP,*) ' XOM1 =',XOM1
180 : C XOM2 = GDOT(PTLD,QTLD)
181 : C WRITE(IOIMP,*) ' XOM2 =',XOM2
182 : IF (ABS(OMEGA).LT.OMETOL) THEN
183 : OMEGA=SIGN(XZPREC,OMEGA)
184 : IWARN=IWARN+1
185 : ENDIF
186 : ALPHA = RHO / OMEGA
187 : *
188 : * Compute current solution vector x
189 : *
190 : C x(i)=x(i-1)+alpha(i)p(i)
190 : CALL GAXPY(ALPHA,P,X)
191 : *
192 : * Compute residual vector R, find norm,
193 : * then check for tolerance
194 : *
195 : C r(i)=r(i-1)-alpha(i)q(i)
195 : CALL GAXPY(-ALPHA,Q,R)
196 : *
197 : RESID = GNRM2(R) / BNRM2
198 : IF(IMPR.GT.7) THEN
199 : WRITE(IOIMP,*) 'Résidu=',RESID
200 : ENDIF
201 : LRES.PROG(ITER+1)=RESID
202 : LNMV.LECT(ITER+1)=INMV
203 : IF (RESID.GT.RESMAX) THEN
204 : IGRAN=IGRAN+1
205 : ELSE
206 : IGRAN=0
207 : ENDIF
208 : *
209 : * Iteration successful
210 : *
211 : IF (RESID.LE.TOL) THEN
212 : IRET = 0
213 : *
214 : * Iteration fails
215 : *
216 : * ELSE IF (ITER.EQ.MAXIT.OR.IGRAN.GT.MGRAN) THEN
217 : ELSE IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
218 : IRET = 1
219 : *
220 : * Do some more
221 : *
222 : ELSE
223 : C r'(i)=r'(i-1)-alpha(i)q'(i)
223 : CALL GAXPY(-ALPHA,QTLD,RTLD)
224 : RHO1 = RHO
225 : GOTO 10
226 : ENDIF
227 : *
228 : * Normal termination
229 : *
230 : 30 CONTINUE
231 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
232 : WRITE(IOIMP,*) 'BiCG nearly failed ',IWARN
233 : $ ,' times.'
234 : ENDIF
235 : RETURN
236 : *
237 : * End of BICG
238 : *
239 : END
240 :
241 :
242 :
243 :
244 :
245 :
246 :
247 :
248 :
249 :
© Cast3M 2003 - All rights reserved.
Disclaimer