Télécharger kresll.eso

Retour à la liste

kresll
C KRESLL    SOURCE    GOUNAND   25/03/11    21:15:04     12185          
      SUBROUTINE KRESLL()
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
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
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
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 SMCHPOI
      POINTEUR KCLIM.MCHPOI
      POINTEUR KSMBR.MCHPOI
      POINTEUR MCHINI.MCHPOI
      POINTEUR MCHSOL.MCHPOI
-INC SMTABLE
      POINTEUR MTINV.MTABLE
*  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
C
      CHARACTER*8 TYPE
*     Paramètre m du GMRES(m)
      INTEGER IRSTRT
      INTEGER ITER
      REAL*8 BRTOL,RESID,RXMILU
*
      REAL*8 XLFIL,XDTOL
      INTEGER KPREC
      INTEGER NMATRI
      INTEGER IP,KTYPI
      INTEGER IMPINV,IIMPR
      CHARACTER*4 MRENU,MMULAG
      INTEGER IMPR,IRET
C
C     Définition des options
C
C     'IMPR' + un entier : niveau d'impression
C     'TYPI' + une table contenant les informations
C              sur le type de résolution
C     'CLIM' + un chpoint : conditions aux limites de Dirichlet
C     'SMBR' + un chpoint : second membre (forces...)
C     'PREC' + une matrice dont on va utiliser le préconditionnement
C
      INTEGER NBM
      PARAMETER (NBM=4)
      CHARACTER*4 LMOT(NBM)
C     Tableau de correspondance : entier <-> mot
C     pour le type d'inversion
      INTEGER NBTYPI
C      INTEGER LNTYPI
      PARAMETER (NBTYPI=5)
C      PARAMETER (LNTYPI=18)
C      CHARACTER*(LNTYPI) LTYPI(NBTYPI)
      LOGICAL LBID,LVAL,LTIME,LMRENU,LDEPE,LLDEPE
C
C Initialisation des tableaux
C
      DATA LMOT/'IMPR','TYPI','CLIM','SMBR'/
C      DATA LTYPI/'Choleski          ',
C     $     'Conjugate Gradient',
C     $     'BiCGSTAB          ',
C     $     'BiCGSTAB(2)       ',
C     $     'GMRES(m)          '/
C
      TYPE='MATRIK'
      CALL LIROBJ(TYPE,MATRIK,1,IRET)
      IF(IRET.EQ.0) GOTO 9999
c
      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
c Valeurs par défaut reprise de prkres
      KSMBR=0
      MTINV=0
      IMPR=0
      KCLIM=0
      KTYPI=1
      MATASS=MATRIK
      MAPREC=MATRIK
      MRENU='SLOA'
      LMRENU=.FALSE.
      MMULAG='APR2'
      ISCAL=0
      IOUBL=0
      IMPINV=0
      MCHINI=0
      ITER=2000
      RESID=1.D-10
      BRTOL=1.D-40
      IRSTRT=50
      KPREC=3
      RXMILU=1.D0
      RXILUP=0.5D0
      XLFIL=2.D0
      XDTOL=-1.D0
      XSPIV=0.1D0
      LBCG=4
      ICALRS=0
      METASS=5
      LTIME=.FALSE.
      KTIME=0
      LDEPE=.TRUE.
      LLDEPE=.FALSE.
      IDDOT=0
      IMVEC=2
      IORINC=0
C
 1    CONTINUE
      CALL LIRMOT(LMOT,NBM,IP,0)
      IF(IP.EQ.0) GOTO 2
      GOTO (11,12,13,14),IP
C     Lecture du niveau d'impression
C     option 'IMPR'
 11   CONTINUE
c     write(6,*)' option IMPR'
      CALL LIRENT(IMPR,1,IRET)
      IF (IRET.EQ.0) GOTO 9999
      GO TO 1
C     Lecture de la table pour les méthodes d'inversion
C     option 'TYPI'
 12   CONTINUE
      CALL LIRTAB('METHINV',MTINV,1,IRET)
      IF(IRET.EQ.0)RETURN
       CALL LIRENT(KTYPI,1,IRET)
C      KTYPI : Méthode d'inversion du système (cf. LTYPI)
C                   1 : résolution directe (Choleski)
C                   2 : Gradient Conjugué
C                   3 : Bi-Gradient Conjugué Stabilisé (BiCGSTAB)
C                   4 : BiCGSTAB(2)
C                   5 : GMRES(m)
C     CALL ACME(MTINV,'TYPINV',KTYPI)
      IF (IERR.NE.0) GOTO 9999
C     write(6,*)' Ktypi=',ktypi
      IF (KTYPI.LT.1.OR.KTYPI.GT.NBTYPI) THEN
         WRITE(IOIMP,*) 'L''indice TYPINV ',
     $        'is out of range (1..',NBTYPI,') :',KTYPI
         GOTO 9999
      ENDIF
      GO TO 1
