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

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