Télécharger priltp.eso

Retour à la liste

Numérotation des lignes :

  1. C PRILTP SOURCE PV 16/11/17 22:01:10 9180
  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
  76. INTEGER NTT,NJA,NBVA
  77. POINTEUR AMORS.PMORS
  78. POINTEUR AISA.IZA
  79. POINTEUR KMORS.PMORS
  80. POINTEUR KISA.IZA
  81. POINTEUR ILUM.PMORS
  82. POINTEUR ILUI.IZA
  83. POINTEUR IDMATP.IDMAT
  84. POINTEUR IDMATD.IDMAT
  85. POINTEUR INCX.IZA
  86. POINTEUR INCX2.IZA
  87. -INC SMLENTI
  88. INTEGER JG
  89. POINTEUR IWORK.MLENTI
  90. POINTEUR JWORK.MLENTI
  91. POINTEUR KWORK.MLENTI
  92. POINTEUR IPERM.MLENTI
  93. POINTEUR JPERM.MLENTI
  94. -INC SMLREEL
  95. POINTEUR RWORK.MLREEL
  96. POINTEUR QWORK.MLREEL
  97. POINTEUR NORMP.MLREEL
  98. POINTEUR NORMP2.MLREEL
  99. * -INC SMLLOGI
  100. SEGMENT MLLOGI
  101. LOGICAL LOGI(JG)
  102. ENDSEGMENT
  103. POINTEUR LWORK.MLLOGI
  104. *
  105. REAL*8 XLFIL,XDTOL
  106. INTEGER IMPR,IRET
  107. *
  108. INTEGER ILFIL,NNZA,NNZLU,NNZMLU,NTTDDL
  109. REAL*8 XNTLIG
  110. LOGICAL LRACOU
  111. *
  112. IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans priltp.eso'
  113. C On récupère les segments utiles
  114. SEGACT MATRIK
  115. ILUM =KIDMAT(6)
  116. ILUI =KIDMAT(7)
  117. SEGDES MATRIK
  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. SEGDES ILUI
  213. NTT=NTTDDL-1
  214. NJA=NNZLU
  215. SEGADJ,ILUM
  216. SEGDES ILUM
  217. SEGDES AISA
  218. SEGDES AMORS
  219. * SEGPRT,ILUM
  220. * SEGPRT,ILUI
  221. *
  222. * On s'occupe des permutations
  223. *
  224. SEGACT MATRIK*MOD
  225. IF (.NOT.LRACOU) THEN
  226. KMORS=MATRIK.KIDMAT(4)
  227. SEGACT KMORS*MOD
  228. DO IJA=KMORS.IA(1),KMORS.IA(NTTDDL+1)-1
  229. KMORS.JA(IJA)=JPERM.LECT(KMORS.JA(IJA))
  230. ENDDO
  231. SEGDES KMORS
  232. ENDIF
  233. IDMATP=MATRIK.KIDMAT(1)
  234. IDMATD=MATRIK.KIDMAT(2)
  235. IF (IDMATD.EQ.0) THEN
  236. WRITE(IOIMP,*) 'Erreur hyper grave'
  237. ENDIF
  238. IF (IDMATP.EQ.IDMATD) THEN
  239. SEGINI,IDMATP=IDMATD
  240. ELSE
  241. SEGACT IDMATP*MOD
  242. ENDIF
  243. * Copie de NUAN dans NUNA
  244. DO ITTDDL=1,NTTDDL
  245. IDMATP.NUNA(ITTDDL)=IDMATP.NUAN(ITTDDL)
  246. ENDDO
  247. * Mise à jour de NUAN
  248. DO ITTDDL=1,NTTDDL
  249. IDMATP.NUAN(ITTDDL)=JPERM.LECT(IDMATP.NUNA(ITTDDL))
  250. ENDDO
  251. * Mise à jour de NUNA
  252. DO ITTDDL=1,NTTDDL
  253. IDMATP.NUNA(IDMATP.NUAN(ITTDDL))=ITTDDL
  254. ENDDO
  255. INCX2=INCX
  256. SEGACT INCX2
  257. NBVA=NTTDDL
  258. SEGINI,INCX
  259. DO ITTDDL=1,NTTDDL
  260. INCX.A(JPERM.LECT(ITTDDL))=INCX2.A(ITTDDL)
  261. ENDDO
  262. SEGDES INCX
  263. SEGSUP INCX2
  264. IF (NORMP.NE.0) THEN
  265. NORMP2=NORMP
  266. SEGACT NORMP2
  267. JG=NTTDDL
  268. SEGINI,NORMP
  269. DO ITTDDL=1,NTTDDL
  270. NORMP.PROG(JPERM.LECT(ITTDDL))=NORMP2.PROG(ITTDDL)
  271. ENDDO
  272. SEGDES NORMP
  273. SEGSUP NORMP2
  274. ENDIF
  275. SEGSUP JPERM
  276. C SEGDES IDMATP
  277. SEGDES IDMATP
  278. MATRIK.KIDMAT(1)=IDMATP
  279. SEGDES MATRIK
  280. C
  281. C On stocke la factorisation obtenue du préconditionneur
  282. C
  283. SEGACT MATRIK*MOD
  284. KIDMAT(6)=ILUM
  285. KIDMAT(7)=ILUI
  286. SEGDES MATRIK
  287. C
  288. IF (IMPR.GT.6) THEN
  289. WRITE(IOIMP,*) 'création du préconditionneur Morse',
  290. $ ' de pointeurs',ILUM,'et',ILUI
  291. IF (IMPR.GT.8) THEN
  292. CALL ECMORS(ILUM,ILUI,(IMPR-1))
  293. ENDIF
  294. ENDIF
  295. ELSE
  296. IF (IMPR.GT.1) THEN
  297. WRITE(IOIMP,*)
  298. $ 'Préconditionneur deja construit.'
  299. ENDIF
  300. IF (IMPR.GT.6) THEN
  301. WRITE(IOIMP,*) 'Le préconditionneur est déjà construit :',
  302. $ 'ILUM=',ILUM,' et ILUI=',ILUI
  303. IF (IMPR.GT.8) THEN
  304. CALL ECMORS(ILUM,ILUI,(IMPR-1))
  305. ENDIF
  306. ENDIF
  307. ENDIF
  308. *
  309. * Normal termination
  310. *
  311. RETURN
  312. *
  313. * Format handling
  314. *
  315. *
  316. * Error handling
  317. *
  318. 9999 CONTINUE
  319. WRITE(IOIMP,*) 'An error was detected in priltp.eso'
  320. IRET=1
  321. RETURN
  322. *
  323. * End of PRILTP
  324. *
  325. END
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  

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