C     Lecture du chpoint de conditions aux limites
C     option CLIM
 13   CONTINUE
      CALL QUETYP(TYPE,1,IRET)
      IF(IRET.EQ.0) GOTO 9999
      IF(TYPE.EQ.'CHPOINT')THEN
         CALL LIROBJ(TYPE,KCLIM,1,IRET)
      ELSE
         WRITE(IOIMP,*) 'Erreur lecture C.lim'
         GOTO 9999
      ENDIF
      GO TO 1
C     Lecture du chpoint second membre
C     option SMBR
 14   CONTINUE
      CALL QUETYP(TYPE,1,IRET)
      IF(IRET.EQ.0)RETURN
      IF(TYPE.EQ.'CHPOINT')THEN
         CALL LIROBJ('CHPOINT',KSMBR,1,IRET)
c       write(6,*)' lecture 2 membre '
      ELSE
         WRITE(IOIMP,*) 'Erreur lecture second membre'
         GOTO 9999
      ENDIF
      GOTO 1
 2    CONTINUE
C     write(6,*)' On a pas finit les lectures '
c     NAT=2
c           NSOUPO=0
c           SEGINI MCHPO1
c           MCHPO1.JATTRI(1)=2
c     CALL ECROBJ('CHPOINT',MCHPO1)
c     return
C
C Matrice dont on utilise l'assemblage
C (éventuellement identique à la matrice transmise)
C
      TYPE='MATRIK  '
C     CALL ACMO(MTINV,'MATASS',TYPE,MATASS)
C     IF (IERR.NE.0) GOTO 9999
      CALL LIROBJ('MATRIK',MATASS,1,IRET)
      IF(IRET.EQ.0) GOTO 9999
C - Pour l'assemblage : type de renumérotation
C     CALL ACMM(MTINV,'TYRENU',MRENU)
C     IF (IERR.NE.0) GOTO 9999
      CALL LIRCHA(MRENU,1,LCHAR)
      IF(LCHAR.EQ.0)GO TO 9999
C - Pour l'assemblage : prise en compte des mult.lag
C     CALL ACMM(MTINV,'PCMLAG',MMULAG)
C     IF (IERR.NE.0) GOTO 9999
      CALL LIRCHA(MMULAG,1,LCHAR)
      IF(LCHAR.EQ.0)GO TO 9999
C - Pour l'assemblage : scaling de la matrice
C      CALL ACME(MTINV,'SCALING',ISCAL)
C      IF (IERR.NE.0) GOTO 9999
      ISCAL=0
C     Niveau d'impression pour la partie résolution itérative
c     CALL ACME(MTINV,'IMPINV',IMPINV)
c     IF (IERR.NE.0) GOTO 9999
c     IMPR=MAX(IMPR,IMPINV)
      IMPR=0
C
C Assemblage proprement dit
C
      IIMPR=0
