1 : C KRESLL SOURCE GOUNAND 11/10/07 21:15:57 7149 2 : SUBROUTINE KRESLL() 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C************************************************************************* 6 : C Operateur KRES 7 : C 8 : C OBJET : Resoud l'equation AX=B par différentes méthodes 9 : C * directe (factorisation Choleski) 10 : C * itérative (Gradient conjugué, BiCGSTAB, 11 : C GMRES(m)) avec différents préconditionneurs 12 : C diagonal (Jacobi), D-ILU, ILU(0) (Choleski 13 : C incomplet), MILU(0) 14 : C SYNTAXE : CHPO3 = KRES MA1 'TYPI' TAB1 15 : C 'CLIM' CHPO1 16 : C 'SMBR' CHPO2 17 : C 'IMPR' VAL1 ; 18 : C Voir la notice pour plus de précisions. 19 : C 20 : C*********************************************************************** 21 : C APPELES : KRES3, KRES4, KRES5 22 : C APPELES (E/S) : LIROBJ, ECROBJ, ERREUR, LIRMOT, LIRENT, LIRTAB, 23 : C ACME, ACMO, ACMM, ACMF, ECMO 24 : C APPELES (UTIL.) : QUETYP 25 : C APPELE PAR : KOPS 26 : C*********************************************************************** 27 : C*********************************************************************** 28 : C HISTORIQUE : 26/10/98 : prise en compte du cas particulier 29 : C où A est diagonale. C'est en fait assez pénible 30 : C car on utilise le CHPOINT comme structure de 31 : C données pour la matrice A et les vecteurs X,B,CLIM 32 : C HISTORIQUE : 09/02/99 : on peut utiliser le préconditionnement d'une 33 : C autre matrice que celle que l'on inverse 34 : C HISTORIQUE : 01/06/99 : on modifie la partie résolution itérative 35 : C en effet, lors de l'imposition des CL. de 36 : C Dirichlet, on changeait les valeurs de la 37 : C matrice Morse. 38 : C Ceci n'est pas bon lorsqu'on veut utiliser la 39 : C matrice assemblée pour plusieurs pas de temps. 40 : C On travaille maintenant sur une copie. 41 : C HISTORIQUE : 20/12/99 : reprogrammation de l'assemblage 42 : C Le nouvel assemblage n'est, pour l'instant effectif que 43 : C lorsqu'une méthode itérative est sélectionnée (-> fabrication 44 : C d'une matrice Morse). Le nouvel assemblage est plus performant 45 : C (temps de calcul, utilisation de la mémoire) et plus général (cas 46 : C particuliers à peu près supprimés) que le précédent. 47 : C HISTORIQUE : 05/01/00 : On ne supprime plus les 0.D0 de la matrice 48 : C assemblée (cf. clmors appelé par melim). Ceci pour ne plus avoir 49 : C perte éventuelle de symétrie de la matrice assemblée. Cela permet 50 : C aussi de ne plus dupliquer la matrice assemblée. 51 : C HISTORIQUE : 13/01/00 : Mise en conformité du solveur direct avec le 52 : C nouvel assemblage. 53 : C HISTORIQUE : 22/03/00 : Rajout des préconditionneurs ILUT 54 : C HISTORIQUE : 13/04/00 : Séparation Lecture des données 55 : C Ecriture des résultats (kres2) 56 : C Assemblage kres3 57 : C Méthode directe kres4 58 : C Méthodes itératives kres5 59 : C 60 : C*********************************************************************** 61 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 62 : C en cas de modification de ce sous-programme afin de faciliter 63 : C la maintenance ! 64 : C*********************************************************************** 65 : -INC CCOPTIO 66 : -INC SMLREEL 67 : POINTEUR LRES.MLREEL 68 : -INC SMCHPOI 69 : POINTEUR KCLIM.MCHPOI 70 : POINTEUR KSMBR.MCHPOI 71 : POINTEUR MCHINI.MCHPOI 72 : POINTEUR MCHSOL.MCHPOI 73 : -INC SMTABLE 74 : POINTEUR MTINV.MTABLE 75 : * MATRIK est la matrice à inverser 76 : * MAPREC est la matrice dont on utilise le préconditionnement 77 : * MATASS est la matrice dont on utilise l'assemblage 78 : * pour préconditionner celui de MATRIK 79 : -INC SMMATRIK 80 : POINTEUR MAPREC.MATRIK 81 : POINTEUR MATASS.MATRIK 82 : C 83 : CHARACTER*8 TYPE 84 : * Paramètre m du GMRES(m) 85 : INTEGER IRSTRT 86 : INTEGER ITER 87 : REAL*8 BRTOL,RESID,RXMILU 88 : * 89 : REAL*8 XLFIL,XDTOL 90 : INTEGER KPREC 91 : INTEGER NMATRI 92 : INTEGER IP,KTYPI 93 : INTEGER IMPINV,IIMPR 94 : CHARACTER*4 MRENU,MMULAG 95 : INTEGER IMPR,IRET 96 : C 97 : C Définition des options 98 : C 99 : C 'IMPR' + un entier : niveau d'impression 100 : C 'TYPI' + une table contenant les informations 101 : C sur le type de résolution 102 : C 'CLIM' + un chpoint : conditions aux limites de Dirichlet 103 : C 'SMBR' + un chpoint : second membre (forces...) 104 : C 'PREC' + une matrice dont on va utiliser le préconditionnement 105 : C 106 : INTEGER NBM 107 : PARAMETER (NBM=4) 108 : CHARACTER*4 LMOT(NBM) 109 : C Tableau de correspondance : entier <-> mot 110 : C pour le type d'inversion 111 : INTEGER NBTYPI 112 : C INTEGER LNTYPI 113 : PARAMETER (NBTYPI=5) 114 : C PARAMETER (LNTYPI=18) 115 : C CHARACTER*(LNTYPI) LTYPI(NBTYPI) 116 : LOGICAL LBID,LVAL,LTIME,LMRENU,LDEPE,LLDEPE 117 : C 118 : C Initialisation des tableaux 119 : C 120 : DATA LMOT/'IMPR','TYPI','CLIM','SMBR'/ 121 : C DATA LTYPI/'Choleski ', 122 : C $ 'Conjugate Gradient', 123 : C $ 'BiCGSTAB ', 124 : C $ 'BiCGSTAB(2) ', 125 : C $ 'GMRES(m) '/ 126 : C 127 : TYPE='MATRIK' 127 : CALL LIROBJ(TYPE,MATRIK,1,IRET) 128 : IF(IRET.EQ.0) GOTO 9999 129 : c 130 : SEGACT MATRIK 131 : NMATRI=IRIGEL(/2) 132 : IF(NMATRI.EQ.0)THEN 133 : C% Résolution impossible : la matrice de RIGIDITE est vide 133 : CALL ERREUR(727) 134 : RETURN 135 : ENDIF 136 : SEGDES MATRIK 137 : c Valeurs par défaut reprise de prkres 138 : KSMBR=0 139 : MTINV=0 140 : IMPR=0 141 : KCLIM=0 142 : KTYPI=1 143 : MATASS=MATRIK 144 : MAPREC=MATRIK 145 : MRENU='SLOA' 146 : LMRENU=.FALSE. 147 : MMULAG='APR2' 148 : ISCAL=0 149 : IOUBL=0 150 : IMPINV=0 151 : MCHINI=0 152 : ITER=2000 153 : RESID=1.D-10 154 : BRTOL=1.D-40 155 : IRSTRT=50 156 : KPREC=3 157 : RXMILU=1.D0 158 : RXILUP=0.5D0 159 : XLFIL=2.D0 160 : XDTOL=-1.D0 161 : XSPIV=0.1D0 162 : LBCG=4 163 : ICALRS=0 164 : METASS=5 165 : LTIME=.FALSE. 166 : KTIME=0 167 : LDEPE=.TRUE. 168 : LLDEPE=.FALSE. 169 : IDDOT=0 170 : IMVEC=2 171 : C 172 : 1 CONTINUE 172 : CALL LIRMOT(LMOT,NBM,IP,0) 173 : IF(IP.EQ.0) GOTO 2 174 : GOTO (11,12,13,14),IP 175 : C Lecture du niveau d'impression 176 : C option 'IMPR' 177 : 11 CONTINUE 178 : c write(6,*)' option IMPR' 178 : CALL LIRENT(IMPR,1,IRET) 179 : IF (IRET.EQ.0) GOTO 9999 180 : GO TO 1 181 : C Lecture de la table pour les méthodes d'inversion 182 : C option 'TYPI' 183 : 12 CONTINUE 183 : CALL LIRTAB('METHINV',MTINV,1,IRET) 184 : IF(IRET.EQ.0)RETURN 184 : CALL LIRENT(KTYPI,1,IRET) 185 : C KTYPI : Méthode d'inversion du système (cf. LTYPI) 186 : C 1 : résolution directe (Choleski) 187 : C 2 : Gradient Conjugué 188 : C 3 : Bi-Gradient Conjugué Stabilisé (BiCGSTAB) 189 : C 4 : BiCGSTAB(2) 190 : C 5 : GMRES(m) 190 : C CALL ACME(MTINV,'TYPINV',KTYPI) 191 : IF (IERR.NE.0) GOTO 9999 192 : C write(6,*)' Ktypi=',ktypi 193 : IF (KTYPI.LT.1.OR.KTYPI.GT.NBTYPI) THEN 194 : WRITE(IOIMP,*) 'L''indice TYPINV ', 195 : $ 'is out of range (1..',NBTYPI,') :',KTYPI 196 : GOTO 9999 197 : ENDIF 198 : GO TO 1 199 : C Lecture du chpoint de conditions aux limites 200 : C option CLIM 201 : 13 CONTINUE 201 : CALL QUETYP(TYPE,1,IRET) 202 : IF(IRET.EQ.0) GOTO 9999 203 : IF(TYPE.EQ.'CHPOINT')THEN 203 : CALL LIROBJ(TYPE,KCLIM,1,IRET) 204 : ELSE 205 : WRITE(IOIMP,*) 'Erreur lecture C.lim' 206 : GOTO 9999 207 : ENDIF 208 : GO TO 1 209 : C Lecture du chpoint second membre 210 : C option SMBR 211 : 14 CONTINUE 211 : CALL QUETYP(TYPE,1,IRET) 212 : IF(IRET.EQ.0)RETURN 213 : IF(TYPE.EQ.'CHPOINT')THEN 213 : CALL LIROBJ('CHPOINT',KSMBR,1,IRET) 214 : c write(6,*)' lecture 2 membre ' 215 : ELSE 216 : WRITE(IOIMP,*) 'Erreur lecture second membre' 217 : GOTO 9999 218 : ENDIF 219 : GOTO 1 220 : 2 CONTINUE 221 : C write(6,*)' On a pas finit les lectures ' 222 : c NAT=2 223 : c NSOUPO=0 224 : c SEGINI MCHPO1 225 : c MCHPO1.JATTRI(1)=2 225 : c CALL ECROBJ('CHPOINT',MCHPO1) 226 : c return 227 : C 228 : C Matrice dont on utilise l'assemblage 229 : C (éventuellement identique à la matrice transmise) 230 : C 231 : TYPE='MATRIK ' 231 : C CALL ACMO(MTINV,'MATASS',TYPE,MATASS) 232 : C IF (IERR.NE.0) GOTO 9999 232 : CALL LIROBJ('MATRIK',MATASS,1,IRET) 233 : IF(IRET.EQ.0) GOTO 9999 234 : C - Pour l'assemblage : type de renumérotation 234 : C CALL ACMM(MTINV,'TYRENU',MRENU) 235 : C IF (IERR.NE.0) GOTO 9999 235 : CALL LIRCHA(MRENU,1,LCHAR) 236 : IF(LCHAR.EQ.0)GO TO 9999 237 : C - Pour l'assemblage : prise en compte des mult.lag 237 : C CALL ACMM(MTINV,'PCMLAG',MMULAG) 238 : C IF (IERR.NE.0) GOTO 9999 238 : CALL LIRCHA(MMULAG,1,LCHAR) 239 : IF(LCHAR.EQ.0)GO TO 9999 240 : C - Pour l'assemblage : scaling de la matrice 240 : C CALL ACME(MTINV,'SCALING',ISCAL) 241 : C IF (IERR.NE.0) GOTO 9999 242 : ISCAL=0 243 : C Niveau d'impression pour la partie résolution itérative 243 : c CALL ACME(MTINV,'IMPINV',IMPINV) 244 : c IF (IERR.NE.0) GOTO 9999 245 : c IMPR=MAX(IMPR,IMPINV) 246 : IMPR=0 247 : C 248 : C Assemblage proprement dit 249 : C 250 : IIMPR=0 251 : METASS=4 251 : CALL KRES3(MATRIK,MATASS,MRENU,MMULAG,METASS, 252 : $ KTIME,LTIME, 253 : $ IIMPR,IRET) 254 : IF (IRET.NE.0) GOTO 9999 255 : *! WRITE(IOIMP,*) 'Aprés assemblage' 255 : *! CALL ECRENT(5) 255 : *! CALL ECROBJ('MATRIK ',MATRIK) 255 : *! CALL PRLIST 256 : * 257 : * "Oubli" des valeurs des matrice élémentaires 258 : * On met les tableaux de LIZAFM à 0 => à MENAGE de les supprimmer 259 : * si besoin est. 260 : * 260 : CALL OUBIMA(MATRIK,IMPR,IRET) 261 : IF (IRET.NE.0) GOTO 9999 262 : C 263 : C Méthode directe 264 : C 265 : IF (KTYPI.EQ.1) THEN 265 : CALL KRES4(MATRIK,KCLIM,KSMBR, 266 : $ ISCAL, 267 : $ MCHSOL, 268 : $ IMPR,IRET) 269 : IF (IRET.NE.0) GOTO 9999 270 : C 271 : C Méthodes itératives 272 : C 273 : ELSEIF(KTYPI.GE.2.AND.KTYPI.LE.NBTYPI)THEN 274 : C Paramètres pour la méthode itérative proprement dite 275 : C - Champoint d'initialisation de la méthode 276 : C (i.e. estimation de l'inconnue) 277 : TYPE='CHPOINT ' 277 : C CALL ACMO(MTINV,'XINIT',TYPE,MCHINI) 278 : C IF (IERR.NE.0) GOTO 9999 278 : CALL LIROBJ('CHPOINT',MCHINI,1,IRET) 279 : IF(IRET.EQ.0) GOTO 9999 280 : C - Nombre maxi d'itérations à effectuer 280 : C CALL ACME(MTINV,'NITMAX',ITER) 281 : C IF (IERR.NE.0) GOTO 9999 281 : CALL LIRENT(ITER,1,IRET) 282 : IF(IRET.EQ.0) GOTO 9999 283 : C - Norme maxi (L2 normé par le second membre) du résidu 283 : C CALL ACMF(MTINV,'RESID',RESID) 284 : C IF (IERR.NE.0) GOTO 9999 284 : CALL LIRREE(RESID,1,IRET) 285 : IF(IRET.EQ.0) GOTO 9999 286 : C - 'Breakdown tolerance' pour les méthodes de type 287 : C BiCGSTAB. Si un certain produit scalaire de vecteurs 288 : C "direction" est inférieur à cette tolérance, la 289 : C méthode s'arrete. 290 : IF (KTYPI.EQ.3.OR.KTYPI.EQ.4) THEN 290 : C CALL ACMF(MTINV,'BCGSBTOL',BRTOL) 291 : C IF (IERR.NE.0) GOTO 9999 291 : CALL LIRREE(BRTOL,1,IRET) 292 : IF(IRET.EQ.0) GOTO 9999 293 : ENDIF 294 : C - Paramètre m de redémarrage pour GMRES(m) 295 : C i.e. nombre max. de vecteurs par rapport auxquels on 296 : C orthogonalise le résidu 297 : IF (KTYPI.EQ.5) THEN 297 : C CALL ACME(MTINV,'GMRESTRT',IRSTRT) 298 : C IF (IERR.NE.0) GOTO 9999 298 : CALL LIRENT(IRSTRT,1,IRET) 299 : IF(IRET.EQ.0) GOTO 9999 300 : ENDIF 301 : C Paramètres pour les différents préconditionneurs 302 : C Matrice dont on utilise le préconditionneur 303 : C (éventuellement identique à la matrice transmise). 304 : C TYPE='MATRIK ' 304 : C CALL ACMO(MTINV,'MAPREC',TYPE,MAPREC) 305 : C IF (IERR.NE.0) GOTO 9999 305 : CALL LIROBJ('MATRIK',MAPREC,1,IRET) 306 : IF(IRET.EQ.0) GOTO 9999 307 : C - Type de préconditionnement : (cf. LPREC) 308 : C 0 : pas de préconditionnement 309 : C 1 : préconditionnement par la diagonale 310 : C 2 : préconditionnement D-ILU 311 : C 3 : préconditionnement ILU(0) (Choleski) 312 : C 4 : préconditionnement MILU(0) (Choleski modifié) 313 : C 5 : préconditionnement ILUT (dual truncation) 314 : C 6 : préconditionnement ILUT2 (une variante du 315 : C précédent qui remplit mieux la mémoire et 316 : C fonctionne mieux quelquefois) 316 : C CALL ACME(MTINV,'PRECOND',KPREC) 317 : C IF (IERR.NE.0) GOTO 9999 317 : CALL LIRENT(KPREC,1,IRET) 318 : IF(IRET.EQ.0) GOTO 9999 319 : C - Paramètre de relaxation pour le préconditionnement 320 : C MILU(0) compris entre 0.D0 et 1.D0 321 : C S'il est égal à 0, on se ramène à ILU(0) 322 : C S'il est égal à 1, MILU(0) est dit non relaxé 323 : IF (KPREC.EQ.4) THEN 323 : C CALL ACMF(MTINV,'MILURELX',RXMILU) 324 : C IF (IERR.NE.0) GOTO 9999 324 : CALL LIRREE(RXMILU,1,IRET) 325 : IF(IRET.EQ.0) GOTO 9999 326 : ENDIF 327 : C - Pour une méthode ILUT, on a les deux indices suivant : 328 : C * ILUTLFIL : encombrement maximal (approximatif) du 329 : C préconditionneur, par rapport à la matrice. 330 : C * ILUTDTOL : "drop tolerance" pour le préconditionneur. 331 : C i.e. en-dessous de cette valeur relative, les 332 : C termes de la factorisation incomplète seront 333 : C oubliés. 334 : IF (KPREC.GE.5.AND.KPREC.LE.9) THEN 334 : C CALL ACMF(MTINV,'ILUTLFIL',XLFIL) 335 : C IF (IERR.NE.0) GOTO 9999 335 : CALL LIRREE(XLFIL,1,IRET) 336 : IF(IRET.EQ.0) GOTO 9999 336 : C CALL ACMF(MTINV,'ILUTDTOL',XDTOL) 337 : C IF (IERR.NE.0) GOTO 9999 337 : CALL LIRREE(XDTOL,1,IRET) 338 : IF(IRET.EQ.0) GOTO 9999 339 : ENDIF 340 : C - Pour une méthode ILUTP, on a les deux indices suivant : 341 : C * ILUTPPIV (type REEL) : sensibilité du pivoting 342 : C IF (KPREC.GE.7.AND.KPREC.LE.9) THEN 342 : C CALL ACMF(MTINV,'ILUTPPIV',XSPIV) 343 : C ENDIF 343 : CALL KRES5(MATRIK,KCLIM,KSMBR,KTYPI, 344 : $ MCHINI,ITER,RESID, 345 : $ BRTOL,IRSTRT,LBCG,ICALRS, 346 : $ MAPREC,KPREC, 347 : $ RXMILU,RXILUP,XLFIL,XDTOL,XSPIV, 348 : $ ISCAL, 349 : $ KTIME,LTIME,IDDOT,IMVEC, 350 : $ MCHSOL,LRES,LNMV,ICVG, 351 : $ IMPR,IRET) 352 : IF (IRET.NE.0) GOTO 9999 352 : CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES) 353 : ELSE 354 : WRITE(IOIMP,*) 'KTYPI=',KTYPI,'out of range (1..',NBTYPI,')' 355 : GOTO 9999 356 : ENDIF 356 : CALL ECROBJ('CHPOINT ',MCHSOL) 357 : * 358 : * Normal termination 359 : * 360 : RETURN 361 : * 362 : * Format handling 363 : * 364 : * 365 : * Error handling 366 : * 367 : 9999 CONTINUE 368 : WRITE(IOIMP,*) 'An error was detected in kresll.eso' 369 : * 153 2 370 : * Opération illicite dans ce contexte 370 : CALL ERREUR(153) 371 : RETURN 372 : * 373 : * End of KRES2 374 : * 375 : END 376 : 377 : 378 : 379 : 380 : 381 : 382 : 383 : 384 : 385 : 386 : 387 : 388 : 389 : 390 : 391 : 392 : 393 : 394 : 395 : 396 : 397 :
© Cast3M 2003 - All rights reserved.
Disclaimer