Télécharger la plaquette Cast3M


fr en

Télécharger memilu.eso

Retour à la liste

Numérotation des lignes : 

memilu

  1. C MEMILU    SOURCE    GOUNAND   06/03/06    21:18:28     5319
  2.       SUBROUTINE MEMILU(NTT,NNZ,A,JA,IA,
  3.      $     ALU,JLU,JU,
  4.      $     IWORK,
  5.      $     RXMILU,
  6.      $     IMPR,IRET)
  7.       IMPLICIT INTEGER(I-N)
  8.       IMPLICIT REAL*8 (A-H,O-Z)
  9. C***********************************************************************
  10. C NOM         : MEMILU
  11. C DESCRIPTION :
  12. C     Calcul du préconditionneur MILU(0) d'une matrice Morse
  13. C     MILU(0) : Modified Incomplete LU factorization of level 0
  14. C               appelée aussi Choleski ou Crout incomplet modifié !
  15. C     Ici, c'est meme un MILU(0) relaxé (avec RXMILU)...
  16. C     La matrice Morse est au format CSR (Column Sparse Row)
  17. C     Le préconditionneur est au format MSR (Modified Sparse Row,
  18. C     stockage de l'inverse de la diagonale) (ALU, JLU).
  19. C     JU est un tableau supplémentaire pour le repérage des diagonales.
  20. C     (voir la doc. de Sparskit version 2)
  21. C
  22. C LANGAGE     : FORTRAN 77 (sauf E/S)
  23. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  24. C               mél : gounand@semt2.smts.cea.fr
  25. C REFERENCE (bibtex-like) :
  26. C  @BOOK{templates,
  27. C        AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  28. C                J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  29. C                H. Van der Vorst},
  30. C        TITLE={Templates for the Solution of Linear Systems :
  31. C               Building Blocks for Iterative Methods},
  32. C        PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  33. C  -> URL : http://www.netlib.org/templates/Templates.html
  34. C  Sparskit : a basic tool kit for sparse matrix computations
  35. C  Version 2 (Youcef Saad)
  36. C  -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
  37. C
  38. C***********************************************************************
  39. C APPELES          : -
  40. C APPELE PAR       : PRILU0
  41. C***********************************************************************
  42. C ENTREES            : NTT, NNZ, IA, JA, A, RXMILU
  43. C ENTREES/SORTIES    : IWORK (tableau de travail)
  44. C SORTIES            : JU, JLU, ALU
  45. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  46. C     RXMILU : type REAL*8.
  47. C              facteur de relaxation (normalement compris
  48. C              entre 0.D0 et 1.D0)
  49. C              0.D0 : on se ramène à Choleski incomplet
  50. C              1.D0 : on se ramène à MILU(0) non relaxé.
  51. C***********************************************************************
  52. C VERSION    : v1.1, 21/03/00
  53. C HISTORIQUE : v1.1, 22/03/00, gounand
  54. C Amélioration de la détection des pivots nulles.
  55. C HISTORIQUE : v1, 20/12/99, création
  56. C HISTORIQUE :
  57. C HISTORIQUE :
  58. C***********************************************************************
  59. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  60. C en cas de modification de ce sous-programme afin de faciliter
  61. C la maintenance !
  62. C***********************************************************************
  63.  
  64. -INC PPARAM
  65. -INC CCOPTIO
  66. -INC CCREEL
  67. *     .. Scalar Arguments ..
  68.       INTEGER  NTT,NNZ
  69.       REAL*8 RXMILU
  70.       INTEGER  IMPR,IRET
  71. *     ..
  72. *     .. Array Arguments ..
  73.       REAL*8  A(NNZ),ALU(NNZ+1)
  74.       INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
  75.       INTEGER IWORK(NNZ)
  76. *
  77. *     .. Variables locales
  78. *     .. Parameters
  79.       REAL*8     ZERO      ,ONE
  80.       PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  81. *     ..
  82. C     Nombre de pivots petits
  83.       INTEGER NBPIVP
  84. C     Nombre de pivots inférieurs à 0
  85.       INTEGER NBPIVI
  86.       INTEGER ITT,JNZ,KNZ,JCOL,JROW
  87.       INTEGER JSTRT,JMID,JSTOP,JU0,JW
  88.       INTEGER ISTRT,ISTOP
  89.       REAL*8 TL,VALPIV,TNORM,PETIT
  90.       REAL*8 FILVAL
  91. *
  92. * Executable statements
  93. *
  94.       PETIT=XZPREC**(0.75D0)
  95.       NBPIVP=0
  96.       NBPIVI=0
  97.       JU0=NTT+2
  98.       JLU(1)=JU0
  99. *
  100. * On boucle sur les ddl
  101. *
  102.       DO 1 ITT=1,NTT
  103.          JSTRT=JU0
  104. *
  105. * Création de la ligne ITT pour L et U
  106. *
  107.          ISTRT=IA(ITT)
  108.          ISTOP=IA(ITT+1)-1
  109.          TNORM=ZERO
  110.          DO 12 JNZ=ISTRT,ISTOP
  111. *
  112. * Copie de la ligne ITT du format CSR au format MSR
  113. *
  114.             JCOL=JA(JNZ)
  115.             TNORM=TNORM+ABS(A(JNZ))
  116.             IF (JCOL.EQ.ITT) THEN
  117.                ALU(ITT)   =A(JNZ)
  118.                IWORK(JCOL)=ITT
  119.                JU(ITT)    =JU0
  120.             ELSE
  121.                ALU(JU0)   =A(JNZ)
  122.                JLU(JU0)   =JA(JNZ)
  123.                IWORK(JCOL)=JU0
  124.                JU0 = JU0+1
  125.             ENDIF
  126.  12      CONTINUE
  127.          TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  128.          IF (TNORM.LE.XPETIT) THEN
  129.             WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  130.             GOTO 9999
  131.          ENDIF
  132.          JLU(ITT+1)=JU0
  133.          JSTOP=JU0-1
  134.          JMID=JU(ITT)-1
  135.          FILVAL=ZERO
  136.          DO 14 JNZ=JSTRT,JMID
  137.             JROW=JLU(JNZ)
  138.             TL=ALU(JNZ)*ALU(JROW)
  139.             ALU(JNZ)=TL
  140. *
  141. * On effectue la combinaison linéaire
  142. *
  143.             DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
  144.                JW=IWORK(JLU(KNZ))
  145.                IF (JW.NE.0) THEN
  146.                   ALU(JW)=ALU(JW)-TL*ALU(KNZ)
  147.                ELSE
  148.                   FILVAL=FILVAL+TL*ALU(KNZ)
  149.                ENDIF
  150.  142        CONTINUE
  151.  14      CONTINUE
  152. *
  153. *  On s'occupe de l'élément diagonal
  154. *
  155.          ALU(ITT)=ALU(ITT)-(FILVAL*RXMILU)
  156.          VALPIV=ALU(ITT)
  157.          IF (VALPIV.EQ.ZERO) THEN
  158.             WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
  159.      $           'le préconditionnement MILU(0) est impossible.'
  160.             IRET=-2
  161.             GOTO 9999
  162.          ELSE
  163.             IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  164.                NBPIVP=NBPIVP+1
  165.             ENDIF
  166.             IF (VALPIV.LT.ZERO) THEN
  167.                NBPIVI=NBPIVI+1
  168.             ENDIF
  169.          ENDIF
  170.          ALU(ITT)=ONE/ALU(ITT)
  171. *
  172. * On remet le tableau de travail à zéro
  173. *
  174.          IWORK(ITT)=0
  175.          DO 16 JNZ=JSTRT,JSTOP
  176.             IWORK(JLU(JNZ))=0
  177.  16      CONTINUE
  178.  1    CONTINUE
  179. *
  180. * Warning(s)
  181. *
  182.       IF (IMPR.GT.0) THEN
  183.          IF (NBPIVP.GT.0) THEN
  184.             WRITE(IOIMP,*) 'WARNING !',
  185.      $           ' memilu :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  186.          ENDIF
  187. *         IF (NBPIVI.GT.0) THEN
  188. *            WRITE(IOIMP,*) 'WARNING !',
  189. *     $           ' memilu :',NBPIVI,' pivot(s) < 0'
  190. *         ENDIF
  191.       ENDIF
  192. *
  193. * Normal termination
  194. *
  195.       IRET=0
  196.       RETURN
  197. *
  198. * Error handling
  199. *
  200.  9999 CONTINUE
  201.       WRITE(IOIMP,*) 'An error was detected in memilu.eso'
  202.       RETURN
  203. *
  204. *     End of MEMILU
  205. *
  206.       END
  207.  
  208.  
  209.  
  210.  
  211.  
  212.  
  213.  
  214.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales