1 : C PRCGS SOURCE GOUNAND 11/07/21 21:15:40 7046
2 : SUBROUTINE PRCGS(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,
3 : $ INCX,ITER,INMV,
4 : $ RESID,KPREC,
5 : $ BRTOL,ICALRS,IMPR,IRET)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : IMPLICIT INTEGER (I-N)
8 : C***********************************************************************
9 : C NOM : PRCGS
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)
13 : C préconditionnée ou non.
14 : C
15 : C LANGAGE : ESOPE
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***********************************************************************
29 : C VERSION : v1, 10/02/06, version initiale
30 : C HISTORIQUE : 10/02/06 : Crétation
31 : C HISTORIQUE :
32 : C***********************************************************************
33 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
34 : C en cas de modification de ce sous-programme afin de faciliter
35 : C la maintenance !
36 : C***********************************************************************
37 : *
38 : * .. Includes et pointeurs associés ..
39 : -INC CCOPTIO
40 : -INC SMLREEL
41 : INTEGER JG
42 : POINTEUR LRES.MLREEL
43 : -INC SMLENTI
44 : POINTEUR LNMV.MLENTI
45 : -INC SMMATRIK
46 : POINTEUR MAPREC.MATRIK
47 : POINTEUR KMORS.PMORS
48 : POINTEUR KISA.IZA
49 : POINTEUR KS2B.IZA
50 : POINTEUR INCX.IZA
51 : POINTEUR INVDIA.IZA
52 : POINTEUR ILUM.PMORS
53 : POINTEUR ILUI.IZA
54 : * .. Parameters
55 : * This is a breakdown tolerance for BiCGSTAB-type method
56 : REAL*8 BRTOL
57 : * .. Work Vectors for CGS..
58 : POINTEUR IR.IZA,IRTLD.IZA,IP.IZA,IPHAT.IZA
59 : POINTEUR IQ.IZA,IQHAT.IZA,IUHAT.IZA
60 : *STAT-INC SMSTAT
61 : * .. Scalar Arguments ..
62 : INTEGER ITER, KPREC, IMPR, IRET
63 : REAL*8 RESID
64 : * ..
65 : * Vars reqd for STOPTEST2
66 : * REAL*8 TOL, BNRM2
67 : * ..
68 : * .. External subroutines ..
69 : * EXTERNAL STOPTEST2
70 : INTEGER NBVA,NJA,NTT,NTTPRE
71 : * ..
72 : * .. Executable Statements ..
73 : *
74 : IRET = 0
75 : *
76 : * On récupère les paramètres
77 : *
78 :
79 : SEGACT MATRIK
80 : SEGACT MAPREC
81 : IF (KSYM.EQ.0) THEN
82 : IF (IMPR.GT.0) THEN
83 : WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
84 : $ 'Use a Conjugate Gradient instead !'
85 : C IRET=-2
86 : C GOTO 9999
87 : ENDIF
88 : ENDIF
89 : * Pour le préconditionneur
90 : ILUM =MAPREC.KIDMAT(6)
91 : ILUI =MAPREC.KIDMAT(7)
92 : IDMAT=MAPREC.KIDMAT(1)
93 : SEGACT IDMAT
94 : INVDIA=IDIAG
95 : SEGDES IDMAT
96 : SEGACT KMORS
97 : NTT =KMORS.IA(/1)-1
98 : * NJA =KMORS.JA(/1)
99 : SEGACT KISA
100 : SEGACT KS2B
101 : SEGACT INCX*MOD
102 : IF (KPREC.NE.0) THEN
103 : C Paramètres des préconditionnements diagonaux et D-ILU
104 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
105 : C Est-il compatible avec la matrice ?
106 : SEGACT INVDIA
107 : NTTPRE=INVDIA.A(/1)
108 : IF (NTTPRE.NE.NTT) THEN
109 : WRITE(IOIMP,*) 'La matrice et son préconditionnement'
110 : WRITE(IOIMP,*) 'ne sont pas compatibles...'
111 : WRITE(IOIMP,*) 'NTT=',NTT
112 : WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
113 : IRET=-2
114 : GOTO 9999
115 : ENDIF
116 : C Paramètres des préconditionnements ILU(0), MILU(0), ILUT, ILUT2
117 : C ,ilutp, ilutpg, ilutpg2
118 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
119 : SEGACT ILUM
120 : SEGACT ILUI
121 : * SEGPRT,ILUM
122 : * SEGPRT,ILUI
123 : NTTPRE=ILUM.IA(/1)
124 : IF (NTTPRE.NE.NTT) THEN
125 : WRITE(IOIMP,*) 'La matrice et son préconditionnement',
126 : $ 'ne sont pas compatibles...'
127 : WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
128 : IRET=-2
129 : GOTO 9999
130 : ENDIF
131 : ELSE
132 : INVDIA=0
133 : WRITE(IOIMP,*) 'Erreur de programmation'
134 : GOTO 9999
135 : ENDIF
136 : ENDIF
137 : C
138 : C Initialisation de l'historique de convergence
139 : C
140 : JG=ITER+1
141 : SEGINI LNMV
142 : SEGINI LRES
143 : C
144 : C
145 : C Initialisation des vecteurs de travail pour BiCGSTAB
146 : C
146 : *STAT CALL INMSTA(MSTAT,0)
147 : NBVA=NTT
148 : SEGINI IR,IRTLD,IP,IPHAT
149 : SEGINI IQ,IQHAT,IUHAT
149 : *STAT CALL PRMSTA(' PRCGS : objets temporaires',MSTAT,1)
149 : *STAT CALL SUMSTA(MSTAT,0)
150 : C
151 : C Appel de CGS
152 : C
152 : CALL CGS(KMORS,KISA,KS2B,INCX,
153 : $ KPREC,INVDIA,ILUM,ILUI,
154 : $ LRES,LNMV,
155 : $ IR,IRTLD,IP,IPHAT,
156 : $ IQ,IQHAT,IUHAT,
157 : $ ITER,INMV,BRTOL,RESID,ICALRS,IMPR,IRET)
158 : C
159 : C Désactivation
160 : C
161 : SEGSUP IR,IRTLD,IP,IPHAT
162 : SEGSUP IQ,IQHAT,IUHAT
163 : JG=ITER+1
164 : SEGADJ LRES
165 : SEGDES LRES
166 : SEGADJ LNMV
167 : SEGDES LNMV
168 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
169 : SEGDES INVDIA
170 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
171 : SEGDES ILUI
172 : SEGDES ILUM
173 : ENDIF
174 : SEGDES INCX
175 : SEGDES KS2B
176 : SEGDES KISA
177 : SEGDES KMORS
178 : SEGDES MAPREC
179 : SEGDES MATRIK
180 : C
181 : C A breakdown has occured in the CGS method
182 : C
183 : IF (IRET.LT.0) GOTO 9999
184 : *
185 : * Normal termination
186 : *
187 : RETURN
188 : *
189 : * Format handling
190 : *
191 : * 1002 FORMAT(10(1X,1PE11.4))
192 : *
193 : * Error handling
194 : *
195 : 9999 CONTINUE
196 : WRITE(IOIMP,*) 'An error was detected in prcgs.eso'
197 : RETURN
198 : *
199 : * End of PRCGS
200 : *
201 : END
202 :
203 :
204 :
205 :
206 :
207 :
208 :
209 :
210 :
211 :
212 :
© Cast3M 2003 - All rights reserved.
Disclaimer