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 <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