Sources Esope | Cast3M

prilut
C PRILUT    SOURCE    PV        20/09/26    21:19:30     10724          
      SUBROUTINE PRILUT(KMORS,KISA,MATRIK,XLFIL,XDTOL,IVARI,
     $     IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : PRILUT
C DESCRIPTION :
C     Calcul du préconditionneur ILUT d'une matrice Morse (IVARI=0)
C     ou d'une variante qui remplit mieux la mémoire et fonctionne mieux
C     quelquefois (IVARI=1)
C     Le préconditionneur est une matrice stockée
C     au format MSR (Modified Sparse Row, stockage de l'inverse de la
C     diagonale) de meme profil que la matrice Morse (format CSR) qu'il
C     préconditionne.
C     Le profil et les valeurs du préconditionneur sont
C     stockés dans KIDMAT(6 et 7) (réutilisation de l'existant).
C
C     Ce sous-programme est en fait une interface à :
C       meilut
C     qui est en Fortran presque pur (pour raison de rapidité)
C     et effectue la construction proprement dite du
C     préconditionneur.
C
C  ATTENTION : pour une matrice A quelconque, la factorisation
C  ---------   ILUT peut ne pas exister (pivot nul) ou avoir
C              des pivots négatifs MEME SI la factorisation
C              complète de A existe et n'a que des pivots
C              positifs.
C
C
C LANGAGE     : ESOPE
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 REMARQUES :
C On pourrait effectuer l'ordonnancement des colonnes du préconditionneur
C (peut-être que la montée-descente va plus vite après, il faudrait
C  tester)
C***********************************************************************
C APPELES       : MEILUT
C APPELE PAR    : KRES2
C***********************************************************************
C ENTREES            : MATRIK, XLFIL, XDTOL, IMPR
C SORTIES            : ILUM, ILUI (KIDMAT(6-7) dans MATRIK), IRET
C CODE RETOUR (IRET) : 0 si ok
C                     &lt;0 si problème
C     MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
C              on pioche dedans les informations nécessaires
C              (différents pointeurs, nb. de ddl...)
C     IMPR   : niveau d'impression
C     ILUM   : pointeur sur segment PMORS de l'include SMMATRIK
C              profil morse du préconditionneur ILU(0)
C              =KIDMAT(6)=KMORS=KIDMAT(4) dans MATRIK
C     ILUI   : pointeur sur segment IZA de l'include SMMATRIK
C              valeur du préconditionneur ILU(0)
C              =KIDMAT(7) dans MATRIK
C***********************************************************************
C VERSION    : v1, 23/02/2000, version initiale
C HISTORIQUE : v1, 23/02/2000, création
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 SMMATRIK
      INTEGER NTT,NJA,NBVA
      POINTEUR KMORS.PMORS
      POINTEUR KISA.IZA
      POINTEUR ILUM.PMORS
      POINTEUR ILUI.IZA
-INC SMLENTI
      INTEGER JG
      POINTEUR IWORK.MLENTI
      POINTEUR JWORK.MLENTI
-INC SMLREEL
      POINTEUR RWORK.MLREEL
*
      REAL*8 XLFIL,XDTOL
      INTEGER IVARI,IMPR,IRET
*
      INTEGER ILFIL,NNZA,NNZLU,NNZMLU,NTTDDL
      REAL*8 XNTLIG
*
      IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans prilut.eso'
C     On récupère les segments utiles
      SEGACT MATRIK
      ILUM =KIDMAT(6)
      ILUI =KIDMAT(7)
      SEGDES MATRIK
C     Le préconditionneur est-il déjà construit ?
      IF ((ILUM.EQ.0).OR.(ILUI.EQ.0)) THEN
C     La matrice Morse et son préconditionneur ont le meme
C     profil
         SEGACT KMORS
         NTTDDL=KMORS.IA(/1)-1
         NNZA=KMORS.JA(/1)
         SEGACT KISA
* Pour faire comme dans priltp
*         XNTLIG=DBLE(NNZA)/DBLE(NTTDDL)
*         ILFIL=INT(XNTLIG*XLFIL/2.D0)
*         NNZMLU=NTTDDL*(2*ILFIL+1)
*
         NNZMLU = INT(DBLE(NNZA)*XLFIL) + 1
*
         NTT=NTTDDL-1
         NJA=NNZMLU+1
         SEGINI ILUM
         NBVA=NNZMLU+1
         SEGINI ILUI
         JG=NTTDDL
         SEGINI IWORK
         JG=NTTDDL
         SEGINI JWORK
         JG=NTTDDL
         SEGINI RWORK
*         WRITE(IOIMP,*) 'NTTDDL=',NTTDDL,'NNZA=',NNZA,
*     $        'NNZMLU=',NNZMLU
*         WRITE(IOIMP,*) 'ILFIL=',ILFIL,'XLFIL=',XLFIL,'XDTOL=',XDTOL
C     Les boucles sont en Fortran pur
         CALL MEILUT(NTTDDL,NNZA,KISA.A,KMORS.JA,KMORS.IA,
     $        XLFIL,XDTOL,NNZMLU,
     $        ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
     $        IWORK.LECT,JWORK.LECT,RWORK.PROG,
     $        IVARI,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         IF (IMPR.GT.1) THEN
            WRITE(IOIMP,*) 'Préconditionneur : nb.termesstockés=',
     $           NNZLU-1
         ENDIF
*         SEGPRT,ILUM
*         SEGPRT,ILUI
         SEGSUP RWORK
         SEGSUP JWORK
         SEGSUP IWORK
         NBVA=NNZLU
         SEGADJ,ILUI
         SEGDES ILUI
         NTT=NTTDDL-1
         NJA=NNZLU
         SEGADJ,ILUM
         SEGDES ILUM
         SEGDES KISA
         SEGDES KMORS
*         SEGPRT,ILUM
*         SEGPRT,ILUI
* Ordonnancement (peut-être que la montée-descente va plus vite, il
*                 faudrait tester)
*
*         SEGACT ILUM*MOD
*         SEGACT ILUI*MOD
*         JG=MAX(NTTDDL+1,2*(NNZLU-(NTTDDL+1)))
*         SEGINI IWORK
*         CALL CSORT2(NTTDDL,ILUI.A(NTTDDL+2),ILUM.JA(NTTDDL+2),
*     $        ILUM.JA(1),IWORK.LECT,.TRUE.)
*         SEGSUP IWORK
*         SEGDES ILUI
*         SEGDES ILUM
C
C     On stocke la factorisation obtenue du préconditionneur
C
         SEGACT MATRIK*MOD
         KIDMAT(6)=ILUM
         KIDMAT(7)=ILUI
         SEGDES MATRIK
C
         IF (IMPR.GT.6) THEN
            WRITE(IOIMP,*) 'création du préconditionneur Morse',
     $           ' de pointeurs',ILUM,'et',ILUI
            IF (IMPR.GT.8) THEN
               CALL ECMORS(ILUM,ILUI,(IMPR-1))
            ENDIF
         ENDIF
      ELSE
         IF (IMPR.GT.6) THEN
            WRITE(IOIMP,*) 'Le préconditionneur est déjà construit :',
     $           'ILUM=',ILUM,' et ILUI=',ILUI
            IF (IMPR.GT.8) THEN
               CALL ECMORS(ILUM,ILUI,(IMPR-1))
            ENDIF
         ENDIF
      ENDIF
*
*     Normal termination
*
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in prilut.eso'
      RETURN
*
*     End of PRILUT
*
      END