C KRES2 SOURCE GOUNAND 25/11/12 21:15:34 12399 SUBROUTINE KRES2() IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C************************************************************************* C Operateur KRES C C OBJET : Resoud l'equation AX=B par différentes méthodes C * directe (factorisation Choleski) C * itérative (Gradient conjugué, BiCGSTAB, C GMRES(m)) avec différents préconditionneurs C diagonal (Jacobi), D-ILU, ILU(0) (Choleski C incomplet), MILU(0) C SYNTAXE : CHPO3 = KRES MA1 'TYPI' TAB1 C 'CLIM' CHPO1 C 'SMBR' CHPO2 C 'IMPR' VAL1 ; C Voir la notice pour plus de précisions. C C*********************************************************************** C APPELES : KRES3, KRES4, KRES5 C APPELES (E/S) : LIROBJ, ECROBJ, ERREUR, LIRMOT, LIRENT, LIRTAB, C ACME, ACMO, ACMM, ACMF, ECMO, ECME C APPELES (UTIL.) : QUETYP C APPELE PAR : KOPS C*********************************************************************** C*********************************************************************** C HISTORIQUE : 26/10/98 : prise en compte du cas particulier C où A est diagonale. C'est en fait assez pénible C car on utilise le CHPOINT comme structure de C données pour la matrice A et les vecteurs X,B,CLIM C HISTORIQUE : 09/02/99 : on peut utiliser le préconditionnement d'une C autre matrice que celle que l'on inverse C HISTORIQUE : 01/06/99 : on modifie la partie résolution itérative C en effet, lors de l'imposition des CL. de C Dirichlet, on changeait les valeurs de la C matrice Morse. C Ceci n'est pas bon lorsqu'on veut utiliser la C matrice assemblée pour plusieurs pas de temps. C On travaille maintenant sur une copie. C HISTORIQUE : 20/12/99 : reprogrammation de l'assemblage C Le nouvel assemblage n'est, pour l'instant effectif que C lorsqu'une méthode itérative est sélectionnée (-> fabrication C d'une matrice Morse). Le nouvel assemblage est plus performant C (temps de calcul, utilisation de la mémoire) et plus général (cas C particuliers à peu près supprimés) que le précédent. C HISTORIQUE : 05/01/00 : On ne supprime plus les 0.D0 de la matrice C assemblée (cf. clmors appelé par melim). Ceci pour ne plus avoir C perte éventuelle de symétrie de la matrice assemblée. Cela permet C aussi de ne plus dupliquer la matrice assemblée. C HISTORIQUE : 13/01/00 : Mise en conformité du solveur direct avec le C nouvel assemblage. C HISTORIQUE : 22/03/00 : Rajout des préconditionneurs ILUT C HISTORIQUE : 13/04/00 : Séparation Lecture des données C Ecriture des résultats (kres2) C Assemblage kres3 C Méthode directe kres4 C Méthodes itératives kres5 C HISTORIQUE : 06/04/04 : Renumérotation (double mult.) C HISTORIQUE : 06/04/04 : Scaling C HISTORIQUE : 08/04/04 : ajout ILUTP C HISTORIQUE : 09/02/06 : ajout nb prod matrice vecteur (NMATVEC) C simplification du code C HISTORIQUE : 22/02/06 : Gros nettoyage au niveau de l'entrée des C arguments (Nouvelle syntaxe) C HISTORIQUE : 08/2011 : En vue de la suppression de l'objet MATRIK C on utilise l'assemblage de RESOU. C HISTORIQUE : 04/2019 : remplacement de NOEL par NELIM C Idéalement, il faudrait reprendre ce que Pierre a fait dans C RESOU dans les fiches 10[0,1]?? et avec MREM.En vue de la C suppression de l'objet MATRIK C HISTORIQUE : 04/2025 : Grand menage pour la numerotation et AGMG : C Elimination recursive via resouc et resour C Utilisation systematique assemblage via RESOU C Placement des multiplicateurs type pression revus C AGMG Stokes et Navier-Stokes C En vue de la suppression de l'objet MATRIK C C*********************************************************************** C Prière de PRENDRE LE TEMPS de compléter les commentaires C en cas de modification de ce sous-programme afin de faciliter C la maintenance ! C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC SMLREEL POINTEUR LRES.MLREEL -INC SMLENTI POINTEUR LNMV.MLENTI -INC SMCHPOI POINTEUR KCLIM.MCHPOI POINTEUR KSMBR.MCHPOI POINTEUR MCHINI.MCHPOI POINTEUR MCHSOL.MCHPOI -INC SMTABLE POINTEUR MTINV.MTABLE POINTEUR KTIME.MTABLE DIMENSION ITTIME(4) CHARACTER*16 CHARI * MATRIK est la matrice à inverser * MAPREC est la matrice dont on utilise le préconditionnement * MATASS est la matrice dont on utilise l'assemblage * pour préconditionner celui de MATRIK -INC SMMATRIK POINTEUR MAPREC.MATRIK POINTEUR MATASS.MATRIK * STAT -INC SMSTAT PARAMETER (NELMAX=30) SEGMENT IDEMEM(0) segment ideme0(idemem(/1),NELMAX) segment ideme1(idemem(/1),NELMAX) C CHARACTER*8 TYPE,CHARR * Paramètre m du GMRES(m) INTEGER RESTRT INTEGER ITER REAL*8 BRTOL,RESID,RXMILU,RXILUP * REAL*8 XLFIL,XDTOL INTEGER KPREC INTEGER NMATRI INTEGER IP,KTYPI INTEGER IMPINV,IIMPR CHARACTER*4 MRENU,MMULAG LOGICAL LRIG,LTIME,LDETR,LASRIG,LOGII,LDUMP,LOGIR INTEGER IMPR,IRET C C Lecture des arguments et mise à défaut des optionnels () C * IMPR=4 C Lecture des arguments et mise à défaut des optionnels () C C MATRIK : La matrice lue en entrée au format MATRIK C MTINV : L'éventuelle table d'inversion (obsolète) C IMPR : Niveau d'impression solveur direct C KCLIM : Chpoint éventuel de conditions aux limites de Dirichlet C KSMBR : Chpoint second membre C KTYPI : Type de méthode de résolution C MATASS : Matrice utilisée pour préconditionner l'assemblage C MAPREC : Matrice utilisée pour préconditionner la construction du C préconditionneur C MRENU : Type de renumérotation C MMULAG : Placement des multiplicateurs de Lagrange C ISCAL : Scaling de la matrice C INORMU : Mise à l'échelle des multiplicateurs de Lagrange C IOUBL : Oubli des matrices élémentaires ? C IMPINV : Niveau d'impression solveur itératif C MCHINI : Chpoint estimation de l'inconnue C ITER : Nombre maxi d'itérations à effectuer C RESID : Norme L2 maxi du résidu C BRTOL : Breakdown tolerance pour les méthodes de type Bi-CG C IRSTRT : Paramètre m de redémarrage pour GMRES C KPREC : Type du préconditionneur C RXMILU : Paramètre de relaxation pour MILU(0) C RXILUP : Paramètre de filtre pour ILU(0)-PV C XLFIL : Paramètre de remplissage pour les préconditionneurs ILUT C XDTOL : Drop tolerance pour les préconditionneurs ILUT C XSPIV : Sensibilité du pivoting pour les préconditionneurs ILUTP C LBCG : le l dans BicgStab(l) C ICALRS : façon de calculer le résidu C METASS : méthode d'assemblage C LTIME : construit une table avec des statistiques temporelles C si vrai C NELIM : nombre de passes d'elimination C Par defaut : 2 C IDDOT : 0 => utilisation du produit scalaire normal dans les C méthodes itératives C 1 => utilisation du produit scalaire compensé C IMVEC : 0, pas de parallélisme pour les produits matrice-vecteur C 1, parallélisme stratégie 1, entrelace les lignes. C 2, parallélisme stratégie 2, groupe les lignes. C Par defaut : 2 C IORINC : pointeur sur une LISTE de LISTMOTS indiquant quels noms C d'inconnues pour chaque bloc pour AGMG par blocs C Par defaut : 0 C MLAG1 : pointeur sur un LISTMOTS indiquant les noms d'inconnues a C traiter comme des multiplicateurs de Lagrange, à placer après les C inconnues en relation C Par defaut : 0 C MLAG2 : pointeur sur un LISTMOTS indiquant les noms d'inconnues a C traiter comme des multiplicateurs de Lagrange, à placer avant les C inconnues en relation C Par defaut : 0 IVALI=0 XVALI=0.D0 LOGII=.FALSE. IRETI=0 XVALR=0.D0 *inutile IOBRE=0 IRETR=0 IMPR=0 * WRITE(IOIMP,*) 'coucou kres2' * *STAT CALL INMSTA(MSTAT,IMPR) * CALL PRKRES(MATRIK,MTINV,IMPR,KCLIM,KSMBR,KTYPI,MATASS,MAPREC, $ MRENU,MMULAG,ISCAL,INORMU,IOUBL,IMPINV,MCHINI,ITER,RESID $ ,BRTOL,IRSTRT,KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,LBCG $ ,ICALRS,METASS,LTIME,NELIM,MRIGID,IDDOT,IMVEC,IORINC,MLAG1 $ ,MLAG2,LDUMP,ISMOOT,IRET) IF (IRET.NE.0) GOTO 9999 IMPR=MAX(IMPR,IMPINV) * * Préparation de la matrice et du second membre * suivant les cas * * LASRIG=.TRUE. on utilise l'assemblage de RESOU * LASRIG=(METASS.EQ.6) LASRIG=.TRUE. * Pour l'instant, il faut toujours dédoubler les multiplicateurs * quand on assemble avec l'assemblage de RESO car le traitement des * multiplicateurs dans ldmt1 l'impose (simple multiplicateur non * prévu) MRIGI0=MRIGID *dbg write(ioimp,*) 'LASRIG=',LASRIG,' NELIM=',NELIM IF (MRIGID.NE.0) THEN IF (IIMPI.NE.0) THEN WRITE(IOIMP,*) $ '*** ELIMINATION DES MULTIPLICATEURS DE LAGRANGE (LX) KRES6B' ENDIF SEGINI IDEMEM IDEMEM(**)=KSMBR SEGINI IDEME0,IDEME1 CALL KRES6B(MRIGID,IDEMEM,IDEME0,IDEME1,NELIM, $ MRIGIC,ICOND,NPASS) IF (IERR.NE.0) RETURN KSMBRC=IDEMEM(1) IF (LASRIG) THEN IF (IIMPI.NE.0) THEN WRITE(IOIMP,*) $ '*** ASSEMBLAGE RENUMEROTATION RESOU KRES8' ENDIF * Gestion de la normalisation NORICO=NORINC NORIDO=NORIND IF (ISCAL.EQ.0) THEN NORINC=0 NORIND=0 ELSE NORINC=-1 NORIND=0 ENDIF * Gestion de la renumérotation NUCROO=NUCROU IF (MRENU.EQ.'RIEN') THEN NUCROU=-1 * La renumérotation sera en fait Reverse Cuthill-McKee dans NUMOPT ELSEIF (MRENU.EQ.'RCM ') THEN NUCROU=1 ELSEIF (MRENU.EQ.'NEDI') THEN NUCROU=2 * La renumérotation sera en fait Nested Dissection dans NUMOPT ELSEIF (MRENU.EQ.'SLOA') THEN NUCROU=3 * ELSEIF (MRENU.EQ.'GIPR'.OR.MRENU.EQ.'GIBA') THEN * NUCROU=4 ELSE WRITE(IOIMP,*) 'MRENU=',MRENU CALL ERREUR(5) RETURN ENDIF CALL KRES8(MRIGIC,KSMBRC,INORMU, $ KTYPI,ITER,RESID,ICALRS,IRSTRT,LBCG,BRTOL,IDDOT,IMVEC, $ IORINC,MLAG1,MLAG2, $ KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV, $ KTIME,LTIME,LDUMP,ISMOOT, $ MCHSOC,LRES,LNMV,ICVG,IMPR) IF (IERR.NE.0) RETURN IDEMEM(1)=MCHSOC * Gestion de la normalisation NORINC=NORICO NORIND=NORIDO NUCROU=NUCROO IF (LTIME) CALL ECROBJ('TABLE ',KTIME) IF (MTINV.NE.0) THEN CALL ECME(MTINV,'CVGOK',ICVG) IF (LRES.NE.0) CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES) IF (LNMV.NE.0) CALL ECMO(MTINV,'NMATVEC','LISTENTI',LNMV) ENDIF * On décondense si nécessaire * * write (6,*) ' resou - mchsol ' * call ecchpo(mchsol,0) * call mucpri(mchsol,mrigid,iresi) * write (6,*) ' kres - iresi ' * call ecchpo(iresi,0) CALL KRES7B(MRIGIC,IDEMEM,IDEME0,IDEME1,NPASS) IF (IERR.NE.0) RETURN MCHSOL=IDEMEM(1) SEGSUP IDEMEM SEGSUP IDEME0,IDEME1 SEGACT MCHSOL*MOD MCHSOL.JATTRI(1)=1 CALL ACTOBJ('CHPOINT ',MCHSOL,1) CALL ECROBJ('CHPOINT ',MCHSOL) RETURN ELSE IF (IIMPI.NE.0) THEN WRITE(IOIMP,*) $ '*** TRANSFORMATION RIGIDITE -> MATRIK' WRITE(IOIMP,*) $ '*** ASSEMBLAGE RENUMEROTATION KRES2' ENDIF CALL ECROBJ('RIGIDITE',MRIGIC) CALL RIMA IF (IERR.NE.0) GOTO 9999 CALL MACHI2(1) IF (IERR.NE.0) GOTO 9999 CALL LIROBJ('MATRIK',MATRIK,1,IRET) IF(IRET.EQ.0) GOTO 9999 * Changement des noms d'inconnues du second membre IF (KSMBRC.NE.0) THEN CALL ECROBJ('CHPOINT ',KSMBRC) CALL MACHI2(1) CALL LIROBJ('CHPOINT ',KSMBR,1,IRET) IF (IRET.EQ.0) GOTO 9999 ENDIF ENDIF * write (6,*) ' le vecteur 2' * call ecchpo(ksmbr,0) * write (6,*) ' la matrice 2' * call ecrent(5) * call ecmatk(matrik) ENDIF * SEGACT MATRIK NMATRI=IRIGEL(/2) IF(NMATRI.EQ.0)THEN C% Résolution impossible : la matrice de RIGIDITE est vide CALL ERREUR(727) RETURN ENDIF SEGDES MATRIK IF (MATASS.EQ.0) MATASS=MATRIK IF (MAPREC.EQ.0) MAPREC=MATRIK * WRITE(IOIMP,*) 'Sortie de prkres' * WRITE(IOIMP,*) 'IOUBL=',IOUBL C IF (LTIME) THEN CALL CRTABL(KTIME) call timespv(ittime,oothrd) ITI1=(ITTIME(1)+ITTIME(2))/10 ELSE KTIME=0 ENDIF *STAT CALL PRMSTA('Lectures',MSTAT,IMPR) * C C Assemblage proprement dit C IIMPR=0 CALL KRES3(MATRIK,MATASS,MRENU,MMULAG,METASS, $ KTYPI,IORINC,MLAG1,MLAG2,IPBLOC, $ KTIME,LTIME, $ IIMPR,IRET) * Gestion du CTRL-C if (ierr.NE.0) return IF (IRET.NE.0) GOTO 9999 *! WRITE(IOIMP,*) 'Aprés assemblage' *! CALL ECRENT(5) *! CALL ECROBJ('MATRIK ',MATRIK) *! CALL PRLIST IF (LTIME) THEN call timespv(ittime,oothrd) ITI2=(ITTIME(1)+ITTIME(2))/10 ENDIF *STAT CALL PRMSTA('Assemblage',MSTAT,IMPR) * * "Oubli" des valeurs des matrice élémentaires * On met les tableaux de LIZAFM à 0 => à MENAGE de les supprimmer * si besoin est. * IOUBD=MOD(IOUBL,10) *! WRITE(IOIMP,*) 'IOUBD=',IOUBD IF (IOUBD.EQ.1) THEN CALL OUBIMA(MATRIK,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 IF (IMPR.GT.2) THEN WRITE(IOIMP,*) 'Oubli des mat. elem.' ENDIF ELSEIF (IOUBD.EQ.2) THEN call ooohor(0) SEGACT MATRIK*MOD LDETR=.FALSE. NMATE=IRIGEL(/2) DO IMATE=1,NMATE IMATRI=IRIGEL(4,IMATE) SEGACT IMATRI*MOD NSOUM =LIZAFM(/1) NTOTIN=LIZAFM(/2) DO ITOTIN=1,NTOTIN DO ISOUM=1,NSOUM IZAFM=LIZAFM(ISOUM,ITOTIN) IF (IZAFM.NE.0) THEN LDETR=.TRUE. SEGSUP IZAFM LIZAFM(ISOUM,ITOTIN)=0 ENDIF ENDDO ENDDO SEGDES IMATRI ENDDO IF (IMPR.GT.2.AND.LDETR) THEN WRITE(IOIMP,*) 'Destruction des mat. elem.' ENDIF ELSEIF (IOUBD.NE.0) THEN WRITE(IOIMP,*) 'IOUBL=',IOUBL, ' non prevu' GOTO 9999 ENDIF *STAT CALL PRMSTA('Oubli',MSTAT,IMPR) *! WRITE(IOIMP,*) 'Aprés oubli' C C Méthode directe C IF (KTYPI.EQ.1) THEN CALL KRES4(MATRIK,KCLIM,KSMBR, $ ISCAL, $ MCHSOL, $ IMPR,IRET) if (ierr.ne.0) return IF (IRET.NE.0) GOTO 9999 *STAT CALL PRMSTA('Methode directe',MSTAT,IMPR) C C Méthodes itératives C ELSE C CALL KRES5(MATRIK,KCLIM,KSMBR,KTYPI, $ MCHINI,ITER,RESID, $ BRTOL,IRSTRT,LBCG,ICALRS, $ MAPREC,KPREC, $ RXMILU,RXILUP,XLFIL,XDTOL,XSPIV, $ ISCAL, $ KTIME,LTIME,LDUMP,ISMOOT,IDDOT,IMVEC,IPBLOC, $ MCHSOL,LRES,LNMV,ICVG, $ IMPR,IRET) if (ierr.ne.0) return IF (IRET.NE.0) GOTO 9999 *STAT CALL PRMSTA('Methode iterative',MSTAT,IMPR) IF (MTINV.NE.0) THEN CALL ECME(MTINV,'CVGOK',ICVG) CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES) CALL ECMO(MTINV,'NMATVEC','LISTENTI',LNMV) ENDIF ENDIF IF (LTIME) THEN call timespv(ittime,oothrd) ITI3=(ITTIME(1)+ITTIME(2))/10 CHARI='KRES ASS+RENU' CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, $ 'ENTIER ',ITI2-ITI1,XVALR,CHARR,LOGIR,IRETR) CHARI='KRES PRE+RESO' CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, $ 'ENTIER ',ITI3-ITI2,XVALR,CHARR,LOGIR,IRETR) CHARI='KRES TOTAL ' CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, $ 'ENTIER ',ITI3-ITI1,XVALR,CHARR,LOGIR,IRETR) SEGDES KTIME CALL ECROBJ('TABLE ',KTIME) ENDIF IOUBE=IOUBL/10 *! WRITE(IOIMP,*) 'IOUBE=',IOUBE IF (IOUBE.GE.1) THEN call ooohor(0) SEGACT MATRIK*MOD IF (IOUBE.EQ.2) THEN PMORS=KIDMAT(4) IF (PMORS.NE.0) THEN IF (IMPR.GT.2) THEN WRITE(IOIMP,*) 'Destruction du profil morse' ENDIF SEGSUP PMORS KIDMAT(4)=0 ENDIF ENDIF IZA=KIDMAT(5) IF (IZA.NE.0) THEN IF (IMPR.GT.2) THEN WRITE(IOIMP,*) 'Destruction des valeurs morses' ENDIF SEGSUP IZA KIDMAT(5)=0 ENDIF PMORS=KIDMAT(6) IF (PMORS.NE.0) THEN IF (IMPR.GT.2) THEN WRITE(IOIMP,*) 'Destruction du profil du precon' ENDIF SEGSUP PMORS KIDMAT(6)=0 ENDIF IZA=KIDMAT(7) IF (IZA.NE.0) THEN IF (IMPR.GT.2) THEN WRITE(IOIMP,*) 'Destruction des valeurs du precon' ENDIF SEGSUP IZA KIDMAT(7)=0 ENDIF SEGDES MATRIK ELSEIF (IOUBE.NE.0) THEN WRITE(IOIMP,*) 'IOUBL=',IOUBL, ' non prevu' GOTO 9999 ENDIF * * On décondense si nécessaire * * write (6,*) ' resou - mchsol ' * call ecchpo(mchsol,0) * call mucpri(mchsol,mrigid,iresi) * write (6,*) ' kres - iresi ' * call ecchpo(iresi,0) * WRITE(IOIMP,*) 'Avant KRES7' IF (MRIGI0.NE.0) THEN IDEMEM(1)=MCHSOL CALL KRES7B(MRIGIC,IDEMEM,IDEME0,IDEME1,NPASS) IF (IERR.NE.0) RETURN MCHSOL=IDEMEM(1) SEGSUP IDEMEM SEGSUP IDEME0,IDEME1 ENDIF SEGACT MCHSOL*MOD MCHSOL.JATTRI(1)=1 CALL ACTOBJ('CHPOINT ',MCHSOL,1) CALL ECROBJ('CHPOINT ',MCHSOL) *STAT CALL SUMSTA(MSTAT,IMPR) * * Normal termination * RETURN * * Format handling * * * Error handling * 9999 CONTINUE WRITE(IOIMP,*) 'An error was detected in kres2.eso' * 153 2 * Opération illicite dans ce contexte CALL ERREUR(153) RETURN * * End of KRES2 * END