Télécharger la plaquette Cast3M


fr en

Télécharger graco6.eso

Retour à la liste

Numérotation des lignes : 

graco6

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

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