1 : C PRCGSN SOURCE GOUNAND 11/10/07 21:16:08 7149
2 : SUBROUTINE PRCGSN(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,
3 : $ INCX,ITER,INMV,
4 : $ RESID,KPREC,
5 : $ BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : IMPLICIT INTEGER (I-N)
8 : C***********************************************************************
9 : C NOM : PRCGSN
10 : C DESCRIPTION :
11 : C Préparation de la résolution d'un système linéaire Ax=b
12 : C par une méthode de gradient conjugué au carré (CGS) avec la
13 : C stratégie de stabilisation de Neumaier préconditionnée ou non.
14 : C On stocke deux vecteurs en plus de ceux de CGS normal
15 : C
16 : C LANGAGE : ESOPE
17 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
18 : C mél : gounand@semt2.smts.cea.fr
19 : C REFERENCE (bibtex-like) :
20 : C @BOOK{templates,
21 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
22 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
23 : C H. Van der Vorst},
24 : C TITLE={Templates for the Solution of Linear Systems :
25 : C Building Blocks for Iterative Methods},
26 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
27 : C -> URL : http://www.netlib.org/templates/Templates.html
28 : C pour CGS
29 : C
30 : C@TechReport{vorstsleipjen,
31 : C author = {GLG Sleijpen and HA. Van der Vorst},
32 : C title = {Hybrid Bi-Conjugate Gradient Methods for CFD Problems},
33 : C institution = {Universiteit Utrecht},
34 : C year = {1995},
35 : C type = {Preprint},
36 : C number = {902}}
37 : C pour la stratégie de Neumaier
38 : C***********************************************************************
39 : C***********************************************************************
40 : C VERSION : v1, 10/02/06, version initiale
41 : C HISTORIQUE : 10/02/06 : Crétation
42 : C HISTORIQUE :
43 : C***********************************************************************
44 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
45 : C en cas de modification de ce sous-programme afin de faciliter
46 : C la maintenance !
47 : C***********************************************************************
48 : *
49 : * .. Includes et pointeurs associés ..
50 : -INC CCOPTIO
51 : -INC SMLREEL
52 : INTEGER JG
53 : POINTEUR LRES.MLREEL
54 : -INC SMLENTI
55 : POINTEUR LNMV.MLENTI
56 : -INC SMMATRIK
57 : POINTEUR MAPREC.MATRIK
58 : POINTEUR KMORS.PMORS
59 : POINTEUR KISA.IZA
60 : POINTEUR KS2B.IZA
61 : POINTEUR INCX.IZA
62 : POINTEUR INVDIA.IZA
63 : POINTEUR ILUM.PMORS
64 : POINTEUR ILUI.IZA
65 : * .. Parameters
66 : * This is a breakdown tolerance for BiCGSTAB-type method
67 : REAL*8 BRTOL
68 : * .. Work Vectors for CGS..
69 : POINTEUR IR.IZA,IRTLD.IZA,IP.IZA,IPHAT.IZA
70 : POINTEUR IQ.IZA,IQHAT.IZA,IUHAT.IZA
71 : POINTEUR IXP.IZA,IBP.IZA
72 : *STAT-INC SMSTAT
73 : * .. Scalar Arguments ..
74 : INTEGER ITER, KPREC, IMPR, IRET
75 : REAL*8 RESID
76 : * ..
77 : * Vars reqd for STOPTEST2
78 : * REAL*8 TOL, BNRM2
79 : * ..
80 : * .. External subroutines ..
81 : * EXTERNAL STOPTEST2
82 : INTEGER NBVA,NJA,NTT,NTTPRE
83 : * ..
84 : * .. Executable Statements ..
85 : *
86 : IRET = 0
87 : *
88 : * On récupère les paramètres
89 : *
90 :
91 : SEGACT MATRIK
92 : SEGACT MAPREC
93 : IF (KSYM.EQ.0) THEN
94 : IF (IMPR.GT.0) THEN
95 : WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
96 : $ 'Use a Conjugate Gradient instead !'
97 : C IRET=-2
98 : C GOTO 9999
99 : ENDIF
100 : ENDIF
101 : * Pour le préconditionneur
102 : ILUM =MAPREC.KIDMAT(6)
103 : ILUI =MAPREC.KIDMAT(7)
104 : IDMAT=MAPREC.KIDMAT(1)
105 : SEGACT IDMAT
106 : INVDIA=IDIAG
107 : SEGDES IDMAT
108 : SEGACT KMORS
109 : NTT =KMORS.IA(/1)-1
110 : * NJA =KMORS.JA(/1)
111 : SEGACT KISA
112 : SEGACT KS2B
113 : SEGACT INCX*MOD
114 : IF (KPREC.NE.0) THEN
115 : C Paramètres des préconditionnements diagonaux et D-ILU
116 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
117 : C Est-il compatible avec la matrice ?
118 : SEGACT INVDIA
119 : NTTPRE=INVDIA.A(/1)
120 : IF (NTTPRE.NE.NTT) THEN
121 : WRITE(IOIMP,*) 'La matrice et son préconditionnement'
122 : WRITE(IOIMP,*) 'ne sont pas compatibles...'
123 : WRITE(IOIMP,*) 'NTT=',NTT
124 : WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
125 : IRET=-2
126 : GOTO 9999
127 : ENDIF
128 : C Paramètres des préconditionnements ILU(0), MILU(0), ILUT, ILUT2
129 : C ,ilutp, ilutpg, ilutpg2
130 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
131 : SEGACT ILUM
132 : SEGACT ILUI
133 : * SEGPRT,ILUM
134 : * SEGPRT,ILUI
135 : NTTPRE=ILUM.IA(/1)
136 : IF (NTTPRE.NE.NTT) THEN
137 : WRITE(IOIMP,*) 'La matrice et son préconditionnement',
138 : $ 'ne sont pas compatibles...'
139 : WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
140 : IRET=-2
141 : GOTO 9999
142 : ENDIF
143 : ELSE
144 : INVDIA=0
145 : WRITE(IOIMP,*) 'Erreur de programmation'
146 : GOTO 9999
147 : ENDIF
148 : ENDIF
149 : C
150 : C Initialisation de l'historique de convergence
151 : C
152 : JG=ITER+1
153 : SEGINI LNMV
154 : SEGINI LRES
155 : C
156 : C
157 : C Initialisation des vecteurs de travail pour BiCGSTAB
158 : C
158 : *STAT CALL INMSTA(MSTAT,0)
159 : NBVA=NTT
160 : SEGINI IR,IRTLD,IP,IPHAT
161 : SEGINI IQ,IQHAT,IUHAT
162 : SEGINI IXP,IBP
162 : *STAT CALL PRMSTA(' PRCGSN : objets temporaires',MSTAT,1)
162 : *STAT CALL SUMSTA(MSTAT,0)
163 : C
164 : C Appel de CGSN
165 : C
165 : CALL CGSN(KMORS,KISA,KS2B,INCX,
166 : $ KPREC,INVDIA,ILUM,ILUI,
167 : $ LRES,LNMV,
168 : $ IR,IRTLD,IP,IPHAT,
169 : $ IQ,IQHAT,IUHAT,
170 : $ IXP,IBP,
171 : $ ITER,INMV,BRTOL,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
172 : C
173 : C Désactivation
174 : C
175 : SEGSUP IXP,IBP
176 : SEGSUP IR,IRTLD,IP,IPHAT
177 : SEGSUP IQ,IQHAT,IUHAT
178 : JG=ITER+1
179 : SEGADJ LRES
180 : SEGDES LRES
181 : SEGADJ LNMV
182 : SEGDES LNMV
183 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
184 : SEGDES INVDIA
185 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
186 : SEGDES ILUI
187 : SEGDES ILUM
188 : ENDIF
189 : SEGDES INCX
190 : SEGDES KS2B
191 : SEGDES KISA
192 : SEGDES KMORS
193 : SEGDES MAPREC
194 : SEGDES MATRIK
195 : C
196 : C A breakdown has occured in the CGS method
197 : C
198 : IF (IRET.LT.0) GOTO 9999
199 : *
200 : * Normal termination
201 : *
202 : RETURN
203 : *
204 : * Format handling
205 : *
206 : * 1002 FORMAT(10(1X,1PE11.4))
207 : *
208 : * Error handling
209 : *
210 : 9999 CONTINUE
211 : WRITE(IOIMP,*) 'An error was detected in prcgsn.eso'
212 : RETURN
213 : *
214 : * End of PRCGSN
215 : *
216 : END
217 :
218 :
219 :
220 :
221 :
222 :
223 :
224 :
225 :
226 :
227 :
228 :
© Cast3M 2003 - All rights reserved.
Disclaimer