Télécharger graco6.eso

Retour à la liste

Numérotation des lignes :

  1. C GRACO6 SOURCE FANDEUR 16/11/28 21:15:06 9208
  2.  
  3. SUBROUTINE GRACO6 ( ICHOLX, MSECO ,NOEN,MSOL,lenb)
  4. *
  5. * EXECUTE LES ITERATIONS DU GRADIENTS CONJUGUE.
  6. *
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8(A-H,O-Z)
  9.  
  10. -INC CCOPTIO
  11. -INC CCASSIS
  12.  
  13. -INC SMMATRI
  14. -INC SMVECTD
  15.  
  16. segment ilicre
  17. * stockage matrice initiale en creux
  18. * ilideb position debut de ligne dans ligcre
  19. integer ilideb(nbinc+1)
  20. integer iliinc(ino+1)
  21. integer ligcrp
  22. endsegment
  23. segment ligcre
  24. * lmatr: longueur reelle ligne
  25. * posm: numero inconnue
  26. * valm: valeur terme
  27. integer posm(lmat)
  28. real*8 valm(lmat)
  29. endsegment
  30. POINTEUR MVECT4.MVECTD,MVECT5.MVECTD,MVECT6.MVECTD,
  31. $ MVECT7.MVECTD
  32.  
  33. nbthr=nbthrs
  34. nbthr=min(nbthr,64)
  35. if (nbthr.gt.1) call threadii
  36.  
  37. MMATRI=ICHOLX
  38. * activation de la matrice une fois pour toute.
  39. SEGACT,MMATRI
  40. MILIGN=IASLIG
  41. SEGACT,MILIGN
  42. INO=ILIGN(/1)
  43. MDNOR=IDNORM
  44. SEGACT MDNOR
  45. DO I=1,INO
  46. LLIGN=ILIGN(I)
  47. SEGACT LLIGN
  48. enddo
  49. MVECTD = MSECO
  50. SEGACT MVECTD
  51. INC = VECTBB(/1)
  52. MAXIT=INC-1
  53. MILIGN=IILIGN
  54. SEGACT MILIGN
  55. SEGINI MVECT1, MVECT3,MVECT4
  56. segini mvect2,mvect6
  57. * conversion matrice en structure ligne complete creuse
  58. call graco11(mmatri,ilicre)
  59. * idem matrice factorisee
  60. call graco12(mmatri,ifacre,ifatra)
  61. *
  62. * on fait une premiere pseudo resolution pour avoir un ordre de
  63. * grandeur de l'energie
  64. *
  65. MVECTD=MSECO
  66. DO 31 IB=1,INC
  67. MVECT1.VECTBB(IB)=VECTBB(IB)*DNOR(IB)
  68. 31 CONTINUE
  69. CALL GRACO8 ( ICHOLX,MVECT1,NOEN,ifacre,ifatra)
  70. SEGACT MILIGN
  71. *p DO 38 IB=1,INC
  72. *p 38 CONTINUE
  73. CALL GRACO7(ilicre,MVECT1,MVECT2,inc,nbthr,lenb)
  74. ENEI = DDOT2(INC,MVECT1.VECTBB(1),MVECT2.VECTBB(1))
  75. IF (IIMPI.NE.0) THEN
  76. WRITE(IOIMP,*) ' premiere estimation de l energie',ENEI
  77. ENDIF
  78. IF (ENEI.LE.0.D0) THEN
  79. write(ioimp,*) 'ATTENTION : ENEI = ',ENEI,' !'
  80. * INTERR(1) = 0
  81. * CALL ERREUR(49)
  82. * RETURN
  83. ENDIF
  84. MVECT7=MSECO
  85. DO 3 IB=1,INC
  86. MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT2.VECTBB(IB)
  87. 3 CONTINUE
  88. ** SEGDES MVECT7
  89. DO 59 IB=1,INC
  90. MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB)
  91. 59 CONTINUE
  92. CALL GRACO8(ICHOLX,MVECT3,NOEN,ifacre,ifatra)
  93. DO 60 IB=1,INC
  94. MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB)
  95. 60 CONTINUE
  96. C
  97. C DEBUT DE TOURNER EN ROND : MAXIT = INC-1
  98. C
  99. IF (IIMPI.NE.0) THEN
  100. WRITE(ioimp,*)' Debut des iterations (max =',MAXIT,')'
  101. ENDIF
  102. C
  103. DO 100 iter = 1, MAXIT
  104.  
  105. IF (IERR.NE.0) RETURN
  106. *
  107. * Recalcul du critere et du residu toutes les 100 iterations :
  108. IF (iter/100 * 100 - iter .EQ. 0) then
  109. SEGACT MILIGN
  110. CALL GRACO7( ilicre,MVECT1,MVECT6,inc,nbthr,lenb)
  111. ENEN = DDOT2(INC,MVECT1.VECTBB(1),MVECT6.VECTBB(1))
  112. IF (IIMPI.NE.0) THEN
  113. WRITE(ioimp,*)' Nouvelle energie ',ENEN,crit,'(iter=',iter,')'
  114. ENDIF
  115. IF (ENEN.LE.0.D0) THEN
  116. write(ioimp,*) 'ATTENTION : ENEN = ',ENEN,' !'
  117. * INTERR(1) = 0
  118. * CALL ERREUR(49)
  119. * RETURN
  120. ENDIF
  121. ENEI = ENEN
  122. IF (IIMPI.NE.0) THEN
  123. if (iter/10 * 10 - iter .eq. 0) then
  124. WRITE(IOIMP,FMT='('' ITERATION '',I6,'' CRITERE '',E12.5)')
  125. * ITER ,CRIT
  126. endif
  127. ENDIF
  128. * recalcul residu reel
  129. * IF (iter/1000 * 1000 - iter .EQ. 0) then
  130. * DO IB=1,INC
  131. * MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT6.VECTBB(IB)
  132. * ENDDO
  133. * ENDIF
  134. ENDIF
  135. * fin recalcul
  136. CALL GRACO7(Ilicre,MVECT4,MVECT6,inc,nbthr,lenb)
  137. if (iter.EQ.1) then
  138. R0RR0 = DDOT2(INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1))
  139. endif
  140. P0AP0 = DDOT2(INC,MVECT4.VECTBB(1),MVECT6.VECTBB(1))
  141. IF (P0AP0.EQ.0.D0) P0AP0=1d-50
  142. IF (ABS(P0AP0).lt.1d-45) write(ioimp,*) ' p0ap0 ',p0ap0
  143. ALP0 = R0RR0 / P0AP0
  144. DO 6 IB = 1, INC
  145. MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)+ALP0*MVECT4.VECTBB(IB)
  146. 6 CONTINUE
  147. DO 7 IB=1,INC
  148. MVECT2.VECTBB(IB)=MVECT2.VECTBB(IB)-ALP0*MVECT6.VECTBB(IB)
  149. 7 CONTINUE
  150. DO 61 IB=1,INC
  151. MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB)
  152. 61 CONTINUE
  153. CALL GRACO8(ICHOLX,MVECT3,NOEN,ifacre,ifatra)
  154. R1RR1 = DDOT2(INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1))
  155. CRIT = R1RR1 / ENEI
  156. IF (IIMPI.GT.1) WRITE(ioimp,*) ' critere ' , crit,iter
  157. IF (ABS(CRIT) .LT. 1.d-20 .OR. ITER.GT.80000) THEN
  158. * on a converge (mais en fait pas toujours si ITER est trop grand !)
  159. GO TO 101
  160. ENDIF
  161. IF (R0RR0.EQ.0.D0) R0RR0=1d-30
  162. IF (ABS(R0RR0).lt.1d-45) write(ioimp,*) ' r0rr0 ',r0rr0
  163. *pv IF (R0RR0.EQ.0.D0) THEN
  164. *pv CALL ERREUR(835)
  165. *pVC WRITE(IOIMP,FMT='('' PROBLEME DE DIVISION PAR ZERO '')')
  166. *pv RETURN
  167. *pv ENDIF
  168. BET1 = R1RR1 / R0RR0
  169. DO 9 IB=1, INC
  170. MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB)+BET1*MVECT4.VECTBB(IB)
  171. 9 CONTINUE
  172. R0RR0 = R1RR1
  173. 100 CONTINUE
  174. * Fin des iterations !
  175. if (nbthr.gt.1) call threadis
  176. CALL ERREUR(460)
  177. C WRITE( IOIMP,FMT='('' PAS DE CONVERGENCE'')')
  178. RETURN
  179. 101 MSOL = MVECT1
  180. if (nbthr.gt.1) call threadis
  181. if (iter.gt.80000) then
  182. write(ioimp,*) 'pas de convergence en ',iter,' iterations',crit
  183. else
  184. write(ioimp,*) ' convergence en ',iter, ' iterations',crit
  185. endif
  186. DO 67 IB=1,INC
  187. MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)*DNOR(IB)
  188. 67 CONTINUE
  189. SEGDES MDNOR,MMATRI,MILIGN
  190. SEGSUP MVECT2,MVECT3,MVECT4,MVECT6
  191. ligcre=ligcrp
  192. segsup ilicre,ligcre
  193. ilicre=ifacre
  194. ligcre=ligcrp
  195. segsup ilicre,ligcre
  196. ilicre=ifatra
  197. ligcre=ligcrp
  198. segsup ilicre,ligcre
  199.  
  200. RETURN
  201. END
  202.  
  203.  
  204.  

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