Sources Esope | Cast3M

mdilun
C MDILUN    SOURCE    GOUNAND   06/03/06    21:18:10     5319
      SUBROUTINE MDILUN(N,NNZ,
     $     ROWPTR,COLIND,VAL,
     $     DINV,
     $     IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : MDILUN
C DESCRIPTION :
C     Calcul du préconditionneur D-ILU d'une matrice Morse
C     NON symétrique.
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 + chouilla 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
C      INTEGER TMPCOL(N)
C      REAL*8  TMPVAL(N)
      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
            ELSE
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
               JCOL=ICOL+1
               JFIN=IFIN
               IF (JCOL.LE.JFIN) THEN
C     On parcourt les valeurs Aij (j>i).
C     Pour chaque valeur, on regarde ligne j s'il y a une
C     valeur Aji correspondante.
C     Si oui, Djj &lt;- Djj - AjiDiiAij
C     Sinon, rien !
 12               CONTINUE
                  JN=COLIND(JCOL)
C     Recherche de Aji
                  KDEB=ROWPTR(JN)
                  KFIN=ROWPTR(JN+1)-1
                  KCOL=KDEB
 121              CONTINUE
                  IF (COLIND(KCOL).LT.IN.AND.KCOL.LT.KFIN) THEN
                     KCOL=KCOL+1
                     GOTO 121
                  ENDIF
C     On a trouve Aji !
                  IF (COLIND(KCOL).EQ.IN) THEN
                     DINV(JN)=DINV(JN)-
     $                    (VAL(KCOL)*VALPIV*VAL(JCOL))
                  ENDIF
C     On boucle...
                  JCOL=JCOL+1
                  IF (JCOL.LE.JFIN) GOTO 12
               ENDIF
            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,*) 'mdilun :',NBPIVP,' |pivot(s)|&lt;',PETIT
      ENDIF
*!      IF (NBPIVI.GT.0) THEN
*!         WRITE(IOIMP,*) 'WARNING !'
*!         WRITE(IOIMP,*) 'mdilun :',NBPIVI,' pivot(s) &lt; 0'
*!      ENDIF
*
*     Normal termination
*
      RETURN
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in mdilun.eso'
      RETURN
*
*     End of MDILUN
*
      END