Télécharger kresll.eso

Retour à la liste

Numérotation des lignes :

  1. C KRESLL SOURCE PV 16/11/17 22:00:26 9180
  2. SUBROUTINE KRESLL()
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C*************************************************************************
  6. C Operateur KRES
  7. C
  8. C OBJET : Resoud l'equation AX=B par différentes méthodes
  9. C * directe (factorisation Choleski)
  10. C * itérative (Gradient conjugué, BiCGSTAB,
  11. C GMRES(m)) avec différents préconditionneurs
  12. C diagonal (Jacobi), D-ILU, ILU(0) (Choleski
  13. C incomplet), MILU(0)
  14. C SYNTAXE : CHPO3 = KRES MA1 'TYPI' TAB1
  15. C 'CLIM' CHPO1
  16. C 'SMBR' CHPO2
  17. C 'IMPR' VAL1 ;
  18. C Voir la notice pour plus de précisions.
  19. C
  20. C***********************************************************************
  21. C APPELES : KRES3, KRES4, KRES5
  22. C APPELES (E/S) : LIROBJ, ECROBJ, ERREUR, LIRMOT, LIRENT, LIRTAB,
  23. C ACME, ACMO, ACMM, ACMF, ECMO
  24. C APPELES (UTIL.) : QUETYP
  25. C APPELE PAR : KOPS
  26. C***********************************************************************
  27. C***********************************************************************
  28. C HISTORIQUE : 26/10/98 : prise en compte du cas particulier
  29. C où A est diagonale. C'est en fait assez pénible
  30. C car on utilise le CHPOINT comme structure de
  31. C données pour la matrice A et les vecteurs X,B,CLIM
  32. C HISTORIQUE : 09/02/99 : on peut utiliser le préconditionnement d'une
  33. C autre matrice que celle que l'on inverse
  34. C HISTORIQUE : 01/06/99 : on modifie la partie résolution itérative
  35. C en effet, lors de l'imposition des CL. de
  36. C Dirichlet, on changeait les valeurs de la
  37. C matrice Morse.
  38. C Ceci n'est pas bon lorsqu'on veut utiliser la
  39. C matrice assemblée pour plusieurs pas de temps.
  40. C On travaille maintenant sur une copie.
  41. C HISTORIQUE : 20/12/99 : reprogrammation de l'assemblage
  42. C Le nouvel assemblage n'est, pour l'instant effectif que
  43. C lorsqu'une méthode itérative est sélectionnée (-> fabrication
  44. C d'une matrice Morse). Le nouvel assemblage est plus performant
  45. C (temps de calcul, utilisation de la mémoire) et plus général (cas
  46. C particuliers à peu près supprimés) que le précédent.
  47. C HISTORIQUE : 05/01/00 : On ne supprime plus les 0.D0 de la matrice
  48. C assemblée (cf. clmors appelé par melim). Ceci pour ne plus avoir
  49. C perte éventuelle de symétrie de la matrice assemblée. Cela permet
  50. C aussi de ne plus dupliquer la matrice assemblée.
  51. C HISTORIQUE : 13/01/00 : Mise en conformité du solveur direct avec le
  52. C nouvel assemblage.
  53. C HISTORIQUE : 22/03/00 : Rajout des préconditionneurs ILUT
  54. C HISTORIQUE : 13/04/00 : Séparation Lecture des données
  55. C Ecriture des résultats (kres2)
  56. C Assemblage kres3
  57. C Méthode directe kres4
  58. C Méthodes itératives kres5
  59. C
  60. C***********************************************************************
  61. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  62. C en cas de modification de ce sous-programme afin de faciliter
  63. C la maintenance !
  64. C***********************************************************************
  65. -INC CCOPTIO
  66. -INC SMLREEL
  67. POINTEUR LRES.MLREEL
  68. -INC SMCHPOI
  69. POINTEUR KCLIM.MCHPOI
  70. POINTEUR KSMBR.MCHPOI
  71. POINTEUR MCHINI.MCHPOI
  72. POINTEUR MCHSOL.MCHPOI
  73. -INC SMTABLE
  74. POINTEUR MTINV.MTABLE
  75. * MATRIK est la matrice à inverser
  76. * MAPREC est la matrice dont on utilise le préconditionnement
  77. * MATASS est la matrice dont on utilise l'assemblage
  78. * pour préconditionner celui de MATRIK
  79. POINTEUR MAPREC.MATRIK
  80. POINTEUR MATASS.MATRIK
  81. C
  82. CHARACTER*8 TYPE
  83. * Paramètre m du GMRES(m)
  84. INTEGER IRSTRT
  85. INTEGER ITER
  86. REAL*8 BRTOL,RESID,RXMILU
  87. *
  88. REAL*8 XLFIL,XDTOL
  89. INTEGER KPREC
  90. INTEGER NMATRI
  91. INTEGER IP,KTYPI
  92. INTEGER IMPINV,IIMPR
  93. CHARACTER*4 MRENU,MMULAG
  94. INTEGER IMPR,IRET
  95. C
  96. C Définition des options
  97. C
  98. C 'IMPR' + un entier : niveau d'impression
  99. C 'TYPI' + une table contenant les informations
  100. C sur le type de résolution
  101. C 'CLIM' + un chpoint : conditions aux limites de Dirichlet
  102. C 'SMBR' + un chpoint : second membre (forces...)
  103. C 'PREC' + une matrice dont on va utiliser le préconditionnement
  104. C
  105. INTEGER NBM
  106. PARAMETER (NBM=4)
  107. CHARACTER*4 LMOT(NBM)
  108. C Tableau de correspondance : entier <-> mot
  109. C pour le type d'inversion
  110. INTEGER NBTYPI
  111. C INTEGER LNTYPI
  112. PARAMETER (NBTYPI=5)
  113. C PARAMETER (LNTYPI=18)
  114. C CHARACTER*(LNTYPI) LTYPI(NBTYPI)
  115. LOGICAL LBID,LVAL,LTIME,LMRENU,LDEPE,LLDEPE
  116. C
  117. C Initialisation des tableaux
  118. C
  119. DATA LMOT/'IMPR','TYPI','CLIM','SMBR'/
  120. C DATA LTYPI/'Choleski ',
  121. C $ 'Conjugate Gradient',
  122. C $ 'BiCGSTAB ',
  123. C $ 'BiCGSTAB(2) ',
  124. C $ 'GMRES(m) '/
  125. C
  126. TYPE='MATRIK'
  127. CALL LIROBJ(TYPE,MATRIK,1,IRET)
  128. IF(IRET.EQ.0) GOTO 9999
  129. c
  130. SEGACT MATRIK
  131. NMATRI=IRIGEL(/2)
  132. IF(NMATRI.EQ.0)THEN
  133. C% Résolution impossible : la matrice de RIGIDITE est vide
  134. CALL ERREUR(727)
  135. RETURN
  136. ENDIF
  137. SEGDES MATRIK
  138. c Valeurs par défaut reprise de prkres
  139. KSMBR=0
  140. MTINV=0
  141. IMPR=0
  142. KCLIM=0
  143. KTYPI=1
  144. MATASS=MATRIK
  145. MAPREC=MATRIK
  146. MRENU='SLOA'
  147. LMRENU=.FALSE.
  148. MMULAG='APR2'
  149. ISCAL=0
  150. IOUBL=0
  151. IMPINV=0
  152. MCHINI=0
  153. ITER=2000
  154. RESID=1.D-10
  155. BRTOL=1.D-40
  156. IRSTRT=50
  157. KPREC=3
  158. RXMILU=1.D0
  159. RXILUP=0.5D0
  160. XLFIL=2.D0
  161. XDTOL=-1.D0
  162. XSPIV=0.1D0
  163. LBCG=4
  164. ICALRS=0
  165. METASS=5
  166. LTIME=.FALSE.
  167. KTIME=0
  168. LDEPE=.TRUE.
  169. LLDEPE=.FALSE.
  170. IDDOT=0
  171. IMVEC=2
  172. C
  173. 1 CONTINUE
  174. CALL LIRMOT(LMOT,NBM,IP,0)
  175. IF(IP.EQ.0) GOTO 2
  176. GOTO (11,12,13,14),IP
  177. C Lecture du niveau d'impression
  178. C option 'IMPR'
  179. 11 CONTINUE
  180. c write(6,*)' option IMPR'
  181. CALL LIRENT(IMPR,1,IRET)
  182. IF (IRET.EQ.0) GOTO 9999
  183. GO TO 1
  184. C Lecture de la table pour les méthodes d'inversion
  185. C option 'TYPI'
  186. 12 CONTINUE
  187. CALL LIRTAB('METHINV',MTINV,1,IRET)
  188. IF(IRET.EQ.0)RETURN
  189. CALL LIRENT(KTYPI,1,IRET)
  190. C KTYPI : Méthode d'inversion du système (cf. LTYPI)
  191. C 1 : résolution directe (Choleski)
  192. C 2 : Gradient Conjugué
  193. C 3 : Bi-Gradient Conjugué Stabilisé (BiCGSTAB)
  194. C 4 : BiCGSTAB(2)
  195. C 5 : GMRES(m)
  196. C CALL ACME(MTINV,'TYPINV',KTYPI)
  197. IF (IERR.NE.0) GOTO 9999
  198. C write(6,*)' Ktypi=',ktypi
  199. IF (KTYPI.LT.1.OR.KTYPI.GT.NBTYPI) THEN
  200. WRITE(IOIMP,*) 'L''indice TYPINV ',
  201. $ 'is out of range (1..',NBTYPI,') :',KTYPI
  202. GOTO 9999
  203. ENDIF
  204. GO TO 1
  205. C Lecture du chpoint de conditions aux limites
  206. C option CLIM
  207. 13 CONTINUE
  208. CALL QUETYP(TYPE,1,IRET)
  209. IF(IRET.EQ.0) GOTO 9999
  210. IF(TYPE.EQ.'CHPOINT')THEN
  211. CALL LIROBJ(TYPE,KCLIM,1,IRET)
  212. ELSE
  213. WRITE(IOIMP,*) 'Erreur lecture C.lim'
  214. GOTO 9999
  215. ENDIF
  216. GO TO 1
  217. C Lecture du chpoint second membre
  218. C option SMBR
  219. 14 CONTINUE
  220. CALL QUETYP(TYPE,1,IRET)
  221. IF(IRET.EQ.0)RETURN
  222. IF(TYPE.EQ.'CHPOINT')THEN
  223. CALL LIROBJ('CHPOINT',KSMBR,1,IRET)
  224. c write(6,*)' lecture 2 membre '
  225. ELSE
  226. WRITE(IOIMP,*) 'Erreur lecture second membre'
  227. GOTO 9999
  228. ENDIF
  229. GOTO 1
  230. 2 CONTINUE
  231. C write(6,*)' On a pas finit les lectures '
  232. c NAT=2
  233. c NSOUPO=0
  234. c SEGINI MCHPO1
  235. c MCHPO1.JATTRI(1)=2
  236. c CALL ECROBJ('CHPOINT',MCHPO1)
  237. c return
  238. C
  239. C Matrice dont on utilise l'assemblage
  240. C (éventuellement identique à la matrice transmise)
  241. C
  242. TYPE='MATRIK '
  243. C CALL ACMO(MTINV,'MATASS',TYPE,MATASS)
  244. C IF (IERR.NE.0) GOTO 9999
  245. CALL LIROBJ('MATRIK',MATASS,1,IRET)
  246. IF(IRET.EQ.0) GOTO 9999
  247. C - Pour l'assemblage : type de renumérotation
  248. C CALL ACMM(MTINV,'TYRENU',MRENU)
  249. C IF (IERR.NE.0) GOTO 9999
  250. CALL LIRCHA(MRENU,1,LCHAR)
  251. IF(LCHAR.EQ.0)GO TO 9999
  252. C - Pour l'assemblage : prise en compte des mult.lag
  253. C CALL ACMM(MTINV,'PCMLAG',MMULAG)
  254. C IF (IERR.NE.0) GOTO 9999
  255. CALL LIRCHA(MMULAG,1,LCHAR)
  256. IF(LCHAR.EQ.0)GO TO 9999
  257. C - Pour l'assemblage : scaling de la matrice
  258. C CALL ACME(MTINV,'SCALING',ISCAL)
  259. C IF (IERR.NE.0) GOTO 9999
  260. ISCAL=0
  261. C Niveau d'impression pour la partie résolution itérative
  262. c CALL ACME(MTINV,'IMPINV',IMPINV)
  263. c IF (IERR.NE.0) GOTO 9999
  264. c IMPR=MAX(IMPR,IMPINV)
  265. IMPR=0
  266. C
  267. C Assemblage proprement dit
  268. C
  269. IIMPR=0
  270. * SG 2016/02/09 : non à la ligne suivante : il faut que METASS soit
  271. * égale à 5 (voir remarque dans makpr2)
  272. * METASS=4
  273. CALL KRES3(MATRIK,MATASS,MRENU,MMULAG,METASS,
  274. $ KTIME,LTIME,
  275. $ IIMPR,IRET)
  276. IF (IRET.NE.0) GOTO 9999
  277. *! WRITE(IOIMP,*) 'Aprés assemblage'
  278. *! CALL ECRENT(5)
  279. *! CALL ECROBJ('MATRIK ',MATRIK)
  280. *! CALL PRLIST
  281. *
  282. * "Oubli" des valeurs des matrice élémentaires
  283. * On met les tableaux de LIZAFM à 0 => à MENAGE de les supprimmer
  284. * si besoin est.
  285. *
  286. CALL OUBIMA(MATRIK,IMPR,IRET)
  287. IF (IRET.NE.0) GOTO 9999
  288. C
  289. C Méthode directe
  290. C
  291. IF (KTYPI.EQ.1) THEN
  292. CALL KRES4(MATRIK,KCLIM,KSMBR,
  293. $ ISCAL,
  294. $ MCHSOL,
  295. $ IMPR,IRET)
  296. IF (IRET.NE.0) GOTO 9999
  297. C
  298. C Méthodes itératives
  299. C
  300. ELSEIF(KTYPI.GE.2.AND.KTYPI.LE.NBTYPI)THEN
  301. C Paramètres pour la méthode itérative proprement dite
  302. C - Champoint d'initialisation de la méthode
  303. C (i.e. estimation de l'inconnue)
  304. TYPE='CHPOINT '
  305. C CALL ACMO(MTINV,'XINIT',TYPE,MCHINI)
  306. C IF (IERR.NE.0) GOTO 9999
  307. CALL LIROBJ('CHPOINT',MCHINI,1,IRET)
  308. IF(IRET.EQ.0) GOTO 9999
  309. C - Nombre maxi d'itérations à effectuer
  310. C CALL ACME(MTINV,'NITMAX',ITER)
  311. C IF (IERR.NE.0) GOTO 9999
  312. CALL LIRENT(ITER,1,IRET)
  313. IF(IRET.EQ.0) GOTO 9999
  314. C - Norme maxi (L2 normé par le second membre) du résidu
  315. C CALL ACMF(MTINV,'RESID',RESID)
  316. C IF (IERR.NE.0) GOTO 9999
  317. CALL LIRREE(RESID,1,IRET)
  318. IF(IRET.EQ.0) GOTO 9999
  319. C - 'Breakdown tolerance' pour les méthodes de type
  320. C BiCGSTAB. Si un certain produit scalaire de vecteurs
  321. C "direction" est inférieur à cette tolérance, la
  322. C méthode s'arrete.
  323. IF (KTYPI.EQ.3.OR.KTYPI.EQ.4) THEN
  324. C CALL ACMF(MTINV,'BCGSBTOL',BRTOL)
  325. C IF (IERR.NE.0) GOTO 9999
  326. CALL LIRREE(BRTOL,1,IRET)
  327. IF(IRET.EQ.0) GOTO 9999
  328. ENDIF
  329. C - Paramètre m de redémarrage pour GMRES(m)
  330. C i.e. nombre max. de vecteurs par rapport auxquels on
  331. C orthogonalise le résidu
  332. IF (KTYPI.EQ.5) THEN
  333. C CALL ACME(MTINV,'GMRESTRT',IRSTRT)
  334. C IF (IERR.NE.0) GOTO 9999
  335. CALL LIRENT(IRSTRT,1,IRET)
  336. IF(IRET.EQ.0) GOTO 9999
  337. ENDIF
  338. C Paramètres pour les différents préconditionneurs
  339. C Matrice dont on utilise le préconditionneur
  340. C (éventuellement identique à la matrice transmise).
  341. C TYPE='MATRIK '
  342. C CALL ACMO(MTINV,'MAPREC',TYPE,MAPREC)
  343. C IF (IERR.NE.0) GOTO 9999
  344. CALL LIROBJ('MATRIK',MAPREC,1,IRET)
  345. IF(IRET.EQ.0) GOTO 9999
  346. C - Type de préconditionnement : (cf. LPREC)
  347. C 0 : pas de préconditionnement
  348. C 1 : préconditionnement par la diagonale
  349. C 2 : préconditionnement D-ILU
  350. C 3 : préconditionnement ILU(0) (Choleski)
  351. C 4 : préconditionnement MILU(0) (Choleski modifié)
  352. C 5 : préconditionnement ILUT (dual truncation)
  353. C 6 : préconditionnement ILUT2 (une variante du
  354. C précédent qui remplit mieux la mémoire et
  355. C fonctionne mieux quelquefois)
  356. C CALL ACME(MTINV,'PRECOND',KPREC)
  357. C IF (IERR.NE.0) GOTO 9999
  358. CALL LIRENT(KPREC,1,IRET)
  359. IF(IRET.EQ.0) GOTO 9999
  360. C - Paramètre de relaxation pour le préconditionnement
  361. C MILU(0) compris entre 0.D0 et 1.D0
  362. C S'il est égal à 0, on se ramène à ILU(0)
  363. C S'il est égal à 1, MILU(0) est dit non relaxé
  364. IF (KPREC.EQ.4) THEN
  365. C CALL ACMF(MTINV,'MILURELX',RXMILU)
  366. C IF (IERR.NE.0) GOTO 9999
  367. CALL LIRREE(RXMILU,1,IRET)
  368. IF(IRET.EQ.0) GOTO 9999
  369. ENDIF
  370. C - Pour une méthode ILUT, on a les deux indices suivant :
  371. C * ILUTLFIL : encombrement maximal (approximatif) du
  372. C préconditionneur, par rapport à la matrice.
  373. C * ILUTDTOL : "drop tolerance" pour le préconditionneur.
  374. C i.e. en-dessous de cette valeur relative, les
  375. C termes de la factorisation incomplète seront
  376. C oubliés.
  377. IF (KPREC.GE.5.AND.KPREC.LE.9) THEN
  378. C CALL ACMF(MTINV,'ILUTLFIL',XLFIL)
  379. C IF (IERR.NE.0) GOTO 9999
  380. CALL LIRREE(XLFIL,1,IRET)
  381. IF(IRET.EQ.0) GOTO 9999
  382. C CALL ACMF(MTINV,'ILUTDTOL',XDTOL)
  383. C IF (IERR.NE.0) GOTO 9999
  384. CALL LIRREE(XDTOL,1,IRET)
  385. IF(IRET.EQ.0) GOTO 9999
  386. ENDIF
  387. C - Pour une méthode ILUTP, on a les deux indices suivant :
  388. C * ILUTPPIV (type REEL) : sensibilité du pivoting
  389. C IF (KPREC.GE.7.AND.KPREC.LE.9) THEN
  390. C CALL ACMF(MTINV,'ILUTPPIV',XSPIV)
  391. C ENDIF
  392. CALL KRES5(MATRIK,KCLIM,KSMBR,KTYPI,
  393. $ MCHINI,ITER,RESID,
  394. $ BRTOL,IRSTRT,LBCG,ICALRS,
  395. $ MAPREC,KPREC,
  396. $ RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,
  397. $ ISCAL,
  398. $ KTIME,LTIME,IDDOT,IMVEC,
  399. $ MCHSOL,LRES,LNMV,ICVG,
  400. $ IMPR,IRET)
  401. IF (IRET.NE.0) GOTO 9999
  402. CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES)
  403. ELSE
  404. WRITE(IOIMP,*) 'KTYPI=',KTYPI,'out of range (1..',NBTYPI,')'
  405. GOTO 9999
  406. ENDIF
  407. CALL ECROBJ('CHPOINT ',MCHSOL)
  408. *
  409. * Normal termination
  410. *
  411. RETURN
  412. *
  413. * Format handling
  414. *
  415. *
  416. * Error handling
  417. *
  418. 9999 CONTINUE
  419. WRITE(IOIMP,*) 'An error was detected in kresll.eso'
  420. * 153 2
  421. * Opération illicite dans ce contexte
  422. CALL ERREUR(153)
  423. RETURN
  424. *
  425. * End of KRES2
  426. *
  427. END
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  

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