Télécharger la plaquette Cast3M


fr en

Télécharger priltp.eso

Retour à la liste

Numérotation des lignes : 

priltp

  1. C PRILTP    SOURCE    GOUNAND   25/04/30    21:15:32     12258          
  2.       SUBROUTINE PRILTP(AMORS,AISA,MATRIK,XLFIL,XDTOL,XSPIV,
  3.      $     IVARI,
  4.      $     INCX,NORMP,LRACOU,
  5.      $     IMPR,IRET)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7.       IMPLICIT INTEGER (I-N)
  8. C***********************************************************************
  9. C NOM         : PRILTP
  10. C DESCRIPTION :
  11. C     Calcul du préconditionneur ILUTP d'une matrice Morse
  12. C     Le préconditionneur est une matrice stockée
  13. C     au format MSR (Modified Sparse Row, stockage de l'inverse de la
  14. C     diagonale) de meme profil que la matrice Morse (format CSR) qu'il
  15. C     préconditionne.
  16. C     Le profil et les valeurs du préconditionneur sont
  17. C     stockés dans KIDMAT(6 et 7) (réutilisation de l'existant).
  18. C
  19. C     Ce sous-programme est en fait une interface à :
  20. C       meiltp
  21. C     qui est en Fortran presque pur (pour raison de rapidité)
  22. C     et effectue la construction proprement dite du
  23. C     préconditionneur.
  24. C
  25. C  ATTENTION : pour une matrice A quelconque, la factorisation
  26. C  ---------   ILUTP peut ne pas exister (pivot nul) ou avoir
  27. C              des pivots négatifs MEME SI la factorisation
  28. C              complète de A existe et n'a que des pivots
  29. C              positifs.
  30. C
  31. C
  32. C LANGAGE     : ESOPE
  33. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  34. C               mél : gounand@semt2.smts.cea.fr
  35. C REFERENCE (bibtex-like) :
  36. C  @BOOK{templates,
  37. C        AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  38. C                J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  39. C                H. Van der Vorst},
  40. C        TITLE={Templates for the Solution of Linear Systems :
  41. C               Building Blocks for Iterative Methods},
  42. C        PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  43. C  -> URL : http://www.netlib.org/templates/Templates.html
  44. C  Sparskit : a basic tool kit for sparse matrix computations
  45. C  Version 2 (Youcef Saad)
  46. C  -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
  47. C REMARQUES :
  48. C On pourrait effectuer l'ordonnancement des colonnes du préconditionneur
  49. C (peut-être que la montée-descente va plus vite après, il faudrait
  50. C  tester)
  51. C***********************************************************************
  52. C APPELES       : MEILUT
  53. C APPELE PAR    : KRES2
  54. C***********************************************************************
  55. C ENTREES            : MATRIK, XLFIL, XDTOL, IMPR
  56. C SORTIES            : ILUM, ILUI (KIDMAT(6-7) dans MATRIK), IRET
  57. C CODE RETOUR (IRET) : 0 si ok
  58. C                     <0 si problème
  59. C     MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
  60. C              on pioche dedans les informations nécessaires
  61. C              (différents pointeurs, nb. de ddl...)
  62. C     IMPR   : niveau d'impression
  63. C***********************************************************************
  64. C VERSION    : v1, 08/04/2004, version initiale
  65. C HISTORIQUE : v1, 08/04/2004, création
  66. C HISTORIQUE :
  67. C***********************************************************************
  68. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  69. C en cas de modification de ce sous-programme afin de faciliter
  70. C la maintenance !
  71. C***********************************************************************
  72.  
  73. -INC PPARAM
  74. -INC CCOPTIO
  75. -INC CCREEL
  77.       INTEGER NTT,NJA,NBVA
  78.       POINTEUR AMORS.PMORS
  79.       POINTEUR AISA.IZA
  80.       POINTEUR KMORS.PMORS
  81.       POINTEUR KISA.IZA
  82.       POINTEUR ILUM.PMORS
  83.       POINTEUR ILUI.IZA
  84.       POINTEUR IDMATP.IDMAT
  85.       POINTEUR IDMATD.IDMAT
  86.       POINTEUR INCX.IZA
  87.       POINTEUR INCX2.IZA
  88. -INC SMLENTI
  89.       INTEGER JG
  90.       POINTEUR IWORK.MLENTI
  91.       POINTEUR JWORK.MLENTI
  92.       POINTEUR KWORK.MLENTI
  93.       POINTEUR IPERM.MLENTI
  94.       POINTEUR JPERM.MLENTI
  95. -INC SMLREEL
  96.       POINTEUR RWORK.MLREEL
  97.       POINTEUR QWORK.MLREEL
  98.       POINTEUR NORMP.MLREEL
  99.       POINTEUR NORMP2.MLREEL
  100. * -INC SMLLOGI
  101.       SEGMENT MLLOGI
  102.       LOGICAL LOGI(JG)
  103.       ENDSEGMENT
  104.       POINTEUR LWORK.MLLOGI
  105. *
  106.       REAL*8 XLFIL,XDTOL
  107.       INTEGER IMPR,IRET
  108. *
  109.       INTEGER ILFIL,NNZA,NNZLU,NNZMLU,NTTDDL
  110.       REAL*8 XNTLIG
  111.       LOGICAL LRACOU
  112. *
  113.       IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans priltp.eso'
  114. C     On récupère les segments utiles
  115.       SEGACT MATRIK
  116.       ILUM =KIDMAT(6)
  117.       ILUI =KIDMAT(7)
  118. C     Le préconditionneur est-il déjà construit ?
  119.       IF ((ILUM.EQ.0).OR.(ILUI.EQ.0)) THEN
  120. C     La matrice Morse et son préconditionneur ont le meme
  121. C     profil
  122. *! Attention, on modifie la matrice dans meiltp
  123.          SEGACT AMORS*MOD
  124.          NTTDDL=AMORS.IA(/1)-1
  125.          NNZA=AMORS.JA(/1)
  126.          SEGACT AISA*MOD
  127. *
  128. * Ceci bugge avec ILUTPGO2 !!!!!
  129. *
  130. *
  131. *         XNTLIG=DBLE(NNZA)/DBLE(NTTDDL)
  132. *         ILFIL=INT(XNTLIG*XLFIL/2.D0)
  133. *         NNZMLU=NTTDDL*(2*ILFIL+1)
  134. *
  135.          NNZMLU = INT(DBLE(NNZA)*XLFIL) + 1
  136. *
  137.          NTT=NTTDDL-1
  138.          NJA=NNZMLU+1
  139.          SEGINI ILUM
  140.          NBVA=NNZMLU+1
  141.          SEGINI ILUI
  142.          JG=NTTDDL
  143.          SEGINI IWORK
  144.          JG=NTTDDL
  145.          SEGINI JWORK
  146.          JG=NTTDDL
  147.          SEGINI RWORK
  148.          JG=NTTDDL
  149.          SEGINI IPERM
  150.          JG=NTTDDL
  151.          SEGINI JPERM
  152. *         WRITE(IOIMP,*) 'NTTDDL=',NTTDDL,'NNZA=',NNZA,
  153. *     $        'NNZMLU=',NNZMLU
  154. *         WRITE(IOIMP,*) 'ILFIL=',ILFIL,'XLFIL=',XLFIL,'XDTOL=',XDTOL
  155. *         IF (XSPIV.LT.XZERO.OR.XSPIV.GT.1.D0) THEN
  156. *            WRITE(IOIMP,*) 'Error : XSPIV=',XSPIV
  157. *            GOTO 9999
  158. *         ENDIF
  159.          IF (IVARI.EQ.0) THEN
  160. C     Les boucles sont en Fortran pur
  161.             CALL MEILTP(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
  162.      $           XLFIL,XDTOL,XSPIV,NNZMLU,
  163.      $           ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
  164.      $           IWORK.LECT,JWORK.LECT,RWORK.PROG,
  165.      $           IPERM.LECT,JPERM.LECT,
  166.      $           IMPR,IRET)
  167.             IF (IRET.NE.0) GOTO 9999
  168.          ELSEIF (IVARI.EQ.1) THEN
  169.             JG=NTTDDL
  170.             SEGINI LWORK
  171.             CALL MILTP2(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
  172.      $           XLFIL,XDTOL,XSPIV,NNZMLU,
  173.      $           ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
  174.      $           IWORK.LECT,JWORK.LECT,RWORK.PROG,LWORK.LOGI,
  175.      $           IPERM.LECT,JPERM.LECT,
  176.      $           IMPR,IRET)
  177.             IF (IRET.NE.0) GOTO 9999
  178.             SEGSUP LWORK
  179.          ELSEIF (IVARI.EQ.2) THEN
  180.             JG=NTTDDL
  181.             SEGINI LWORK
  182.             JG=NTTDDL
  183.             SEGINI KWORK
  184.             JG=NTTDDL
  185.             SEGINI QWORK
  186.             CALL MILTP3(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
  187.      $           XLFIL,XDTOL,XSPIV,NNZMLU,
  188.      $           ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
  189.      $           IWORK.LECT,JWORK.LECT,RWORK.PROG,LWORK.LOGI,
  190.      $           KWORK.LECT,QWORK.PROG,
  191.      $           IPERM.LECT,JPERM.LECT,
  192.      $           IMPR,IRET)
  193.             IF (IRET.NE.0) GOTO 9999
  194.             SEGSUP QWORK
  195.             SEGSUP KWORK
  196.             SEGSUP LWORK
  197.          ELSE
  198.             WRITE(IOIMP,*) 'IVARI=',IVARI,' non prevu'
  199.          ENDIF
  200.          IF (IMPR.GT.1) THEN
  201.             WRITE(IOIMP,*) 'Préconditionneur : nb.termesstockés=',
  202.      $           NNZLU-1
  203.          ENDIF
  204. *
  205. *         SEGPRT,JPERM
  206.          SEGSUP IPERM
  207.          SEGSUP RWORK
  208.          SEGSUP JWORK
  209.          SEGSUP IWORK
  210.          NBVA=NNZLU
  211.          SEGADJ,ILUI
  212.          NTT=NTTDDL-1
  213.          NJA=NNZLU
  214.          SEGADJ,ILUM
  215. *         SEGPRT,ILUM
  216. *         SEGPRT,ILUI
  217. *
  218. * On s'occupe des permutations
  219. *
  220.          SEGACT MATRIK*MOD
  221.          IF (.NOT.LRACOU) THEN
  222.             KMORS=MATRIK.KIDMAT(4)
  223.             SEGACT KMORS*MOD
  224.             DO IJA=KMORS.IA(1),KMORS.IA(NTTDDL+1)-1
  225.                KMORS.JA(IJA)=JPERM.LECT(KMORS.JA(IJA))
  226.             ENDDO
  227.          ENDIF
  228.          IDMATP=MATRIK.KIDMAT(1)
  229.          IDMATD=MATRIK.KIDMAT(2)
  230.          IF (IDMATD.EQ.0) THEN
  231.             WRITE(IOIMP,*) 'Erreur hyper grave'
  232.          ENDIF
  233.          IF (IDMATP.EQ.IDMATD) THEN
  234.             SEGINI,IDMATP=IDMATD
  235.          ELSE
  236.             SEGACT IDMATP*MOD
  237.          ENDIF
  238. * Copie de NUAN dans NUNA
  239.          DO ITTDDL=1,NTTDDL
  240.             IDMATP.NUNA(ITTDDL)=IDMATP.NUAN(ITTDDL)
  241.          ENDDO
  242. * Mise à jour de NUAN
  243.          DO ITTDDL=1,NTTDDL
  244.             IDMATP.NUAN(ITTDDL)=JPERM.LECT(IDMATP.NUNA(ITTDDL))
  245.          ENDDO
  246. * Mise à jour de NUNA
  247.          DO ITTDDL=1,NTTDDL
  248.             IDMATP.NUNA(IDMATP.NUAN(ITTDDL))=ITTDDL
  249.          ENDDO
  250.          INCX2=INCX
  251.          SEGACT INCX2
  252.          NBVA=NTTDDL
  253.          SEGINI,INCX
  254.          DO ITTDDL=1,NTTDDL
  255.             INCX.A(JPERM.LECT(ITTDDL))=INCX2.A(ITTDDL)
  256.          ENDDO
  257.          SEGSUP INCX2
  258.          IF (NORMP.NE.0) THEN
  259.             NORMP2=NORMP
  260.             SEGACT NORMP2
  261.             JG=NTTDDL
  262.             SEGINI,NORMP
  263.             DO ITTDDL=1,NTTDDL
  264.                NORMP.PROG(JPERM.LECT(ITTDDL))=NORMP2.PROG(ITTDDL)
  265.             ENDDO
  266.             SEGSUP NORMP2
  267.          ENDIF
  268.          SEGSUP JPERM
  269.          MATRIK.KIDMAT(1)=IDMATP
  270. C
  271. C     On stocke la factorisation obtenue du préconditionneur
  272. C
  273.          SEGACT MATRIK*MOD
  274.          KIDMAT(6)=ILUM
  275.          KIDMAT(7)=ILUI
  276. C
  277.          IF (IMPR.GT.6) THEN
  278.             WRITE(IOIMP,*) 'création du préconditionneur Morse',
  279.      $           ' de pointeurs',ILUM,'et',ILUI
  280.             IF (IMPR.GT.8) THEN
  281.                CALL ECMORS(ILUM,ILUI,(IMPR-1))
  282.             ENDIF
  283.          ENDIF
  284.       ELSE
  285.          IF (IMPR.GT.1) THEN
  286.             WRITE(IOIMP,*)
  287.      $           'Préconditionneur deja construit.'
  288.          ENDIF
  289.          IF (IMPR.GT.6) THEN
  290.             WRITE(IOIMP,*) 'Le préconditionneur est déjà construit :',
  291.      $           'ILUM=',ILUM,' et ILUI=',ILUI
  292.             IF (IMPR.GT.8) THEN
  293.                CALL ECMORS(ILUM,ILUI,(IMPR-1))
  294.             ENDIF
  295.          ENDIF
  296.       ENDIF
  297. *
  298. *     Normal termination
  299. *
  300.       RETURN
  301. *
  302. * Format handling
  303. *
  304. *
  305. * Error handling
  306. *
  307.  9999 CONTINUE
  308.       WRITE(IOIMP,*) 'An error was detected in priltp.eso'
  309.       IRET=1
  310.       RETURN
  311. *
  312. *     End of PRILTP
  313. *
  314.       END
  315.  
  316.  

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