1 : C KRES2 SOURCE CHAT 12/04/06 21:15:20 7345 2 : SUBROUTINE KRES2() 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 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 HISTORIQUE : 06/04/04 : Renumérotation (double mult.) 60 : C HISTORIQUE : 06/04/04 : Scaling 61 : C HISTORIQUE : 08/04/04 : ajout ILUTP 62 : C HISTORIQUE : 09/02/06 : ajout nb prod matrice vecteur (NMATVEC) 63 : C simplification du code 64 : C HISTORIQUE : 22/02/06 : Gros nettoyage au niveau de l'entrée des 65 : C arguments (Nouvelle syntaxe) 66 : C HISTORIQUE : 08/2011 : En vue de la suppression de l'objet MATRIK 67 : C on utilise l'assemblage de RESOU. 68 : C 69 : C*********************************************************************** 70 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 71 : C en cas de modification de ce sous-programme afin de faciliter 72 : C la maintenance ! 73 : C*********************************************************************** 74 : -INC CCOPTIO 75 : -INC SMLREEL 76 : POINTEUR LRES.MLREEL 77 : -INC SMLENTI 78 : POINTEUR LNMV.MLENTI 79 : -INC SMCHPOI 80 : POINTEUR KCLIM.MCHPOI 81 : POINTEUR KSMBR.MCHPOI 82 : POINTEUR MCHINI.MCHPOI 83 : POINTEUR MCHSOL.MCHPOI 84 : -INC SMTABLE 85 : POINTEUR MTINV.MTABLE 86 : POINTEUR KTIME.MTABLE 87 : DIMENSION ITTIME(6) 88 : CHARACTER*8 CHARI 89 : * MATRIK est la matrice à inverser 90 : * MAPREC est la matrice dont on utilise le préconditionnement 91 : * MATASS est la matrice dont on utilise l'assemblage 92 : * pour préconditionner celui de MATRIK 93 : -INC SMMATRIK 94 : POINTEUR MAPREC.MATRIK 95 : POINTEUR MATASS.MATRIK 96 : *STAT -INC SMSTAT 97 : C 98 : CHARACTER*8 TYPE 99 : * Paramètre m du GMRES(m) 100 : INTEGER RESTRT 101 : INTEGER ITER 102 : REAL*8 BRTOL,RESID,RXMILU,RXILUP 103 : * 104 : REAL*8 XLFIL,XDTOL 105 : INTEGER KPREC 106 : INTEGER NMATRI 107 : INTEGER IP,KTYPI 108 : INTEGER IMPINV,IIMPR 109 : CHARACTER*4 MRENU,MMULAG 110 : LOGICAL LRIG,LTIME,LDETR,LDEPE,LASRIG,LDMULT 111 : INTEGER IMPR,IRET 112 : C 113 : C Lecture des arguments et mise à défaut des optionnels () 114 : C 115 : C MATRIK : La matrice lue en entrée au format MATRIK 116 : C MTINV : L'éventuelle table d'inversion (obsolète) 117 : C IMPR : Niveau d'impression solveur direct 118 : C KCLIM : Chpoint éventuel de conditions aux limites de Dirichlet 119 : C KSMBR : Chpoint second membre 120 : C KTYPI : Type de méthode de résolution 121 : C MATASS : Matrice utilisée pour préconditionner l'assemblage 122 : C MAPREC : Matrice utilisée pour préconditionner la construction du 123 : C préconditionneur 124 : C MRENU : Type de renumérotation 125 : C MMULAG : Placement des multiplicateurs de Lagrange 126 : C ISCAL : Scaling de la matrice 127 : C IOUBL : Oubli des matrices élémentaires ? 128 : C IMPINV : Niveau d'impression solveur itératif 129 : C MCHINI : Chpoint estimation de l'inconnue 130 : C ITER : Nombre maxi d'itérations à effectuer 131 : C RESID : Norme L2 maxi du résidu 132 : C BRTOL : Breakdown tolerance pour les méthodes de type Bi-CG 133 : C IRSTRT : Paramètre m de redémarrage pour GMRES 134 : C KPREC : Type du préconditionneur 135 : C RXMILU : Paramètre de relaxation pour MILU(0) 136 : C RXILUP : Paramètre de filtre pour ILU(0)-PV 137 : C XLFIL : Paramètre de remplissage pour les préconditionneurs ILUT 138 : C XDTOL : Drop tolerance pour les préconditionneurs ILUT 139 : C XSPIV : Sensibilité du pivoting pour les préconditionneurs ILUTP 140 : C LBCG : le l dans BicgStab(l) 141 : C ICALRS : façon de calculer le résidu 142 : C METASS : Méthode d'assemblage 143 : C LTIME : construit une table avec des statistiques temporelles 144 : C si vrai 145 : C LDEPE : élimine les dépendances si VRAI 146 : C et matrice d'entrée RIGIDITE 147 : C IDDOT : 0 => utilisation du produit scalaire normal dans les 148 : C méthodes itératives 149 : C 1 => utilisation du produit scalaire compensé 150 : * IMPR=4 151 : IMPR=0 152 : * WRITE(IOIMP,*) 'coucou kres2' 153 : * 153 : *STAT CALL INMSTA(MSTAT,IMPR) 154 : * 154 : CALL PRKRES(MATRIK,MTINV,IMPR,KCLIM,KSMBR,KTYPI,MATASS,MAPREC, 155 : $ MRENU,MMULAG,ISCAL,IOUBL,IMPINV,MCHINI,ITER,RESID,BRTOL, 156 : $ IRSTRT,KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,LBCG,ICALRS, 157 : $ METASS,LTIME,LDEPE,MRIGID,IDDOT,IMVEC,IRET) 158 : IF (IRET.NE.0) GOTO 9999 159 : IMPR=MAX(IMPR,IMPINV) 160 : * 161 : * Préparation de la matrice et du second membre 162 : * suivant les cas 163 : * 164 : * LASRIG=.TRUE. on utilise l'assemblage de RESOU 165 : * LDMULT=.TRUE. on dédouble les multiplicateurs de Lagrange 166 : * 167 : 168 : LASRIG=(METASS.EQ.6) 169 : LDMULT=LASRIG 170 : MRIGI0=MRIGID 171 : IF (MRIGID.NE.0) THEN 172 : IF (LDEPE) THEN 173 : KSMBR0=KSMBR 173 : * CALL KRES6(MRIGID,KSMBR,IDEPE, 174 : * $ MRIGIC,KSMBRC,KSMBR0,KSMBR1) 174 : CALL KRES6(MRIGID,KSMBR,LDMULT, 175 : $ MRIGIC,KSMBRC,KSMBR1) 176 : IF (IERR.NE.0) RETURN 177 : MRIGID=MRIGIC 178 : KSMBR=KSMBRC 179 : ENDIF 180 : IF (LASRIG) THEN 181 : * Gestion de la normalisation 182 : NORICO=NORINC 183 : NORIDO=NORIND 184 : IF (ISCAL.EQ.0) THEN 185 : NORINC=0 186 : NORIND=0 187 : ELSE 188 : NORINC=-1 189 : NORIND=0 190 : ENDIF 191 : * Gestion de la renumérotation 192 : NUCROO=NUCROU 193 : IF (MRENU.EQ.'RIEN') THEN 194 : NUCROU=-1 195 : ELSEIF (MRENU.EQ.'SLOA') THEN 196 : NUCROU=1 197 : ELSEIF (MRENU.EQ.'GIPR'.OR.MRENU.EQ.'GIBA') THEN 198 : NUCROU=2 199 : ELSE 200 : WRITE(IOIMP,*) 'MRENU=',MRENU 200 : CALL ERREUR(5) 201 : RETURN 202 : ENDIF 202 : CALL KRES8(MRIGID,KSMBR, 203 : $ KTYPI,ITER,RESID,ICALRS,IRSTRT,LBCG,BRTOL,IDDOT,IMVEC, 204 : $ KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV, 205 : $ KTIME,LTIME, 206 : $ MCHSOL,LRES,LNMV,ICVG,IMPR) 207 : IF (IERR.NE.0) RETURN 208 : * Gestion de la normalisation 209 : NORINC=NORICO 210 : NORIND=NORIDO 211 : NUCROU=NUCROO 211 : IF (LTIME) CALL ECROBJ('TABLE ',KTIME) 212 : IF (MTINV.NE.0) THEN 212 : CALL ECMO(MTINV,'CVGOK','ENTIER',ICVG) 212 : IF (LRES.NE.0) CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES) 212 : IF (LNMV.NE.0) CALL ECMO(MTINV,'NMATVEC','LISTENTI',LNMV) 213 : ENDIF 214 : * On décondense si nécessaire 215 : * 216 : * write (6,*) ' resou - mchsol ' 216 : * call ecchpo(mchsol,0) 216 : * call mucpri(mchsol,mrigid,iresi) 217 : * write (6,*) ' kres - iresi ' 217 : * call ecchpo(iresi,0) 218 : * WRITE(IOIMP,*) 'Avant KRES7' 219 : IF (MRIGI0.NE.0.AND.LDEPE) THEN 220 : MSOLC=MCHSOL 220 : CALL KRES7(MSOLC,MRIGI0,KSMBR0,KSMBR1,1, 221 : $ MCHSOL) 222 : IF (IERR.NE.0) RETURN 223 : ENDIF 223 : CALL ECROBJ('CHPOINT ',MCHSOL) 224 : RETURN 225 : ELSE 225 : CALL ECROBJ('RIGIDITE',MRIGID) 225 : CALL RIMA 226 : IF (IERR.NE.0) GOTO 9999 226 : CALL MACHI2(1) 227 : IF (IERR.NE.0) GOTO 9999 227 : CALL LIROBJ('MATRIK',MATRIK,1,IRET) 228 : IF(IRET.EQ.0) GOTO 9999 229 : 230 : * Changement des noms d'inconnues du second membre 231 : IF (KSMBR.NE.0) THEN 231 : CALL ECROBJ('CHPOINT ',KSMBR) 231 : CALL MACHI2(1) 231 : CALL LIROBJ('CHPOINT ',KSMBR,1,IRET) 232 : IF (IRET.EQ.0) GOTO 9999 233 : ENDIF 234 : ENDIF 235 : * write (6,*) ' le vecteur 2' 235 : * call ecchpo(ksmbr,0) 236 : * write (6,*) ' la matrice 2' 236 : * call ecrent(5) 236 : * call ecmatk(matrik) 237 : ENDIF 238 : * 239 : SEGACT MATRIK 240 : NMATRI=IRIGEL(/2) 241 : IF(NMATRI.EQ.0)THEN 242 : C% Résolution impossible : la matrice de RIGIDITE est vide 242 : CALL ERREUR(727) 243 : RETURN 244 : ENDIF 245 : SEGDES MATRIK 246 : IF (MATASS.EQ.0) MATASS=MATRIK 247 : IF (MAPREC.EQ.0) MAPREC=MATRIK 248 : * WRITE(IOIMP,*) 'Sortie de prkres' 249 : * WRITE(IOIMP,*) 'IOUBL=',IOUBL 250 : C 251 : IF (LTIME) THEN 251 : CALL CRTABL(KTIME) 251 : CALL TIMESPV(ITTIME) 252 : ITI1=(ITTIME(1)+ITTIME(2)) 253 : ELSE 254 : KTIME=0 255 : ENDIF 255 : *STAT CALL PRMSTA('Lectures',MSTAT,IMPR) 256 : * 257 : C 258 : C Assemblage proprement dit 259 : C 260 : IIMPR=0 260 : CALL KRES3(MATRIK,MATASS,MRENU,MMULAG,METASS, 261 : $ KTIME,LTIME, 262 : $ IIMPR,IRET) 263 : IF (IRET.NE.0) GOTO 9999 264 : *! WRITE(IOIMP,*) 'Aprés assemblage' 264 : *! CALL ECRENT(5) 264 : *! CALL ECROBJ('MATRIK ',MATRIK) 264 : *! CALL PRLIST 265 : IF (LTIME) THEN 265 : CALL TIMESPV(ITTIME) 266 : ITI2=(ITTIME(1)+ITTIME(2)) 267 : ENDIF 267 : *STAT CALL PRMSTA('Assemblage',MSTAT,IMPR) 268 : * 269 : * "Oubli" des valeurs des matrice élémentaires 270 : * On met les tableaux de LIZAFM à 0 => à MENAGE de les supprimmer 271 : * si besoin est. 272 : * 273 : IOUBD=MOD(IOUBL,10) 274 : *! WRITE(IOIMP,*) 'IOUBD=',IOUBD 275 : IF (IOUBD.EQ.1) THEN 275 : CALL OUBIMA(MATRIK,IMPR,IRET) 276 : IF (IRET.NE.0) GOTO 9999 277 : IF (IMPR.GT.2) THEN 278 : WRITE(IOIMP,*) 'Oubli des mat. elem.' 279 : ENDIF 280 : ELSEIF (IOUBD.EQ.2) THEN 281 : ith=0 281 : call ooonth(ith) 281 : call ooohor(0,ith) 282 : SEGACT MATRIK*MOD 283 : LDETR=.FALSE. 284 : NMATE=IRIGEL(/2) 285 : DO IMATE=1,NMATE 286 : IMATRI=IRIGEL(4,IMATE) 287 : SEGACT IMATRI*MOD 288 : NSOUM =LIZAFM(/1) 289 : NTOTIN=LIZAFM(/2) 290 : DO ITOTIN=1,NTOTIN 291 : DO ISOUM=1,NSOUM 292 : IZAFM=LIZAFM(ISOUM,ITOTIN) 293 : IF (IZAFM.NE.0) THEN 294 : LDETR=.TRUE. 295 : SEGSUP IZAFM 296 : LIZAFM(ISOUM,ITOTIN)=0 297 : ENDIF 298 : ENDDO 299 : ENDDO 300 : SEGDES IMATRI 301 : ENDDO 302 : IF (IMPR.GT.2.AND.LDETR) THEN 303 : WRITE(IOIMP,*) 'Destruction des mat. elem.' 304 : ENDIF 305 : ELSEIF (IOUBD.NE.0) THEN 306 : WRITE(IOIMP,*) 'IOUBL=',IOUBL, ' non prevu' 307 : GOTO 9999 308 : ENDIF 308 : *STAT CALL PRMSTA('Oubli',MSTAT,IMPR) 309 : *! WRITE(IOIMP,*) 'Aprés oubli' 310 : C 311 : C Méthode directe 312 : C 313 : IF (KTYPI.EQ.1) THEN 313 : CALL KRES4(MATRIK,KCLIM,KSMBR, 314 : $ ISCAL, 315 : $ MCHSOL, 316 : $ IMPR,IRET) 317 : IF (IRET.NE.0) GOTO 9999 317 : *STAT CALL PRMSTA('Methode directe',MSTAT,IMPR) 318 : C 319 : C Méthodes itératives 320 : C 321 : ELSE 322 : C 322 : CALL KRES5(MATRIK,KCLIM,KSMBR,KTYPI, 323 : $ MCHINI,ITER,RESID, 324 : $ BRTOL,IRSTRT,LBCG,ICALRS, 325 : $ MAPREC,KPREC, 326 : $ RXMILU,RXILUP,XLFIL,XDTOL,XSPIV, 327 : $ ISCAL, 328 : $ KTIME,LTIME,IDDOT,IMVEC, 329 : $ MCHSOL,LRES,LNMV,ICVG, 330 : $ IMPR,IRET) 331 : IF (IRET.NE.0) GOTO 9999 331 : *STAT CALL PRMSTA('Methode iterative',MSTAT,IMPR) 332 : IF (MTINV.NE.0) THEN 332 : CALL ECMO(MTINV,'CVGOK','ENTIER',ICVG) 332 : CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES) 332 : CALL ECMO(MTINV,'NMATVEC','LISTENTI',LNMV) 333 : ENDIF 334 : ENDIF 335 : IF (LTIME) THEN 335 : CALL TIMESPV(ITTIME) 336 : ITI3=(ITTIME(1)+ITTIME(2)) 337 : CHARI='ASS+RENU' 337 : CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, 338 : $ 'ENTIER ',ITI2-ITI1,XVALR,CHARR,LOGIR,IRETR) 339 : CHARI='PRE+RESO' 339 : CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, 340 : $ 'ENTIER ',ITI3-ITI2,XVALR,CHARR,LOGIR,IRETR) 341 : CHARI='TOTAL ' 341 : CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, 342 : $ 'ENTIER ',ITI3-ITI1,XVALR,CHARR,LOGIR,IRETR) 343 : SEGDES KTIME 343 : CALL ECROBJ('TABLE ',KTIME) 344 : ENDIF 345 : IOUBE=IOUBL/10 346 : *! WRITE(IOIMP,*) 'IOUBE=',IOUBE 347 : IF (IOUBE.GE.1) THEN 348 : ith=0 348 : call ooonth(ith) 348 : call ooohor(0,ith) 349 : SEGACT MATRIK*MOD 350 : IF (IOUBE.EQ.2) THEN 351 : PMORS=KIDMAT(4) 352 : IF (PMORS.NE.0) THEN 353 : IF (IMPR.GT.2) THEN 354 : WRITE(IOIMP,*) 'Destruction du profil morse' 355 : ENDIF 356 : SEGSUP PMORS 357 : KIDMAT(4)=0 358 : ENDIF 359 : ENDIF 360 : IZA=KIDMAT(5) 361 : IF (IZA.NE.0) THEN 362 : IF (IMPR.GT.2) THEN 363 : WRITE(IOIMP,*) 'Destruction des valeurs morses' 364 : ENDIF 365 : SEGSUP IZA 366 : KIDMAT(5)=0 367 : ENDIF 368 : PMORS=KIDMAT(6) 369 : IF (PMORS.NE.0) THEN 370 : IF (IMPR.GT.2) THEN 371 : WRITE(IOIMP,*) 'Destruction du profil du precon' 372 : ENDIF 373 : SEGSUP PMORS 374 : KIDMAT(6)=0 375 : ENDIF 376 : IZA=KIDMAT(7) 377 : IF (IZA.NE.0) THEN 378 : IF (IMPR.GT.2) THEN 379 : WRITE(IOIMP,*) 'Destruction des valeurs du precon' 380 : ENDIF 381 : SEGSUP IZA 382 : KIDMAT(7)=0 383 : ENDIF 384 : SEGDES MATRIK 385 : ELSEIF (IOUBE.NE.0) THEN 386 : WRITE(IOIMP,*) 'IOUBL=',IOUBL, ' non prevu' 387 : GOTO 9999 388 : ENDIF 389 : * 390 : * On décondense si nécessaire 391 : * 392 : * write (6,*) ' resou - mchsol ' 392 : * call ecchpo(mchsol,0) 392 : * call mucpri(mchsol,mrigid,iresi) 393 : * write (6,*) ' kres - iresi ' 393 : * call ecchpo(iresi,0) 394 : * WRITE(IOIMP,*) 'Avant KRES7' 395 : IF (MRIGI0.NE.0.AND.LDEPE) THEN 396 : MSOLC=MCHSOL 396 : * CALL KRES7(MRIGID,IDEPE,KSMBR0,KSMBR1, 397 : * $ MSOLC, 398 : * $ MCHSOL) 398 : CALL KRES7(MSOLC,MRIGI0,KSMBR0,KSMBR1,1, 399 : $ MCHSOL) 400 : IF (IERR.NE.0) RETURN 401 : ENDIF 401 : CALL ECROBJ('CHPOINT ',MCHSOL) 401 : *STAT CALL SUMSTA(MSTAT,IMPR) 402 : * 403 : * Normal termination 404 : * 405 : RETURN 406 : * 407 : * Format handling 408 : * 409 : * 410 : * Error handling 411 : * 412 : 9999 CONTINUE 413 : WRITE(IOIMP,*) 'An error was detected in kres2.eso' 414 : * 153 2 415 : * Opération illicite dans ce contexte 415 : CALL ERREUR(153) 416 : RETURN 417 : * 418 : * End of KRES2 419 : * 420 : END 421 : 422 : 423 : 424 : 425 : 426 : 427 : 428 : 429 : 430 : 431 : 432 : 433 : 434 : 435 : 436 : 437 : 438 : 439 : 440 : 441 : 442 : 443 : 444 : 445 : 446 :
© Cast3M 2003 - All rights reserved.
Disclaimer