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

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