Sources Esope | Cast3M

mdilus
C MDILUS    SOURCE    CHAT      05/01/13    01:40:17     5004
      SUBROUTINE MDILUS(N,NNZ,
     $     ROWPTR,COLIND,VAL,
     $     DINV,
     $     IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : MDILUS
C DESCRIPTION :
C     Calcul du préconditionneur D-ILU d'une matrice Morse
C     SYMETRIQUE.
C     D-ILU : Diagonal Incomplete LU factorization
C
C     Ce sous-programme est appelé par prdilu.eso.
C     Il est écrit en Fortran car il est censé etre rapide
C
C LANGAGE     : FORTRAN 77 + chouhia ESOPE (E/S)
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
C               mél : gounand@semt2.smts.cea.fr
C REFERENCE (bibtex-like) :
C  @BOOK{templates,
C        AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
C                J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
C                H. Van der Vorst},
C        TITLE={Templates for the Solution of Linear Systems :
C               Building Blocks for Iterative Methods},
C        PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
C  -> URL : http://www.netlib.org/templates/Templates.html
C***********************************************************************
C APPELES : -
C***********************************************************************
C ENTREES            : N, NNZ,
C                      ROWPTR, COLIND, VAL
C ENTREES/SORTIES    : -
C SORTIES            : DINVPIV, IRET
C CODE RETOUR (IRET) : 0 si ok
C                     &lt;0 si problème (diagonale nulle...)
C     N      : nombre de degrés de liberté
C     NNZ    : nombre de valeurs non nulles de la matrice Morse
C     ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
C                           et valeurs de la matrice Morse
C     DINVPIV : vecteur contenant l'inverse de la factorisation
C               D-ILU de la matrice morse ci-dessus.
C***********************************************************************
C VERSION    : v1, 01/04/98, version initiale
C HISTORIQUE : v1, 01/04/98, création
C HISTORIQUE :
C HISTORIQUE :
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
*     .. Scalar Arguments ..
      INTEGER  N,NNZ
      INTEGER  IRET
*     ..
*     .. Array Arguments ..
*     .. Matrices stockées en Morse
      INTEGER ROWPTR(N+1)
      INTEGER COLIND(NNZ)
      REAL*8  VAL(NNZ)
*     .. Vecteurs
      REAL*8  DINV(N)
*
*     .. Variables locales
*     .. Parameters
      REAL*8     ZERO      ,PETIT        ,ONE
      PARAMETER (ZERO=0.0D0,PETIT=1.0D-14,ONE=1.0D0)
*     ..
C     Nombre de pivots petits
      INTEGER NBPIVP
C     Nombre de pivots inférieurs à 0
      INTEGER NBPIVI
      REAL*8 VALPIV,VALDIA
C***
      IRET=0
      NBPIVP=0
      NBPIVI=0
C     On boucle sur les lignes
      DO 1 IN=1,N
         IDEB=ROWPTR(IN)
         IFIN=ROWPTR(IN+1)-1
         ICOL=IDEB
         IF (IDEB.LE.IFIN) THEN
C     On cherche le terme Aii
 11         CONTINUE
            IF (COLIND(ICOL).LT.IN.AND.ICOL.LT.IFIN) THEN
               ICOL=ICOL+1
               GOTO 11
            ENDIF
C     On ne l'a pas trouvé
            IF (COLIND(ICOL).NE.IN) THEN
               WRITE(IOIMP,*) 'diag.',IN,'inexistante'
               WRITE(IOIMP,*) 'le préconditionnement D-ILU'
               WRITE(IOIMP,*) 'est impossible.'
               IRET=-1
               GOTO 9999
            ENDIF
C     On l'a trouvé
C     Dii &lt;- Dii+Aii
C     Dii &lt;- 1/Dii
            VALDIA=VAL(ICOL)
            VALPIV=DINV(IN)+VALDIA
            IF (VALPIV.EQ.ZERO) THEN
               WRITE(IOIMP,*) 'Pivot',IN,'nul'
               WRITE(IOIMP,*) 'le préconditionnement D-ILU'
               WRITE(IOIMP,*) 'est impossible.'
               IRET=-2
               GOTO 9999
            ELSE
               IF (ABS(VALPIV).LT.PETIT) THEN
                  NBPIVP=NBPIVP+1
               ENDIF
               IF (VALPIV.LT.ZERO) THEN
                  NBPIVI=NBPIVI+1
               ENDIF
            ENDIF
            VALPIV=ONE/VALPIV
            DINV(IN)=VALPIV
C     pour j>i : Djj &lt;- Djj-AjiDiiAij
            JCOL=ICOL+1
            JFIN=IFIN
            IF (JCOL.LE.JFIN) THEN
 12            CONTINUE
               JN=COLIND(JCOL)
               AIJ=VAL(JCOL)
               DINV(JN)=DINV(JN)-(VALPIV*AIJ*AIJ)
               JCOL=JCOL+1
               IF (JCOL.LE.JFIN) GOTO 12
            ENDIF
         ELSE
            WRITE(IOIMP,*) 'Ligne',IN,'vide'
            WRITE(IOIMP,*) 'le préconditionnement D-ILU'
            WRITE(IOIMP,*) 'est impossible.'
            IRET=-3
            GOTO 9999
         ENDIF
 1    CONTINUE
*
*     Warning(s)
*
      IF (NBPIVP.GT.0) THEN
         WRITE(IOIMP,*) 'WARNING !'
         WRITE(IOIMP,*) 'mdilus :',NBPIVP,' |pivot(s)|&lt;',PETIT
      ENDIF
      IF (NBPIVI.GT.0) THEN
         WRITE(IOIMP,*) 'WARNING !'
         WRITE(IOIMP,*) 'mdilus :',NBPIVI,' pivot(s) &lt; 0'
      ENDIF
*
*     Normal termination
*
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in mdilus.eso'
      RETURN
*
*     End of MDILUS
*
      END