Sources Esope | Cast3M

meilu0
C MEILU0    SOURCE    GOUNAND   05/02/16    21:16:06     5029
      SUBROUTINE MEILU0(NTT,NNZ,A,JA,IA,
     $     ALU,JLU,JU,
     $     IWORK,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : MEILU0
C DESCRIPTION :
C     Calcul du préconditionneur ILU(0) d'une matrice Morse.
C     ILU(0) : Incomplete LU factorization of level 0
C              appelée aussi Choleski ou Crout incomplet
C     La matrice Morse est au format CSR (Column Sparse Row)
C     Le préconditionneur est au format MSR (Modified Sparse Row,
C     stockage de l'inverse de la diagonale) (ALU, JLU).
C     JU est un tableau supplémentaire pour le repérage des diagonales.
C     (voir la doc. de Sparskit version 2)
C
C LANGAGE     : FORTRAN 77 (sauf E/S)
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
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  Sparskit : a basic tool kit for sparse matrix computations
C  Version 2 (Youcef Saad)
C  -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
C
C***********************************************************************
C APPELES          : -
C APPELE PAR       : PRILU0
C***********************************************************************
C ENTREES            : NTT, NNZ, IA, JA, A
C ENTREES/SORTIES    : IWORK (tableau de travail)
C SORTIES            : JU, JLU, ALU
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1.1, 21/03/00
C HISTORIQUE : v1.1, 22/03/00, gounand
C Amélioration de la détection des pivots nulles.
C HISTORIQUE : v1,   20/12/99, 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
-INC CCREEL
*     .. Scalar Arguments ..
      INTEGER  NTT,NNZ
      INTEGER  IMPR,IRET
*     ..
*     .. Array Arguments ..
      REAL*8  A(NNZ),ALU(NNZ+1)
      INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
      INTEGER IWORK(NTT)
*
*     .. Variables locales
*     .. Parameters
      REAL*8     ZERO      ,ONE
      PARAMETER (ZERO=0.0D0,ONE=1.0D0)
*     ..
C     Nombre de pivots petits
      INTEGER NBPIVP
C     Nombre de pivots inférieurs à 0
      INTEGER NBPIVI
      INTEGER ITT,JNZ,KNZ,JCOL,JROW
      INTEGER JSTRT,JMID,JSTOP,JU0,JW
      INTEGER ISTRT,ISTOP
      REAL*8 TL,VALPIV,TNORM,PETIT
 
*
* Executable statements
*
      PETIT=XZPREC**(0.75D0)
      NBPIVP=0
      NBPIVI=0
      JU0=NTT+2
      JLU(1)=JU0
*
* On boucle sur les ddl
*
      DO 1 ITT=1,NTT
         JSTRT=JU0
*
* Création de la ligne ITT pour L et U
*
         ISTRT=IA(ITT)
         ISTOP=IA(ITT+1)-1
         TNORM=ZERO
         DO 12 JNZ=ISTRT,ISTOP
*
* Copie de la ligne ITT du format CSR au format MSR
* et calcul de la norme de la ligne
*
            JCOL=JA(JNZ)
            TNORM=TNORM+ABS(A(JNZ))
            IF (JCOL.EQ.ITT) THEN
               ALU(ITT)   =A(JNZ)
               IWORK(JCOL)=ITT
               JU(ITT)    =JU0
            ELSE
               ALU(JU0)   =A(JNZ)
               JLU(JU0)   =JA(JNZ)
               IWORK(JCOL)=JU0
               JU0 = JU0+1
            ENDIF
 12      CONTINUE
         TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
         IF (TNORM.LE.XPETIT) THEN
            WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
            GOTO 9999
         ENDIF
         JLU(ITT+1)=JU0
         JSTOP=JU0-1
         JMID=JU(ITT)-1
         DO 14 JNZ=JSTRT,JMID
            JROW=JLU(JNZ)
            TL=ALU(JNZ)*ALU(JROW)
            ALU(JNZ)=TL
*
* On effectue la combinaison linéaire
*
            DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
               JW=IWORK(JLU(KNZ))
               IF (JW.NE.0) THEN
                  ALU(JW)=ALU(JW)-TL*ALU(KNZ)
               ENDIF
 142        CONTINUE
 14      CONTINUE
*
*  On s'occupe de l'élément diagonal
*
         VALPIV=ALU(ITT)
         IF (VALPIV.EQ.ZERO) THEN
            WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
     $           'le préconditionnement ILU(0) est impossible.'
            IRET=-2
            GOTO 9999
         ELSE
            IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
               NBPIVP=NBPIVP+1
            ENDIF
            IF (VALPIV.LT.ZERO) THEN
               NBPIVI=NBPIVI+1
            ENDIF
         ENDIF
         ALU(ITT)=ONE/ALU(ITT)
*
* On remet le tableau de travail à zéro
*
         IWORK(ITT)=0
         DO 16 JNZ=JSTRT,JSTOP
            IWORK(JLU(JNZ))=0
 16      CONTINUE
 1    CONTINUE
*
* Warning(s)
*
      IF (IMPR.GT.0) THEN
         IF (NBPIVP.GT.0) THEN
            WRITE(IOIMP,*) 'WARNING !',
     $           ' meilu0 :',NBPIVP,' |pivot(s)|/|ligne|&lt;',PETIT
         ENDIF
*!         IF (NBPIVI.GT.0) THEN
*!            WRITE(IOIMP,*) 'WARNING !',
*!     $           ' meilu0 :',NBPIVI,' pivot(s) &lt; 0'
*!         ENDIF
      ENDIF
*
* Normal termination
*
      IRET=0
      RETURN
*
* Error handling
*
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in meilu0.eso'
      RETURN
*
*     End of MEILU0
*
      END