Télécharger kres2.eso

Retour à la liste

Numérotation des lignes :

  1. C KRES2 SOURCE PV 16/11/17 22:00:23 9180
  2. SUBROUTINE KRES2()
  3. IMPLICIT REAL*8 (A-H,O-Z)
  4. IMPLICIT INTEGER (I-N)
  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 HISTORIQUE : 06/04/04 : Renumérotation (double mult.)
  60. C HISTORIQUE : 06/04/04 : Scaling
  61. C HISTORIQUE : 08/04/04 : ajout ILUTP
  62. C HISTORIQUE : 09/02/06 : ajout nb prod matrice vecteur (NMATVEC)
  63. C simplification du code
  64. C HISTORIQUE : 22/02/06 : Gros nettoyage au niveau de l'entrée des
  65. C arguments (Nouvelle syntaxe)
  66. C HISTORIQUE : 08/2011 : En vue de la suppression de l'objet MATRIK
  67. C on utilise l'assemblage de RESOU.
  68. C
  69. C***********************************************************************
  70. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  71. C en cas de modification de ce sous-programme afin de faciliter
  72. C la maintenance !
  73. C***********************************************************************
  74. -INC CCOPTIO
  75. -INC SMLREEL
  76. POINTEUR LRES.MLREEL
  77. -INC SMLENTI
  78. POINTEUR LNMV.MLENTI
  79. -INC SMCHPOI
  80. POINTEUR KCLIM.MCHPOI
  81. POINTEUR KSMBR.MCHPOI
  82. POINTEUR MCHINI.MCHPOI
  83. POINTEUR MCHSOL.MCHPOI
  84. -INC SMTABLE
  85. POINTEUR MTINV.MTABLE
  86. POINTEUR KTIME.MTABLE
  87. DIMENSION ITTIME(4)
  88. CHARACTER*8 CHARI
  89. * MATRIK est la matrice à inverser
  90. * MAPREC est la matrice dont on utilise le préconditionnement
  91. * MATASS est la matrice dont on utilise l'assemblage
  92. * pour préconditionner celui de MATRIK
  93. POINTEUR MAPREC.MATRIK
  94. POINTEUR MATASS.MATRIK
  95. *STAT -INC SMSTAT
  96. C
  97. CHARACTER*8 TYPE
  98. * Paramètre m du GMRES(m)
  99. INTEGER RESTRT
  100. INTEGER ITER
  101. REAL*8 BRTOL,RESID,RXMILU,RXILUP
  102. *
  103. REAL*8 XLFIL,XDTOL
  104. INTEGER KPREC
  105. INTEGER NMATRI
  106. INTEGER IP,KTYPI
  107. INTEGER IMPINV,IIMPR
  108. CHARACTER*4 MRENU,MMULAG
  109. LOGICAL LRIG,LTIME,LDETR,LDEPE,LASRIG,LDMULT,LOGII
  110. INTEGER IMPR,IRET
  111. C
  112. C Lecture des arguments et mise à défaut des optionnels ()
  113. C
  114. C MATRIK : La matrice lue en entrée au format MATRIK
  115. C MTINV : L'éventuelle table d'inversion (obsolète)
  116. C IMPR : Niveau d'impression solveur direct
  117. C KCLIM : Chpoint éventuel de conditions aux limites de Dirichlet
  118. C KSMBR : Chpoint second membre
  119. C KTYPI : Type de méthode de résolution
  120. C MATASS : Matrice utilisée pour préconditionner l'assemblage
  121. C MAPREC : Matrice utilisée pour préconditionner la construction du
  122. C préconditionneur
  123. C MRENU : Type de renumérotation
  124. C MMULAG : Placement des multiplicateurs de Lagrange
  125. C ISCAL : Scaling de la matrice
  126. C IOUBL : Oubli des matrices élémentaires ?
  127. C IMPINV : Niveau d'impression solveur itératif
  128. C MCHINI : Chpoint estimation de l'inconnue
  129. C ITER : Nombre maxi d'itérations à effectuer
  130. C RESID : Norme L2 maxi du résidu
  131. C BRTOL : Breakdown tolerance pour les méthodes de type Bi-CG
  132. C IRSTRT : Paramètre m de redémarrage pour GMRES
  133. C KPREC : Type du préconditionneur
  134. C RXMILU : Paramètre de relaxation pour MILU(0)
  135. C RXILUP : Paramètre de filtre pour ILU(0)-PV
  136. C XLFIL : Paramètre de remplissage pour les préconditionneurs ILUT
  137. C XDTOL : Drop tolerance pour les préconditionneurs ILUT
  138. C XSPIV : Sensibilité du pivoting pour les préconditionneurs ILUTP
  139. C LBCG : le l dans BicgStab(l)
  140. C ICALRS : façon de calculer le résidu
  141. C METASS : Méthode d'assemblage
  142. C LTIME : construit une table avec des statistiques temporelles
  143. C si vrai
  144. C LDEPE : élimine les dépendances si VRAI
  145. C et matrice d'entrée RIGIDITE
  146. C IDDOT : 0 => utilisation du produit scalaire normal dans les
  147. C méthodes itératives
  148. C 1 => utilisation du produit scalaire compensé
  149. * IMPR=4
  150.  
  151. IVALI=0
  152. XVALI=0.D0
  153. LOGII=.FALSE.
  154. IRETI=0
  155. XVALR=0.D0
  156. IOBRE=0
  157. IRETR=0
  158.  
  159. IMPR=0
  160. * WRITE(IOIMP,*) 'coucou kres2'
  161. *
  162. *STAT CALL INMSTA(MSTAT,IMPR)
  163. *
  164. CALL PRKRES(MATRIK,MTINV,IMPR,KCLIM,KSMBR,KTYPI,MATASS,MAPREC,
  165. $ MRENU,MMULAG,ISCAL,IOUBL,IMPINV,MCHINI,ITER,RESID,BRTOL,
  166. $ IRSTRT,KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,LBCG,ICALRS,
  167. $ METASS,LTIME,LDEPE,MRIGID,IDDOT,IMVEC,IRET)
  168. IF (IRET.NE.0) GOTO 9999
  169. IMPR=MAX(IMPR,IMPINV)
  170. *
  171. * Préparation de la matrice et du second membre
  172. * suivant les cas
  173. *
  174. * LASRIG=.TRUE. on utilise l'assemblage de RESOU
  175. * LDMULT=.TRUE. on dédouble les multiplicateurs de Lagrange
  176. *
  177.  
  178. LASRIG=(METASS.EQ.6)
  179. LDMULT=LASRIG
  180. MRIGI0=MRIGID
  181. IF (MRIGID.NE.0) THEN
  182. IF (LDEPE) THEN
  183. KSMBR0=KSMBR
  184. * CALL KRES6(MRIGID,KSMBR,IDEPE,
  185. * $ MRIGIC,KSMBRC,KSMBR0,KSMBR1)
  186. CALL KRES6(MRIGID,KSMBR,LDMULT,
  187. $ MRIGIC,KSMBRC,KSMBR1)
  188. IF (IERR.NE.0) RETURN
  189. MRIGID=MRIGIC
  190. KSMBR=KSMBRC
  191. ENDIF
  192. IF (LASRIG) THEN
  193. * Gestion de la normalisation
  194. NORICO=NORINC
  195. NORIDO=NORIND
  196. IF (ISCAL.EQ.0) THEN
  197. NORINC=0
  198. NORIND=0
  199. ELSE
  200. NORINC=-1
  201. NORIND=0
  202. ENDIF
  203. * Gestion de la renumérotation
  204. NUCROO=NUCROU
  205. IF (MRENU.EQ.'RIEN') THEN
  206. NUCROU=-1
  207. ELSEIF (MRENU.EQ.'SLOA') THEN
  208. NUCROU=1
  209. ELSEIF (MRENU.EQ.'GIPR'.OR.MRENU.EQ.'GIBA') THEN
  210. NUCROU=2
  211. ELSE
  212. WRITE(IOIMP,*) 'MRENU=',MRENU
  213. CALL ERREUR(5)
  214. RETURN
  215. ENDIF
  216. CALL KRES8(MRIGID,KSMBR,
  217. $ KTYPI,ITER,RESID,ICALRS,IRSTRT,LBCG,BRTOL,IDDOT,IMVEC,
  218. $ KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,
  219. $ KTIME,LTIME,
  220. $ MCHSOL,LRES,LNMV,ICVG,IMPR)
  221. IF (IERR.NE.0) RETURN
  222. * Gestion de la normalisation
  223. NORINC=NORICO
  224. NORIND=NORIDO
  225. NUCROU=NUCROO
  226. IF (LTIME) CALL ECROBJ('TABLE ',KTIME)
  227. IF (MTINV.NE.0) THEN
  228. CALL ECMO(MTINV,'CVGOK','ENTIER',ICVG)
  229. IF (LRES.NE.0) CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES)
  230. IF (LNMV.NE.0) CALL ECMO(MTINV,'NMATVEC','LISTENTI',LNMV)
  231. ENDIF
  232. * On décondense si nécessaire
  233. *
  234. * write (6,*) ' resou - mchsol '
  235. * call ecchpo(mchsol,0)
  236. * call mucpri(mchsol,mrigid,iresi)
  237. * write (6,*) ' kres - iresi '
  238. * call ecchpo(iresi,0)
  239. * WRITE(IOIMP,*) 'Avant KRES7'
  240. IF (MRIGI0.NE.0.AND.LDEPE) THEN
  241. MSOLC=MCHSOL
  242. CALL KRES7(MSOLC,MRIGI0,KSMBR0,KSMBR1,1,
  243. $ MCHSOL)
  244. IF (IERR.NE.0) RETURN
  245. ENDIF
  246. CALL ECROBJ('CHPOINT ',MCHSOL)
  247. RETURN
  248. ELSE
  249. CALL ECROBJ('RIGIDITE',MRIGID)
  250. CALL RIMA
  251. IF (IERR.NE.0) GOTO 9999
  252. CALL MACHI2(1)
  253. IF (IERR.NE.0) GOTO 9999
  254. CALL LIROBJ('MATRIK',MATRIK,1,IRET)
  255. IF(IRET.EQ.0) GOTO 9999
  256.  
  257. * Changement des noms d'inconnues du second membre
  258. IF (KSMBR.NE.0) THEN
  259. CALL ECROBJ('CHPOINT ',KSMBR)
  260. CALL MACHI2(1)
  261. CALL LIROBJ('CHPOINT ',KSMBR,1,IRET)
  262. IF (IRET.EQ.0) GOTO 9999
  263. ENDIF
  264. ENDIF
  265. * write (6,*) ' le vecteur 2'
  266. * call ecchpo(ksmbr,0)
  267. * write (6,*) ' la matrice 2'
  268. * call ecrent(5)
  269. * call ecmatk(matrik)
  270. ENDIF
  271. *
  272. SEGACT MATRIK
  273. NMATRI=IRIGEL(/2)
  274. IF(NMATRI.EQ.0)THEN
  275. C% Résolution impossible : la matrice de RIGIDITE est vide
  276. CALL ERREUR(727)
  277. RETURN
  278. ENDIF
  279. SEGDES MATRIK
  280. IF (MATASS.EQ.0) MATASS=MATRIK
  281. IF (MAPREC.EQ.0) MAPREC=MATRIK
  282. * WRITE(IOIMP,*) 'Sortie de prkres'
  283. * WRITE(IOIMP,*) 'IOUBL=',IOUBL
  284. C
  285. IF (LTIME) THEN
  286. CALL CRTABL(KTIME)
  287. CALL TIMESPV(ITTIME)
  288. ITI1=(ITTIME(1)+ITTIME(2))/10
  289. ELSE
  290. KTIME=0
  291. ENDIF
  292. *STAT CALL PRMSTA('Lectures',MSTAT,IMPR)
  293. *
  294. C
  295. C Assemblage proprement dit
  296. C
  297. IIMPR=0
  298. CALL KRES3(MATRIK,MATASS,MRENU,MMULAG,METASS,
  299. $ KTIME,LTIME,
  300. $ IIMPR,IRET)
  301. IF (IRET.NE.0) GOTO 9999
  302. *! WRITE(IOIMP,*) 'Aprés assemblage'
  303. *! CALL ECRENT(5)
  304. *! CALL ECROBJ('MATRIK ',MATRIK)
  305. *! CALL PRLIST
  306. IF (LTIME) THEN
  307. CALL TIMESPV(ITTIME)
  308. ITI2=(ITTIME(1)+ITTIME(2))/10
  309. ENDIF
  310. *STAT CALL PRMSTA('Assemblage',MSTAT,IMPR)
  311. *
  312. * "Oubli" des valeurs des matrice élémentaires
  313. * On met les tableaux de LIZAFM à 0 => à MENAGE de les supprimmer
  314. * si besoin est.
  315. *
  316. IOUBD=MOD(IOUBL,10)
  317. *! WRITE(IOIMP,*) 'IOUBD=',IOUBD
  318. IF (IOUBD.EQ.1) THEN
  319. CALL OUBIMA(MATRIK,IMPR,IRET)
  320. IF (IRET.NE.0) GOTO 9999
  321. IF (IMPR.GT.2) THEN
  322. WRITE(IOIMP,*) 'Oubli des mat. elem.'
  323. ENDIF
  324. ELSEIF (IOUBD.EQ.2) THEN
  325. ith=0
  326. call ooonth(ith)
  327. call ooohor(0,ith)
  328. SEGACT MATRIK*MOD
  329. LDETR=.FALSE.
  330. NMATE=IRIGEL(/2)
  331. DO IMATE=1,NMATE
  332. IMATRI=IRIGEL(4,IMATE)
  333. SEGACT IMATRI*MOD
  334. NSOUM =LIZAFM(/1)
  335. NTOTIN=LIZAFM(/2)
  336. DO ITOTIN=1,NTOTIN
  337. DO ISOUM=1,NSOUM
  338. IZAFM=LIZAFM(ISOUM,ITOTIN)
  339. IF (IZAFM.NE.0) THEN
  340. LDETR=.TRUE.
  341. SEGSUP IZAFM
  342. LIZAFM(ISOUM,ITOTIN)=0
  343. ENDIF
  344. ENDDO
  345. ENDDO
  346. SEGDES IMATRI
  347. ENDDO
  348. IF (IMPR.GT.2.AND.LDETR) THEN
  349. WRITE(IOIMP,*) 'Destruction des mat. elem.'
  350. ENDIF
  351. ELSEIF (IOUBD.NE.0) THEN
  352. WRITE(IOIMP,*) 'IOUBL=',IOUBL, ' non prevu'
  353. GOTO 9999
  354. ENDIF
  355. *STAT CALL PRMSTA('Oubli',MSTAT,IMPR)
  356. *! WRITE(IOIMP,*) 'Aprés oubli'
  357. C
  358. C Méthode directe
  359. C
  360. IF (KTYPI.EQ.1) THEN
  361. CALL KRES4(MATRIK,KCLIM,KSMBR,
  362. $ ISCAL,
  363. $ MCHSOL,
  364. $ IMPR,IRET)
  365. IF (IRET.NE.0) GOTO 9999
  366. *STAT CALL PRMSTA('Methode directe',MSTAT,IMPR)
  367. C
  368. C Méthodes itératives
  369. C
  370. ELSE
  371. C
  372. CALL KRES5(MATRIK,KCLIM,KSMBR,KTYPI,
  373. $ MCHINI,ITER,RESID,
  374. $ BRTOL,IRSTRT,LBCG,ICALRS,
  375. $ MAPREC,KPREC,
  376. $ RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,
  377. $ ISCAL,
  378. $ KTIME,LTIME,IDDOT,IMVEC,
  379. $ MCHSOL,LRES,LNMV,ICVG,
  380. $ IMPR,IRET)
  381. IF (IRET.NE.0) GOTO 9999
  382. *STAT CALL PRMSTA('Methode iterative',MSTAT,IMPR)
  383. IF (MTINV.NE.0) THEN
  384. CALL ECMO(MTINV,'CVGOK','ENTIER',ICVG)
  385. CALL ECMO(MTINV,'CONVINV','LISTREEL',LRES)
  386. CALL ECMO(MTINV,'NMATVEC','LISTENTI',LNMV)
  387. ENDIF
  388. ENDIF
  389. IF (LTIME) THEN
  390. CALL TIMESPV(ITTIME)
  391. ITI3=(ITTIME(1)+ITTIME(2))/10
  392. CHARI='ASS+RENU'
  393. CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI,
  394. $ 'ENTIER ',ITI2-ITI1,XVALR,CHARR,LOGIR,IRETR)
  395. CHARI='PRE+RESO'
  396. CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI,
  397. $ 'ENTIER ',ITI3-ITI2,XVALR,CHARR,LOGIR,IRETR)
  398. CHARI='TOTAL '
  399. CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI,
  400. $ 'ENTIER ',ITI3-ITI1,XVALR,CHARR,LOGIR,IRETR)
  401. SEGDES KTIME
  402. CALL ECROBJ('TABLE ',KTIME)
  403. ENDIF
  404. IOUBE=IOUBL/10
  405. *! WRITE(IOIMP,*) 'IOUBE=',IOUBE
  406. IF (IOUBE.GE.1) THEN
  407. ith=0
  408. call ooonth(ith)
  409. call ooohor(0,ith)
  410. SEGACT MATRIK*MOD
  411. IF (IOUBE.EQ.2) THEN
  412. PMORS=KIDMAT(4)
  413. IF (PMORS.NE.0) THEN
  414. IF (IMPR.GT.2) THEN
  415. WRITE(IOIMP,*) 'Destruction du profil morse'
  416. ENDIF
  417. SEGSUP PMORS
  418. KIDMAT(4)=0
  419. ENDIF
  420. ENDIF
  421. IZA=KIDMAT(5)
  422. IF (IZA.NE.0) THEN
  423. IF (IMPR.GT.2) THEN
  424. WRITE(IOIMP,*) 'Destruction des valeurs morses'
  425. ENDIF
  426. SEGSUP IZA
  427. KIDMAT(5)=0
  428. ENDIF
  429. PMORS=KIDMAT(6)
  430. IF (PMORS.NE.0) THEN
  431. IF (IMPR.GT.2) THEN
  432. WRITE(IOIMP,*) 'Destruction du profil du precon'
  433. ENDIF
  434. SEGSUP PMORS
  435. KIDMAT(6)=0
  436. ENDIF
  437. IZA=KIDMAT(7)
  438. IF (IZA.NE.0) THEN
  439. IF (IMPR.GT.2) THEN
  440. WRITE(IOIMP,*) 'Destruction des valeurs du precon'
  441. ENDIF
  442. SEGSUP IZA
  443. KIDMAT(7)=0
  444. ENDIF
  445. SEGDES MATRIK
  446. ELSEIF (IOUBE.NE.0) THEN
  447. WRITE(IOIMP,*) 'IOUBL=',IOUBL, ' non prevu'
  448. GOTO 9999
  449. ENDIF
  450. *
  451. * On décondense si nécessaire
  452. *
  453. * write (6,*) ' resou - mchsol '
  454. * call ecchpo(mchsol,0)
  455. * call mucpri(mchsol,mrigid,iresi)
  456. * write (6,*) ' kres - iresi '
  457. * call ecchpo(iresi,0)
  458. * WRITE(IOIMP,*) 'Avant KRES7'
  459. IF (MRIGI0.NE.0.AND.LDEPE) THEN
  460. MSOLC=MCHSOL
  461. * CALL KRES7(MRIGID,IDEPE,KSMBR0,KSMBR1,
  462. * $ MSOLC,
  463. * $ MCHSOL)
  464. CALL KRES7(MSOLC,MRIGI0,KSMBR0,KSMBR1,1,
  465. $ MCHSOL)
  466. IF (IERR.NE.0) RETURN
  467. ENDIF
  468. CALL ECROBJ('CHPOINT ',MCHSOL)
  469. *STAT CALL SUMSTA(MSTAT,IMPR)
  470. *
  471. * Normal termination
  472. *
  473. RETURN
  474. *
  475. * Format handling
  476. *
  477. *
  478. * Error handling
  479. *
  480. 9999 CONTINUE
  481. WRITE(IOIMP,*) 'An error was detected in kres2.eso'
  482. * 153 2
  483. * Opération illicite dans ce contexte
  484. CALL ERREUR(153)
  485. RETURN
  486. *
  487. * End of KRES2
  488. *
  489. END
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  
  504.  
  505.  
  506.  
  507.  
  508.  
  509.  
  510.  
  511.  
  512.  
  513.  
  514.  
  515.  
  516.  
  517.  
  518.  
  519.  
  520.  

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