1 : C PRAGMG SOURCE GOUNAND 11/10/07 21:16:03 7149 2 : SUBROUTINE PRAGMG(KMORS,KTYP,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV, 3 : $ INCX,ITER,INMV, 4 : $ RESID,KPREC, 5 : $ NREST,ICALRS,IDDOT,IMVEC,KTIME,LTIME,IMPR,IRET) 6 : IMPLICIT REAL*8 (A-H,O-Z) 7 : IMPLICIT INTEGER (I-N) 8 : C*********************************************************************** 9 : C NOM : PRAGMG 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 Multigrille Algébrique 13 : C 14 : C LANGAGE : ESOPE 15 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF) 16 : C mél : gounand@semt2.smts.cea.fr 17 : C REFERENCE (bibtex-like) : 18 : C*********************************************************************** 19 : C*********************************************************************** 20 : C VERSION : v1, 17/06/08, version initiale 21 : C HISTORIQUE : 17/06/08 : Crétation 22 : C HISTORIQUE : 23 : C*********************************************************************** 24 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 25 : C en cas de modification de ce sous-programme afin de faciliter 26 : C la maintenance ! 27 : C*********************************************************************** 28 : * 29 : * .. Includes et pointeurs associés .. 30 : -INC CCOPTIO 31 : -INC CCREEL 32 : -INC SMLREEL 33 : INTEGER JG 34 : POINTEUR LRES.MLREEL 35 : -INC SMLENTI 36 : POINTEUR LNMV.MLENTI 37 : POINTEUR ATYP.MLENTI 38 : POINTEUR KTYP.MLENTI 39 : -INC SMMATRIK 40 : POINTEUR MAPREC.MATRIK 41 : POINTEUR KMORS.PMORS 42 : POINTEUR KISA.IZA 43 : POINTEUR AMORS.PMORS 44 : POINTEUR AISA.IZA 45 : POINTEUR KS2B.IZA,RES.IZA 46 : POINTEUR INCX.IZA 47 : POINTEUR INVDIA.IZA 48 : POINTEUR ILUM.PMORS 49 : POINTEUR ILUI.IZA 50 : -INC SMTABLE 51 : POINTEUR KTIME.MTABLE 52 : CHARACTER*8 CHARI 53 : LOGICAL LTIME 54 : * .. Parameters 55 : * This is a breakdown tolerance for BiCGSTAB-type method 56 : REAL*8 BRTOL 57 : *STAT-INC SMSTAT 58 : * .. Scalar Arguments .. 59 : INTEGER ITER, KPREC, IMPR, IRET 60 : INTEGER ICNT 61 : REAL*8 RESID,TOL,BNRM2,RNRM2,XFACT 62 : * .. 63 : * Vars reqd for STOPTEST2 64 : * REAL*8 TOL, BNRM2 65 : * .. 66 : * .. External subroutines .. 67 : * EXTERNAL STOPTEST2 68 : INTEGER NBVA,NJA,NTT,NTTPRE 69 : REAL*8 GNRM2 70 : EXTERNAL GNRM2 71 : * .. 72 : * .. Executable Statements .. 73 : * 74 : * WRITE(IOIMP,*) 'Debut de pragmg' 75 : IRET = 0 76 : * 77 : * On récupère les paramètres 78 : * 79 : * 80 : SEGACT MATRIK 81 : SEGINI,AMORS=KMORS 82 : NTT =AMORS.IA(/1)-1 83 : * NJA =KMORS.JA(/1) 84 : SEGINI,ATYP=KTYP 85 : SEGINI,AISA=KISA 86 : SEGACT KS2B 87 : SEGINI,RES=KS2B 88 : SEGACT INCX*MOD 89 : C 90 : C Initialisation de l'historique de convergence 91 : C 92 : JG=ITER+1 93 : SEGINI LNMV 94 : SEGINI LRES 95 : * 96 : * Autres paramètres 97 : * 98 : IJOB=0 99 : IF (IMPR.GT.2) THEN 100 : IPRINT=IOIMP 101 : ELSE 102 : IPRINT=-1 103 : ENDIF 104 : TOL=RESID 105 : ICNT=0 106 : * 107 : * AGMG ne fait que du résidu relatif d'où les petits ajustements 108 : * suivants. 109 : * 110 : * res = b - AX 110 : CALL GMOMV(IMVEC,'N',-1.D0,AMORS,AISA,INCX,1.D0,RES) 111 : BNRM2=GNRM2(KS2B) 112 : RNRM2=GNRM2(RES) 113 : IF (ICALRS.EQ.1) BNRM2=RNRM2 114 : IF (BNRM2.LT.XPETIT) BNRM2=1.D0 115 : RESID=RNRM2 / BNRM2 116 : LRES.PROG(ICNT+1)=RESID 117 : LNMV.LECT(ICNT+1)=1 118 : * 119 : IF (RESID.LE.TOL) THEN 120 : ITER=0 121 : GOTO 30 122 : ENDIF 123 : * 124 : IF (ICALRS.EQ.0) THEN 125 : XFACT=1.D0/RESID 126 : TOL=TOL*XFACT 127 : ENDIF 128 : * 129 : * Changement automatique du signe des lignes de la matrice 130 : * et du second membre si le terme diagonal est négatif. 131 : * 131 : * CALL ECMORS(AMORS,AISA,3) 132 : DO I=1,NTT 133 : IFOUND=0 134 : DO J=AMORS.IA(I),(AMORS.IA(I+1)-1) 135 : * WRITE(IOIMP,*) 'J=',J 136 : * WRITE(IOIMP,*) 'JA(J)=',AMORS.JA(J) 137 : IF (AMORS.JA(J).EQ.I) THEN 138 : IF (AISA.A(J).GT.XZERO) THEN 139 : IFOUND=1 140 : ELSEIF (AISA.A(J).LT.XZERO) THEN 141 : IFOUND=-1 142 : ENDIF 143 : GOTO 10 144 : ENDIF 145 : ENDDO 146 : 10 CONTINUE 147 : * WRITE(IOIMP,*) 'IFOUND=',IFOUND 148 : IF (IFOUND.EQ.0) THEN 149 : WRITE(IOIMP,*) 'The ',I 150 : $ ,'th diagonal term of the matrix is nil' 151 : IRET=-3 152 : GOTO 9999 153 : ELSEIF (IFOUND.EQ.-1) THEN 154 : RES.A(I)=-1.D0*RES.A(I) 155 : DO J=AMORS.IA(I),(AMORS.IA(I+1)-1) 156 : AISA.A(J)=-1.D0*AISA.A(J) 157 : ENDDO 158 : ENDIF 159 : ENDDO 160 : * WRITE(IOIMP,*) '4' 160 : * CALL ECMORS(AMORS,AISA,3) 161 : * 162 : * To use the standard AGMG library, you should make sure that it 163 : * is compiled with same compiler and options than Castem's and that 164 : * the sequential example furnished with AGMG works. 165 : * Then it is sufficient to put agmg_seq.o and agmg.o in the 166 : * current directory and use essai for linking. 167 : * 168 : * No interface to the parallel version of agmg for now 169 : * 170 : WRITE(IOIMP,*) '***********************************' 171 : WRITE(IOIMP,*) ' ' 171 : WRITE(IOIMP,*) 'Please uncomment the CALL AGMG(...', 172 : $ ' in the pragmg.eso subroutine if you wish to use' 173 : WRITE(IOIMP,*) 'the AGMG library by Y. Notay' 174 : WRITE(IOIMP,*) ' ' 175 : WRITE(IOIMP,*) '***********************************' 176 : * 251 2 177 : *Tentative d'utilisation d'une option non implémentée 177 : CALL ERREUR(251) 178 : IRET=-1 179 : GOTO 9999 180 : ********************************************************************* 181 : * 181 : * AGMG v1.1 call (standard sequential version) 182 : * 183 : C ITR2=ITER-1 183 : C CALL AGMG(NTT,AISA.A,AMORS.JA,AMORS.IA,RES.A, 184 : C $ IJOB,IPRINT,NREST,ITR2,TOL) 185 : C* X = X_0 + \delta X 185 : C CALL GAXPY(1.D0,RES,INCX) 186 : C* Résidu final 186 : C CALL GCOPY(KS2B,RES) 186 : C CALL GMOMV(IMVEC,'N',-1.D0,AMORS,AISA,INCX,1.D0,RES) 187 : C RNRM2=GNRM2(RES) 188 : C RESID=RNRM2/BNRM2 189 : C* 190 : C IF (ITR2.LT.0) THEN 191 : C IRET=1 192 : C ITER=-ITR2 193 : C ELSE 194 : C IRET=0 195 : C ITER=ITR2 196 : C ENDIF 197 : C* 198 : C ICNT=ICNT+1 199 : C LRES.PROG(ICNT+1)=RESID 200 : C LNMV.LECT(ICNT+1)=1+ITER 201 : C INMV=1+ITER 202 : * 202 : * End of AGMG call (standard sequential version) 203 : * 204 : ********************************************************************* 205 : ********************************************************************* 206 : * 206 : * AGMG call (sequential version slightly modified 207 : * for better error handling and back transmission of residual 208 : * norm history) 209 : * 210 : C ITR2=ITER-1 210 : C* CALL ECMORS(AMORS,AISA,3) 211 : C* WRITE(IOIMP,*) 'Appel de agmg' 211 : C CALL AGMG(NTT,AISA.A,AMORS.JA,AMORS.IA,ATYP.LECT,RES.A, 212 : C $ IJOB,IPRINT,NREST,ITR2,TOL,LTIME,IAT,IST,LRES.PROG,IDDOT) 213 : C* X = X_0 + \delta X 213 : C CALL GAXPY(1.D0,RES,INCX) 214 : C* 215 : C IF (ITR2.LT.0) THEN 216 : C IRET=1 217 : C ITER=-ITR2 218 : C ELSE 219 : C IRET=0 220 : C ITER=ITR2 221 : C ENDIF 222 : C IF (ICALRS.EQ.0) THEN 223 : C DO I=1,ITER 224 : C LRES.PROG(I+1)=LRES.PROG(I+1)/XFACT 225 : C ENDDO 226 : C ENDIF 227 : C DO I=1,ITER 228 : C LNMV.LECT(I+1)=I+1 229 : C ENDDO 230 : C INMV=1+ITER 231 : C RESID=LRES.PROG(ITER+1) 232 : C ICNT=ITER 233 : * 233 : * End of AGMG call (sequential version slightly modified 234 : * for better error handling and back transmission of residual 235 : * norm history) 236 : * 237 : ********************************************************************* 238 : * 239 : IF (LTIME) THEN 240 : CHARI='MGAGGREG' 240 : CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, 241 : $ 'ENTIER ',IAT,XVALR,CHARR,LOGIR,IRETR) 242 : CHARI='MGSOLUTI' 242 : CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, 243 : $ 'ENTIER ',IST,XVALR,CHARR,LOGIR,IRETR) 244 : ENDIF 245 : * 246 : C 247 : C Désactivation 248 : C 249 : 30 CONTINUE 250 : JG=ICNT+1 251 : SEGADJ LRES 252 : SEGDES LRES 253 : SEGADJ LNMV 254 : SEGDES LNMV 255 : SEGDES INCX 256 : SEGDES KS2B 257 : SEGSUP RES 258 : SEGSUP AISA 259 : SEGSUP ATYP 260 : SEGSUP AMORS 261 : SEGDES KISA 262 : SEGDES KTYP 263 : SEGDES KMORS 264 : SEGDES MATRIK 265 : C 266 : C A breakdown has occured in the CGS method 267 : C 268 : IF (IRET.LT.0) GOTO 9999 269 : * 270 : * Normal termination 271 : * 272 : RETURN 273 : * 274 : * Format handling 275 : * 276 : * 1002 FORMAT(10(1X,1PE11.4)) 277 : * 278 : * Error handling 279 : * 280 : 9999 CONTINUE 281 : WRITE(IOIMP,*) 'An error was detected in pragmg.eso' 282 : RETURN 283 : * 284 : * End of PRAGMG 285 : * 286 : END 287 : 288 : 289 : 290 : 291 : 292 : 293 : 294 : 295 : 296 : 297 : 298 : 299 :
© Cast3M 2003 - All rights reserved.
Disclaimer