1 : C PRBCGF SOURCE GOUNAND 11/07/21 21:15:34 7046
2 : SUBROUTINE PRBCGF(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,INCX,
3 : $ ITER,INMV,RESID,KPREC,
4 : $ BRTOL,LBCG,ICALRS,IMPR,IRET)
5 : IMPLICIT INTEGER(I-N)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : C***********************************************************************
8 : C NOM : PRBCGF
9 : C DESCRIPTION :
10 : C Préparation de la résolution d'un système linéaire Ax=b
11 : C par une méthode BiCGSTAB(l) préconditionnée ou non.
12 : C (+reliable updates + convex combination)
13 : C LANGAGE : ESOPE
14 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
15 : C mél : gounand@semt2.smts.cea.fr
16 : C REFERENCE (bibtex-like) :
17 : C @BOOK{templates,
18 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
19 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
20 : C H. Van der Vorst},
21 : C TITLE={Templates for the Solution of Linear Systems :
22 : C Building Blocks for Iterative Methods},
23 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
24 : C -> URL : http://www.netlib.org/templates/Templates.html
25 : C@TechReport{fokkema,
26 : C author = {DR Fokkema},
27 : C title = {Enhanced implementation of BiCGSTAB(l) for solving
28 : C linear systems of equations},
29 : C institution = {?},
30 : C year = {1996}}
31 : C***********************************************************************
32 : C VERSION : v1, 22/02/06, version initiale
33 : C HISTORIQUE : v1, 22/02/06, création
34 : C HISTORIQUE :
35 : C HISTORIQUE :
36 : C***********************************************************************
37 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
38 : C en cas de modification de ce sous-programme afin de faciliter
39 : C la maintenance !
40 : C***********************************************************************
41 : *
42 : * .. Includes et pointeurs associés ..
43 : -INC CCOPTIO
44 : -INC SMLREEL
45 : INTEGER JG
46 : POINTEUR LRES.MLREEL
47 : -INC SMLENTI
48 : POINTEUR LNMV.MLENTI
49 : POINTEUR ITRAV.MLENTI
50 : -INC SMMATRIK
51 : POINTEUR MAPREC.MATRIK
52 : POINTEUR KMORS.PMORS
53 : POINTEUR KISA.IZA
54 : POINTEUR KS2B.IZA
55 : POINTEUR INCX.IZA
56 : POINTEUR INVDIA.IZA
57 : POINTEUR ILUM.PMORS
58 : POINTEUR ILUI.IZA
59 : * .. Work Vectors for BiCGSTAB(l)
60 : SEGMENT SPACE
61 : REAL*8 IJ(NI,NJ)
62 : ENDSEGMENT
63 : POINTEUR IZ.SPACE,IZ2.SPACE,IZ2M1.SPACE
64 : SEGMENT SPACE2
65 : POINTEUR WRK(NIZA).IZA
66 : ENDSEGMENT
67 : POINTEUR IR.SPACE2,IU.SPACE2
68 : POINTEUR IRTLD.IZA,IXTLD.IZA,IUHAT.IZA,IBP.IZA
69 : POINTEUR IY0.IZA,IYL.IZA
70 : POINTEUR ITMP.IZA
71 : * .. Scalar Arguments ..
72 : INTEGER ITER, KPREC, IMPR, IRET
73 : REAL*8 RESID
74 : INTEGER NBVA,NJA,NTT,NTTPRE
75 : * .. Executable Statements ..
76 : *
77 : IRET = 0
78 : *
79 : * On récupère les paramètres
80 : *
81 : SEGACT MATRIK
82 : SEGACT MAPREC
83 : IF (KSYM.EQ.0) THEN
84 : IF (IMPR.GT.2) THEN
85 : WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
86 : $ 'use a Conjugate Gradient instead !'
87 : ENDIF
88 : C IRET=-2
89 : C GOTO 9999
90 : ENDIF
91 : * Pour le préconditionneur
92 : ILUM =MAPREC.KIDMAT(6)
93 : ILUI =MAPREC.KIDMAT(7)
94 : IDMAT=MAPREC.KIDMAT(1)
95 : SEGACT IDMAT
96 : INVDIA=IDIAG
97 : SEGDES IDMAT
98 :
99 : SEGACT KMORS
100 : NTT =KMORS.IA(/1)-1
101 : * NJA =KMORS.JA(/1)
102 : SEGACT KISA
103 : SEGACT KS2B
104 : SEGACT INCX*MOD
105 : C Paramètres des préconditionnements diagonaux et D-ILU
106 : IF (KPREC.NE.0) THEN
107 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
108 : C Est-il compatible avec la matrice ?
109 : SEGACT INVDIA
110 : NTTPRE=INVDIA.A(/1)
111 : IF (NTTPRE.NE.NTT) THEN
112 : WRITE(IOIMP,*) 'La matrice et son préconditionnement'
113 : WRITE(IOIMP,*) 'ne sont pas compatibles...'
114 : WRITE(IOIMP,*) 'NTT=',NTT
115 : WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
116 : IRET=-2
117 : GOTO 9999
118 : ENDIF
119 : C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
120 : C ilutp, ilutpg, ilutpg2
121 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
122 : SEGACT ILUM
123 : SEGACT ILUI
124 : NTTPRE=ILUM.IA(/1)
125 : IF (NTTPRE.NE.NTT) THEN
126 : WRITE(IOIMP,*) 'La matrice et son préconditionnement',
127 : $ 'ne sont pas compatibles...'
128 : WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
129 : IRET=-2
130 : GOTO 9999
131 : ENDIF
132 : ENDIF
133 : ENDIF
134 : C
135 : C Initialisation de l'historique de convergence
136 : C
137 : JG=ITER+1
138 : SEGINI LNMV
139 : SEGINI LRES
140 : C
141 : C
142 : C Initialisation des vecteurs de travail pour BiCGStab(l)
143 : C
144 : NI=LBCG+1
145 : NJ=LBCG+1
146 : SEGINI IZ
147 : NI=LBCG-1
148 : NJ=LBCG-1
149 : SEGINI IZ2,IZ2M1
150 : NBVA=LBCG+1
151 : SEGINI IY0,IYL
152 : JG=LBCG-1
153 : SEGINI ITRAV
154 : NBVA=NTT
155 : SEGINI IRTLD,IXTLD,IUHAT,IBP
156 : NIZA=LBCG+1
157 : SEGINI IR
158 : DO I=1,NIZA
159 : SEGINI ITMP
160 : IR.WRK(I)=ITMP
161 : ENDDO
162 : SEGINI IU
163 : DO I=1,NIZA
164 : SEGINI ITMP
165 : IU.WRK(I)=ITMP
166 : ENDDO
167 : C
167 : CALL BCGF(KMORS,KISA,KS2B,INCX,
168 : $ KPREC,INVDIA,ILUM,ILUI,
169 : $ LRES,LNMV,
170 : $ IRTLD,IXTLD,IUHAT,IR,IU,IZ,IZ2,IZ2M1,IY0,IYL,IBP,ITRAV,
171 : $ ITER,INMV,BRTOL,LBCG,RESID,ICALRS,IMPR,IRET)
172 : C
173 : C Désactivation-suppression
174 : C
175 : DO I=1,NIZA
176 : ITMP=IU.WRK(I)
177 : SEGSUP ITMP
178 : ENDDO
179 : SEGSUP IU
180 : DO I=1,NIZA
181 : ITMP=IR.WRK(I)
182 : SEGSUP ITMP
183 : ENDDO
184 : SEGSUP IR
185 : SEGSUP IRTLD,IXTLD,IUHAT,IBP
186 : SEGSUP ITRAV
187 : SEGSUP IY0,IYL
188 : SEGSUP IZ2,IZ2M1
189 : SEGSUP IZ
190 : JG=ITER+1
191 : SEGADJ LRES
192 : SEGDES LRES
193 : SEGADJ LNMV
194 : SEGDES LNMV
195 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
196 : SEGDES INVDIA
197 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
198 : SEGDES ILUI
199 : SEGDES ILUM
200 : ENDIF
201 : SEGDES INCX
202 : SEGDES KS2B
203 : SEGDES KISA
204 : SEGDES KMORS
205 : SEGDES MAPREC
206 : SEGDES MATRIK
207 : C
208 : C A problem has occured in the GMRES method
209 : C
210 : IF (IRET.LT.0) GOTO 9999
211 : *
212 : * Normal termination
213 : *
214 : RETURN
215 : *
216 : * Format handling
217 : *
218 : *
219 : * Error handling
220 : *
221 : 9999 CONTINUE
222 : WRITE(IOIMP,*) 'An error was detected in prbcgf.eso'
223 : RETURN
224 : *
225 : * End of prbcgf
226 : *
227 : END
228 :
229 :
230 :
231 :
232 :
233 :
234 :
235 :
236 :
237 :
238 :
239 :
240 :
© Cast3M 2003 - All rights reserved.
Disclaimer