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    PV        20/09/26    21:19:27     10724          
  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.       SEGDES MATRIK
  119. C     Le préconditionneur est-il déjà construit ?
  120.       IF ((ILUM.EQ.0).OR.(ILUI.EQ.0)) THEN
  121. C     La matrice Morse et son préconditionneur ont le meme
  122. C     profil
  123. *! Attention, on modifie la matrice dans meiltp
  124.          SEGACT AMORS*MOD
  125.          NTTDDL=AMORS.IA(/1)-1
  126.          NNZA=AMORS.JA(/1)
  127.          SEGACT AISA*MOD
  128. *
  129. * Ceci bugge avec ILUTPGO2 !!!!!
  130. *
  131. *
  132. *         XNTLIG=DBLE(NNZA)/DBLE(NTTDDL)
  133. *         ILFIL=INT(XNTLIG*XLFIL/2.D0)
  134. *         NNZMLU=NTTDDL*(2*ILFIL+1)
  135. *
  136.          NNZMLU = INT(DBLE(NNZA)*XLFIL) + 1
  137. *
  138.          NTT=NTTDDL-1
  139.          NJA=NNZMLU+1
  140.          SEGINI ILUM
  141.          NBVA=NNZMLU+1
  142.          SEGINI ILUI
  143.          JG=NTTDDL
  144.          SEGINI IWORK
  145.          JG=NTTDDL
  146.          SEGINI JWORK
  147.          JG=NTTDDL
  148.          SEGINI RWORK
  149.          JG=NTTDDL
  150.          SEGINI IPERM
  151.          JG=NTTDDL
  152.          SEGINI JPERM
  153. *         WRITE(IOIMP,*) 'NTTDDL=',NTTDDL,'NNZA=',NNZA,
  154. *     $        'NNZMLU=',NNZMLU
  155. *         WRITE(IOIMP,*) 'ILFIL=',ILFIL,'XLFIL=',XLFIL,'XDTOL=',XDTOL
  156. *         IF (XSPIV.LT.XZERO.OR.XSPIV.GT.1.D0) THEN
  157. *            WRITE(IOIMP,*) 'Error : XSPIV=',XSPIV
  158. *            GOTO 9999
  159. *         ENDIF
  160.          IF (IVARI.EQ.0) THEN
  161. C     Les boucles sont en Fortran pur
  162.             CALL MEILTP(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
  163.      $           XLFIL,XDTOL,XSPIV,NNZMLU,
  164.      $           ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
  165.      $           IWORK.LECT,JWORK.LECT,RWORK.PROG,
  166.      $           IPERM.LECT,JPERM.LECT,
  167.      $           IMPR,IRET)
  168.             IF (IRET.NE.0) GOTO 9999
  169.          ELSEIF (IVARI.EQ.1) THEN
  170.             JG=NTTDDL
  171.             SEGINI LWORK
  172.             CALL MILTP2(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
  173.      $           XLFIL,XDTOL,XSPIV,NNZMLU,
  174.      $           ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
  175.      $           IWORK.LECT,JWORK.LECT,RWORK.PROG,LWORK.LOGI,
  176.      $           IPERM.LECT,JPERM.LECT,
  177.      $           IMPR,IRET)
  178.             IF (IRET.NE.0) GOTO 9999
  179.             SEGSUP LWORK
  180.          ELSEIF (IVARI.EQ.2) THEN
  181.             JG=NTTDDL
  182.             SEGINI LWORK
  183.             JG=NTTDDL
  184.             SEGINI KWORK
  185.             JG=NTTDDL
  186.             SEGINI QWORK
  187.             CALL MILTP3(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
  188.      $           XLFIL,XDTOL,XSPIV,NNZMLU,
  189.      $           ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
  190.      $           IWORK.LECT,JWORK.LECT,RWORK.PROG,LWORK.LOGI,
  191.      $           KWORK.LECT,QWORK.PROG,
  192.      $           IPERM.LECT,JPERM.LECT,
  193.      $           IMPR,IRET)
  194.             IF (IRET.NE.0) GOTO 9999
  195.             SEGSUP QWORK
  196.             SEGSUP KWORK
  197.             SEGSUP LWORK
  198.          ELSE
  199.             WRITE(IOIMP,*) 'IVARI=',IVARI,' non prevu'
  200.          ENDIF
  201.          IF (IMPR.GT.1) THEN
  202.             WRITE(IOIMP,*) 'Préconditionneur : nb.termesstockés=',
  203.      $           NNZLU-1
  204.          ENDIF
  205. *
  206. *         SEGPRT,JPERM
  207.          SEGSUP IPERM
  208.          SEGSUP RWORK
  209.          SEGSUP JWORK
  210.          SEGSUP IWORK
  211.          NBVA=NNZLU
  212.          SEGADJ,ILUI
  213.          SEGDES ILUI
  214.          NTT=NTTDDL-1
  215.          NJA=NNZLU
  216.          SEGADJ,ILUM
  217.          SEGDES ILUM
  218.          SEGDES AISA
  219.          SEGDES AMORS
  220. *         SEGPRT,ILUM
  221. *         SEGPRT,ILUI
  222. *
  223. * On s'occupe des permutations
  224. *
  225.          SEGACT MATRIK*MOD
  226.          IF (.NOT.LRACOU) THEN
  227.             KMORS=MATRIK.KIDMAT(4)
  228.             SEGACT KMORS*MOD
  229.             DO IJA=KMORS.IA(1),KMORS.IA(NTTDDL+1)-1
  230.                KMORS.JA(IJA)=JPERM.LECT(KMORS.JA(IJA))
  231.             ENDDO
  232.             SEGDES KMORS
  233.          ENDIF
  234.          IDMATP=MATRIK.KIDMAT(1)
  235.          IDMATD=MATRIK.KIDMAT(2)
  236.          IF (IDMATD.EQ.0) THEN
  237.             WRITE(IOIMP,*) 'Erreur hyper grave'
  238.          ENDIF
  239.          IF (IDMATP.EQ.IDMATD) THEN
  240.             SEGINI,IDMATP=IDMATD
  241.          ELSE
  242.             SEGACT IDMATP*MOD
  243.          ENDIF
  244. * Copie de NUAN dans NUNA
  245.          DO ITTDDL=1,NTTDDL
  246.             IDMATP.NUNA(ITTDDL)=IDMATP.NUAN(ITTDDL)
  247.          ENDDO
  248. * Mise à jour de NUAN
  249.          DO ITTDDL=1,NTTDDL
  250.             IDMATP.NUAN(ITTDDL)=JPERM.LECT(IDMATP.NUNA(ITTDDL))
  251.          ENDDO
  252. * Mise à jour de NUNA
  253.          DO ITTDDL=1,NTTDDL
  254.             IDMATP.NUNA(IDMATP.NUAN(ITTDDL))=ITTDDL
  255.          ENDDO
  256.          INCX2=INCX
  257.          SEGACT INCX2
  258.          NBVA=NTTDDL
  259.          SEGINI,INCX
  260.          DO ITTDDL=1,NTTDDL
  261.             INCX.A(JPERM.LECT(ITTDDL))=INCX2.A(ITTDDL)
  262.          ENDDO
  263.          SEGDES INCX
  264.          SEGSUP INCX2
  265.          IF (NORMP.NE.0) THEN
  266.             NORMP2=NORMP
  267.             SEGACT NORMP2
  268.             JG=NTTDDL
  269.             SEGINI,NORMP
  270.             DO ITTDDL=1,NTTDDL
  271.                NORMP.PROG(JPERM.LECT(ITTDDL))=NORMP2.PROG(ITTDDL)
  272.             ENDDO
  273.             SEGDES NORMP
  274.             SEGSUP NORMP2
  275.          ENDIF
  276.          SEGSUP JPERM
  277. C         SEGDES IDMATP
  278.          SEGDES IDMATP
  279.          MATRIK.KIDMAT(1)=IDMATP
  280.          SEGDES MATRIK
  281. C
  282. C     On stocke la factorisation obtenue du préconditionneur
  283. C
  284.          SEGACT MATRIK*MOD
  285.          KIDMAT(6)=ILUM
  286.          KIDMAT(7)=ILUI
  287.          SEGDES MATRIK
  288. C
  289.          IF (IMPR.GT.6) THEN
  290.             WRITE(IOIMP,*) 'création du préconditionneur Morse',
  291.      $           ' de pointeurs',ILUM,'et',ILUI
  292.             IF (IMPR.GT.8) THEN
  293.                CALL ECMORS(ILUM,ILUI,(IMPR-1))
  294.             ENDIF
  295.          ENDIF
  296.       ELSE
  297.          IF (IMPR.GT.1) THEN
  298.             WRITE(IOIMP,*)
  299.      $           'Préconditionneur deja construit.'
  300.          ENDIF
  301.          IF (IMPR.GT.6) THEN
  302.             WRITE(IOIMP,*) 'Le préconditionneur est déjà construit :',
  303.      $           'ILUM=',ILUM,' et ILUI=',ILUI
  304.             IF (IMPR.GT.8) THEN
  305.                CALL ECMORS(ILUM,ILUI,(IMPR-1))
  306.             ENDIF
  307.          ENDIF
  308.       ENDIF
  309. *
  310. *     Normal termination
  311. *
  312.       RETURN
  313. *
  314. * Format handling
  315. *
  316. *
  317. * Error handling
  318. *
  319.  9999 CONTINUE
  320.       WRITE(IOIMP,*) 'An error was detected in priltp.eso'
  321.       IRET=1
  322.       RETURN
  323. *
  324. *     End of PRILTP
  325. *
  326.       END
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  

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