Télécharger cgresd.eso

Retour à la liste

Numérotation des lignes :

  1. C CGRESD SOURCE CHAT 05/01/12 21:53:31 5004
  2. SUBROUTINE CGRESD(KSTO,
  3. 1 NL,ILG,
  4. 2 IMAT,IPRC,IA,JA,KA,
  5. 3 IZB,IZP,ICOLD,
  6. 4 NPT,NPITE,NEFF,ICONV,EPI,IPOU,VPOU,
  7. 5 NIMPR,IPAT)
  8. C
  9. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  10. C C
  11. C RESOLUTION ITERATIVE D'UN SYSTEME C
  12. C SYMETRIQUE CREUX PAR LA METHODE DU C
  13. C GRADIENT CONJUGUE PRECONDITIONNE PAR LA C
  14. C DIAGONALE. C
  15. C C
  16. C***** INFOS SUR LA METHODE (LIGNE 0) C
  17. C C
  18. C DEUX MODES DE STOCKAGE SONT PROPOSES C
  19. C MORSE OU COMPRESSE. C
  20. C C
  21. C MORSE (KSTO=0) C
  22. C COMPRESSE (KSTO=1) C
  23. C C
  24. C POUR LE CAS DCG, ON RESOUD L'EQUATION SUIVANTE C
  25. C C
  26. C ( (K-1)*A*(K-1) )*xi = (K-1)*b C
  27. C C
  28. C AVEC K = D**(1/2) C
  29. C pi est la pseudo solution C
  30. C La solution x est obtenue en faisant C
  31. C x = (K-1)*xi C
  32. C C
  33. C***** INFOS SUR LES DIMENSIONS (LIGNE 1) C
  34. C C
  35. C NL (L'ORDRE DU SYSTEME) C
  36. C C
  37. C ILG (LONGUEUR LA POUR LE MORSE) C
  38. C (NOMBRE MAXI DE VALEURS PAR LIGNE DE LA C
  39. C MATRICE POUR LE COMPRESSE) C
  40. C C
  41. C***** INFOS SUR LA MATRICE (LIGNE 2) C
  42. C C
  43. C IMAT (LA MATRICE DANS UN LISTREEL) C
  44. C C
  45. C IPRC (LE PRECONDITIONNEMENT DANS UN LISTREEL) C
  46. C C
  47. C DCG : IPRC EST LE POINTEUR SUR UN LISTREEL C
  48. C CONTENANT D**(-1/2) CE LISTREEL A POUR C
  49. C DIMENSION NL AU MOINS C
  50. C C
  51. C IA (TABLEAU DE POINTEURS DU STOCKAGE MORSE) C
  52. C C
  53. C IL DOIT ETRE DIMENSIONNE A NL+1 AU MINIMUM C
  54. C C
  55. C JA (TABLEAU DE CONNECTIVITES DE LA MATRICE) C
  56. C C
  57. C IL EST STOCKE EN MORSE ET DE LONGUEUR LA C
  58. C C
  59. C KA (TABLEAU DE CONNECTIVITES EN MODE COMPRESSE) C
  60. C C
  61. C IL EST DE LA FORME KA(NL,NNZ) AVEC NNZ C
  62. C LE NOMBRE MAXI DE TERMES NON NULS POUR UNE C
  63. C LIGNE DE LA MATRICE C
  64. C C
  65. C NB : EN MORSE, KA N'EST PAS UTILISE C
  66. C EN COMPRESSE, IA ET JA NE SONT PAS UTILISES C
  67. C C
  68. C IA,JA ET KA SONT DES LISTENTI. C
  69. C C
  70. C***** INFOS SUR LE SECOND MEMBRE, LA SOLUTION INITIALE C
  71. C***** ET LA SOLUTION FINALE. (LIGNE 3) C
  72. C C
  73. C IZB SECOND MEMBRE EN ENTREE C
  74. C PAS MODIFIE. C
  75. C C
  76. C IZP ESTIMATION DE LA SOLUTION EN ENTREE C
  77. C LA SOLUTION EN SORTIE, IL FAUT LA C
  78. C CONSERVER POUR LA RESOLUTION SUIVANTE C
  79. C C
  80. C IZB ET IZP SONT DES LISTREEL DIMENSIONNES A NL C
  81. C C
  82. C ICOLD =1 (DEMARRAGE FROID) C
  83. C ICOLD!=1 (DEMARRAGE CHAUD A PARTIR DE LA C
  84. C SOLUTION INITIALE FOURNIE DANS IZP) C
  85. C C
  86. C***** INFOS SUR LA CONVERGENCE (LIGNE 4) C
  87. C C
  88. C NPT NOMBRE MAXI D'ITERATIONS EN ENTREE C
  89. C SI IL EST < 10, IL EST MIS A 10 C
  90. C NPITE FREQUENCE DES TESTS DE CONVERGENCE C
  91. C NEFF NOMBRE REEL EN SORTIE C
  92. C ICONV 0 NON CONVERGE, 1 CONVERGE C
  93. C EPI PRECISION DEMANDEE C
  94. C IPOU=1 LA PRECISION VA ETRE POUSSE D'UN FACTEUR C
  95. C VPOU. CECI PEUT S'AVERER UTILE LORS D'UN C
  96. C CALCUL ITERATIF OU ON VEUT DES LE DEBUT C
  97. C UNE SOLUTION PRECISE POUR EVITER DE PARTIR C
  98. C SUR UNE MAUVAISE PISTE. VPOU A UNE VALEUR C
  99. C PAR DEFAUT DE 100. C
  100. C C
  101. C***** INFOS SUR LE CONTROLE DES IMPRESSIONS (LIGNE 5) C
  102. C C
  103. C NIMPR =1 (INFOS SUR LA CONVERGENCE) C
  104. C NIMPR!=1 (SILENCIEUX SAUF EN CAS DE DIVERGENCE) C
  105. C C
  106. C IPAT = (FACULTATIF, PAR EXEMPLE UN PAS DE TEMPS)C
  107. C C
  108. C RSETD(B,A) B <- A C
  109. C ADIVEC(A,B,C,alpha) C <- A + alpha*B C
  110. C PIMDV(A,B,C) C <- A*B C
  111. C C
  112. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  113. C
  114. IMPLICIT INTEGER(I-N)
  115. IMPLICIT REAL*8 (A-H,O-Z)
  116. -INC SMLENTI
  117. POINTEUR IA.MLENTI
  118. POINTEUR JA.MLENTI
  119. POINTEUR KA.MLENTI
  120. -INC SMLREEL
  121. POINTEUR IMAT.MLREEL
  122. POINTEUR IPRC.MLREEL
  123. POINTEUR IZB.MLREEL
  124. POINTEUR IZBB.MLREEL
  125. POINTEUR IZP.MLREEL
  126. POINTEUR IZPS.MLREEL
  127. POINTEUR IZZPR.MLREEL
  128. POINTEUR IZPMV.MLREEL
  129. POINTEUR IZTMP.MLREEL
  130. POINTEUR IZDIR.MLREEL
  131.  
  132. DATA VPOU0/100./
  133.  
  134. IF(KSTO.EQ.0) LA=ILG
  135. IF(KSTO.EQ.1) NNZ=ILG
  136. IF(IPOU.EQ.1) THEN
  137. VPOUE=VPOU0
  138. IF(VPOU.NE.VPOU0) VPOUE=VPOU
  139. ENDIF
  140.  
  141. CALL OOOMRU(1)
  142.  
  143. C On sauve le second membre en le dupliquant dans IZBB
  144. C qui lui sera ecrase.
  145. SEGACT IZB
  146. JG=IZB.PROG(/1)
  147. SEGINI IZBB
  148. CALL RSETD(IZBB.PROG,IZB.PROG,JG)
  149. SEGDES IZB
  150.  
  151. SEGACT IZP
  152.  
  153. JG=NL
  154. SEGINI IZPS
  155. SEGINI IZZPR
  156. SEGINI IZPMV
  157. SEGINI IZTMP
  158. SEGINI IZDIR
  159.  
  160. SEGACT IMAT
  161. SEGACT IPRC
  162. IF(KSTO.EQ.0) SEGACT IA
  163. IF(KSTO.EQ.0) SEGACT JA
  164. IF(KSTO.EQ.1) SEGACT KA
  165.  
  166. IF(NPT.LT.10) NPT=10
  167. IF(NPT.GT.NL) NPT=NL
  168.  
  169. C Mettre le second membre et l'estimation initiale
  170. C a l'echelle, B <- D-(1/2)*B, PIM <- PIM/(D**(-1/2))
  171. CALL PIMDV(IZBB.PROG,IPRC.PROG,IZBB.PROG,NL)
  172. CALL DIVISE(NL,IZP.PROG,IZP.PROG,IPRC.PROG)
  173. IF(ICOLD.EQ.1) CALL INITD (IZP.PROG,NL,0.D0)
  174. CALL RSETD(IZPS.PROG,IZP.PROG,NL)
  175. CALL PIMDV(IZP.PROG,IPRC.PROG,IZTMP.PROG,NL)
  176. IF(KSTO.EQ.0) THEN
  177. C CALL PMVM(A,IA,JA,X,Y,NL,LA)
  178. CALL PMVM(IMAT.PROG,IA.LECT,JA.LECT,
  179. 1 IZTMP.PROG,IZPMV.PROG,NL,LA)
  180. ELSE
  181. C CALL PMVC(A,KA,X,Y,NL,NNZ)
  182. CALL PMVC(IMAT.PROG,KA.LECT,
  183. 1 IZTMP.PROG,IZPMV.PROG,NL,NNZ)
  184. ENDIF
  185. CALL PIMDV(IZPMV.PROG,IPRC.PROG,IZPMV.PROG,NL)
  186. CALL ADIVEC (IZBB.PROG,IZPMV.PROG,IZBB.PROG,-1.d0,NL)
  187. ANORB=SQRT(PISCAL(IZBB.PROG,IZBB.PROG,NL))
  188. CALL RSETD(IZZPR.PROG,IZBB.PROG,NL)
  189. CALL RSETD(IZDIR.PROG,IZZPR.PROG,NL)
  190. SCAR1=PISCAL(IZBB.PROG,IZDIR.PROG,NL)
  191. C
  192. C===========DEBUT DE LA BOUCLE ITERATIVE=====600========
  193. C
  194. NEFF=0
  195. NCOMP=0
  196. ICONV=0
  197. DO 600 I=1,NPT
  198. NEFF=NEFF+1
  199. NCOMP=NCOMP+1
  200. CALL PIMDV(IZDIR.PROG,IPRC.PROG,IZTMP.PROG,NL)
  201. IF(KSTO.EQ.0) THEN
  202. C CALL PMVM(A,IA,JA,X,Y,NL,LA)
  203. CALL PMVM(IMAT.PROG,IA.LECT,JA.LECT,
  204. 1 IZTMP.PROG,IZPMV.PROG,NL,LA)
  205. ELSE
  206. C CALL PMVC(A,KA,X,Y,NL,NNZ)
  207. CALL PMVC(IMAT.PROG,KA.LECT,
  208. 1 IZTMP.PROG,IZPMV.PROG,NL,NNZ)
  209. ENDIF
  210. CALL PIMDV(IZPMV.PROG,IPRC.PROG,IZPMV.PROG,NL)
  211. SCADIR=PISCAL(IZDIR.PROG,IZPMV.PROG,NL)
  212. ALPHA=SCAR1/SCADIR
  213. CALL ADIVEC(IZP.PROG,IZDIR.PROG,IZP.PROG,ALPHA,NL)
  214. ALP=-1.D0*ALPHA
  215. CALL ADIVEC (IZBB.PROG,IZPMV.PROG,IZBB.PROG,ALP,NL)
  216. CALL RSETD(IZZPR.PROG,IZBB.PROG,NL)
  217. SCAR2=PISCAL(IZBB.PROG,IZZPR.PROG,NL)
  218. BETA=SCAR2/SCAR1
  219. CALL ADIVEC (IZZPR.PROG,IZDIR.PROG,IZDIR.PROG,BETA,NL)
  220. IF (NCOMP.EQ.NPITE.OR.I.EQ.NPT) THEN
  221. IF (I.EQ.NPT) INDIC=1
  222. C On pousse la precision si necessaire
  223. ZEPS=EPI
  224. IF(IPOU.EQ.1) ZEPS=EPI/VPOUE
  225. CALL COTEST (IZP.PROG,IZPS.PROG,
  226. 1 ZEPS,DELTA,DELTAP,
  227. 2 INDIC,ICONV,NL)
  228. SEPA=DELTA
  229. RSEPA=DELTA/DELTAP
  230. CALL RSETD(IZPS.PROG,IZP.PROG,NL)
  231. NCOMP=0
  232. ENDIF
  233. IF (ICONV.EQ.1) GOTO 601
  234. SCAR1=SCAR2
  235. 600 CONTINUE
  236. 601 CONTINUE
  237. CALL MTABL(IZBB.PROG,RES,IMA,NL)
  238. RRES=RES/ANORB
  239. C Pour le DCG la solution est D**(-1/2)*IZP
  240. C Pour le DCG on remet dans PIM, la solution
  241. CALL PIMDV(IZP.PROG,IPRC.PROG,IZP.PROG,NL)
  242. IF (ICONV.NE.1) THEN
  243. WRITE(6,1111)IPAT,NPT,RSEPA,DELTAP,RES,IZBB.PROG(IMA)
  244. ELSE
  245. IF(NIMPR.EQ.1) WRITE(6,2222)IPAT,EPI,
  246. 1 RES,SEPA,I,
  247. 2 RRES,RSEPA,IMA
  248. ENDIF
  249.  
  250. SEGDES IZP
  251.  
  252. SEGSUP IZBB
  253. SEGSUP IZPS
  254. SEGSUP IZZPR
  255. SEGSUP IZDIR
  256. SEGSUP IZTMP
  257. SEGSUP IZPMV
  258.  
  259. SEGDES IMAT*(NOMOD,MRU)
  260. SEGDES IPRC*(NOMOD,MRU)
  261. IF(KSTO.EQ.0) SEGDES IA*(NOMOD,MRU)
  262. IF(KSTO.EQ.0) SEGDES JA*(NOMOD,MRU)
  263. IF(KSTO.EQ.1) SEGDES KA*(NOMOD,MRU)
  264.  
  265. CALL OOOMRU(0)
  266.  
  267. C
  268. 1111 FORMAT(1X,'*CGRESD* ',I6,' NON CONV EN ',I5,
  269. 1' ITERA, RSEPA = ',D8.2,1X,'DELTAP = ',D8.2,1X,
  270. 2' RES = ',D8.2,' P = ',D8.2)
  271. 2222 FORMAT(1X,'*CGRESD* ',I6,' *PREC ',D8.2,
  272. 1 ' * RESIDU-SEPARA @ NITE ',D8.2,1X,D8.2,1X,I5,
  273. 2 ' * RRES-RSEP ',D8.2,1X,D8.2,
  274. 3 ' * ELEMENT ',I7)
  275. C
  276. C FIN RESOLUTION
  277. C
  278. RETURN
  279. END
  280.  
  281.  
  282.  

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