1 : C PRBCGG SOURCE GOUNAND 11/10/07 21:16:06 7149
2 : SUBROUTINE PRBCGG(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,INCX,
3 : $ ITER,INMV,RESID,KPREC,
4 : $ BRTOL,LBCG,ICALRS,IDDOT,IMVEC,IMPR,IRET)
5 : IMPLICIT INTEGER(I-N)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : C***********************************************************************
8 : C NOM : PRBCGG
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)
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,IZM1.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 IY.IZA,IYP.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
145 : NJ=LBCG
146 : SEGINI IZ,IZM1
147 : NBVA=LBCG
148 : SEGINI IY,IYP
149 : JG=LBCG
150 : SEGINI ITRAV
151 : NBVA=NTT
152 : SEGINI IRTLD,IXTLD,IUHAT,IBP
153 : NIZA=LBCG+1
154 : SEGINI IR
155 : DO I=1,NIZA
156 : SEGINI ITMP
157 : IR.WRK(I)=ITMP
158 : ENDDO
159 : SEGINI IU
160 : DO I=1,NIZA
161 : SEGINI ITMP
162 : IU.WRK(I)=ITMP
163 : ENDDO
164 : C
164 : CALL BCGG(KMORS,KISA,KS2B,INCX,
165 : $ KPREC,INVDIA,ILUM,ILUI,
166 : $ LRES,LNMV,
167 : $ IRTLD,IXTLD,IUHAT,IR,IU,IZ,IZM1,IY,IYP,IBP,ITRAV,
168 : $ ITER,INMV,BRTOL,LBCG,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
169 : C
170 : C Désactivation-suppression
171 : C
172 : DO I=1,NIZA
173 : ITMP=IU.WRK(I)
174 : SEGSUP ITMP
175 : ENDDO
176 : SEGSUP IU
177 : DO I=1,NIZA
178 : ITMP=IR.WRK(I)
179 : SEGSUP ITMP
180 : ENDDO
181 : SEGSUP IR
182 : SEGSUP IRTLD,IXTLD,IUHAT,IBP
183 : SEGSUP ITRAV
184 : SEGSUP IY,IYP
185 : SEGSUP IZ,IZM1
186 : JG=ITER+1
187 : SEGADJ LRES
188 : SEGDES LRES
189 : SEGADJ LNMV
190 : SEGDES LNMV
191 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
192 : SEGDES INVDIA
193 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
194 : SEGDES ILUI
195 : SEGDES ILUM
196 : ENDIF
197 : SEGDES INCX
198 : SEGDES KS2B
199 : SEGDES KISA
200 : SEGDES KMORS
201 : SEGDES MAPREC
202 : SEGDES MATRIK
203 : C
204 : C A problem has occured in the GMRES method
205 : C
206 : IF (IRET.LT.0) GOTO 9999
207 : *
208 : * Normal termination
209 : *
210 : RETURN
211 : *
212 : * Format handling
213 : *
214 : *
215 : * Error handling
216 : *
217 : 9999 CONTINUE
218 : WRITE(IOIMP,*) 'An error was detected in prbcgg.eso'
219 : RETURN
220 : *
221 : * End of prbcgg
222 : *
223 : END
224 :
225 :
226 :
227 :
228 :
229 :
230 :
231 :
232 :
233 :
234 :
235 :
236 :
237 :
© Cast3M 2003 - All rights reserved.
Disclaimer