*     SG 2016/02/09 : non à la ligne suivante : il faut que METASS soit
*     égale à 5 (voir remarque dans makpr2)
*      METASS=4
      CALL KRES3(MATRIK,MATASS,MRENU,MMULAG,METASS,
     $     KTIME,LTIME,
     $     IIMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
*!      WRITE(IOIMP,*) 'Aprés assemblage'
*!      CALL ECRENT(5)
*!      CALL ECROBJ('MATRIK  ',MATRIK)
*!      CALL PRLIST
*
* "Oubli" des valeurs des matrice élémentaires
* On met les tableaux de LIZAFM à 0 => à MENAGE de les supprimmer
* si besoin est.
*
      CALL OUBIMA(MATRIK,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C
C Méthode directe
C
      IF (KTYPI.EQ.1) THEN
         CALL KRES4(MATRIK,KCLIM,KSMBR,
     $        ISCAL,
     $        MCHSOL,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
C
C Méthodes itératives
C
      ELSEIF(KTYPI.GE.2.AND.KTYPI.LE.NBTYPI)THEN
C Paramètres pour la méthode itérative proprement dite
C        - Champoint d'initialisation de la méthode
C            (i.e. estimation de l'inconnue)
         TYPE='CHPOINT '
C        CALL ACMO(MTINV,'XINIT',TYPE,MCHINI)
C        IF (IERR.NE.0) GOTO 9999
         CALL LIROBJ('CHPOINT',MCHINI,1,IRET)
         IF(IRET.EQ.0) GOTO 9999
C        - Nombre maxi d'itérations à effectuer
C        CALL ACME(MTINV,'NITMAX',ITER)
C        IF (IERR.NE.0) GOTO 9999
         CALL LIRENT(ITER,1,IRET)
         IF(IRET.EQ.0) GOTO 9999
C        - Norme maxi (L2 normé par le second membre) du résidu
C        CALL ACMF(MTINV,'RESID',RESID)
C        IF (IERR.NE.0) GOTO 9999
         CALL LIRREE(RESID,1,IRET)
         IF(IRET.EQ.0) GOTO 9999
C        - 'Breakdown tolerance' pour les méthodes de type
C           BiCGSTAB. Si un certain produit scalaire de vecteurs
C           "direction" est inférieur à cette tolérance, la
C           méthode s'arrete.
         IF (KTYPI.EQ.3.OR.KTYPI.EQ.4) THEN
C           CALL ACMF(MTINV,'BCGSBTOL',BRTOL)
C           IF (IERR.NE.0) GOTO 9999
         CALL LIRREE(BRTOL,1,IRET)
         IF(IRET.EQ.0) GOTO 9999
         ENDIF
C        - Paramètre m de redémarrage pour GMRES(m)
C          i.e. nombre max. de vecteurs par rapport auxquels on
C          orthogonalise le résidu
         IF (KTYPI.EQ.5) THEN
C            CALL ACME(MTINV,'GMRESTRT',IRSTRT)
C            IF (IERR.NE.0) GOTO 9999
            CALL LIRENT(IRSTRT,1,IRET)
            IF(IRET.EQ.0) GOTO 9999
         ENDIF
C Paramètres pour les différents préconditionneurs
C        Matrice dont on utilise le préconditionneur
C        (éventuellement identique à la matrice transmise).
C        TYPE='MATRIK  '
C        CALL ACMO(MTINV,'MAPREC',TYPE,MAPREC)
C        IF (IERR.NE.0) GOTO 9999
         CALL LIROBJ('MATRIK',MAPREC,1,IRET)
         IF(IRET.EQ.0) GOTO 9999
C        - Type de préconditionnement : (cf. LPREC)
C                   0 : pas de préconditionnement
C                   1 : préconditionnement par la diagonale
C                   2 : préconditionnement D-ILU
C                   3 : préconditionnement ILU(0) (Choleski)
C                   4 : préconditionnement MILU(0) (Choleski modifié)
C                   5 : préconditionnement ILUT (dual truncation)
C                   6 : préconditionnement ILUT2 (une variante du
C                       précédent qui remplit mieux la mémoire et
C                       fonctionne mieux quelquefois)
C        CALL ACME(MTINV,'PRECOND',KPREC)
C        IF (IERR.NE.0) GOTO 9999
         CALL LIRENT(KPREC,1,IRET)
         IF(IRET.EQ.0) GOTO 9999
C        - Paramètre de relaxation pour le préconditionnement
C          MILU(0) compris entre 0.D0 et 1.D0
C          S'il est égal à 0, on se ramène à ILU(0)
C          S'il est égal à 1, MILU(0) est dit non relaxé
         IF (KPREC.EQ.4) THEN
C           CALL ACMF(MTINV,'MILURELX',RXMILU)
C           IF (IERR.NE.0) GOTO 9999
            CALL LIRREE(RXMILU,1,IRET)
            IF(IRET.EQ.0) GOTO 9999
         ENDIF
C        - Pour une méthode ILUT, on a les deux indices suivant :
C          * ILUTLFIL : encombrement maximal (approximatif) du
C                       préconditionneur, par rapport à la matrice.
C          * ILUTDTOL : "drop tolerance" pour le préconditionneur.
C                       i.e. en-dessous de cette valeur relative, les
C                       termes de la factorisation incomplète seront
C                       oubliés.
         IF (KPREC.GE.5.AND.KPREC.LE.9) THEN
C           CALL ACMF(MTINV,'ILUTLFIL',XLFIL)
C           IF (IERR.NE.0) GOTO 9999
            CALL LIRREE(XLFIL,1,IRET)
            IF(IRET.EQ.0) GOTO 9999
C           CALL ACMF(MTINV,'ILUTDTOL',XDTOL)
C           IF (IERR.NE.0) GOTO 9999
            CALL LIRREE(XDTOL,1,IRET)
            IF(IRET.EQ.0) GOTO 9999
         ENDIF
C        - Pour une méthode ILUTP, on a les deux indices suivant :
C        * ILUTPPIV (type REEL) : sensibilité du pivoting
C         IF (KPREC.GE.7.AND.KPREC.LE.9) THEN
C            CALL ACMF(MTINV,'ILUTPPIV',XSPIV)
C         ENDIF
         CALL KRES5(MATRIK,KCLIM,KSMBR,KTYPI,
     $        MCHINI,ITER,RESID,
     $        BRTOL,IRSTRT,LBCG,ICALRS,
     $        MAPREC,KPREC,
     $        RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,
     $        ISCAL,
     $        KTIME,LTIME,IDDOT,IMVEC,IORINC,
     $        MCHSOL,LRES,LNMV,ICVG,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES)
      ELSE
         WRITE(IOIMP,*) 'KTYPI=',KTYPI,'out of range (1..',NBTYPI,')'
         GOTO 9999
      ENDIF
      CALL ECROBJ('CHPOINT ',MCHSOL)
*
* Normal termination
*
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in kresll.eso'
*  153 2
* Opération illicite dans ce contexte
      CALL ERREUR(153)
      RETURN
*
* End of KRES2
*
      END
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales