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    PV        22/02/24    21:15:03     11299          
  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.       MAXIT=INC-1
  55.       MILIGN=IILIGN
  56.       SEGACT MILIGN
  57.       SEGINI MVECT1, MVECT3,MVECT4
  58.       segini mvect2,mvect6
  59. *  conversion matrice en structure ligne complete creuse
  60.       call graco11(mmatri,ilicre,0)
  61. *  idem matrice factorisee
  62.       call graco12(mmatri,ifacre,ifatra)
  63. *
  64. *  on fait une premiere pseudo resolution pour avoir un ordre de
  65. *  grandeur de l'energie
  66. *
  67.       MVECTD=MSECO
  68.       DO 31 IB=1,INC
  69.         MVECT1.VECTBB(IB)=VECTBB(IB)*DNOR(IB)
  70.   31  CONTINUE
  71.       CALL GRACO8 ( ICHOLX,MVECT1,NOEN,ifacre,ifatra)
  72.       SEGACT MILIGN
  73. *p    DO 38 IB=1,INC
  74. *p 38 CONTINUE
  75.       CALL GRACO7(ilicre,MVECT1,MVECT2,inc,nbthr,lenb)
  76.       ENEI = DDOT2(INC,MVECT1.VECTBB(1),MVECT2.VECTBB(1))
  77.       IF (IIMPI.NE.0) THEN
  78.         WRITE(IOIMP,*) ' premiere estimation de l energie',ENEI
  79.       ENDIF
  80.       IF (ENEI.LE.0.D0) THEN
  81.         write(ioimp,*) 'ATTENTION : ENEI = ',ENEI,' !'
  82. *       INTERR(1) = 0
  83. *       CALL ERREUR(49)
  84. *      if (nbthr.gt.1) call threadis
  85. *       RETURN
  86.       ENDIF
  87.       MVECT7=MSECO
  88.       DO 3 IB=1,INC
  89.         MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT2.VECTBB(IB)
  90.   3   CONTINUE
  91. **    SEGDES MVECT7
  92.       DO 59 IB=1,INC
  93.         MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB)
  94.   59  CONTINUE
  95.       CALL  GRACO8(ICHOLX,MVECT3,NOEN,ifacre,ifatra)
  96.       DO 60 IB=1,INC
  97.         MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB)
  98.   60  CONTINUE
  99. C
  100. C  DEBUT DE TOURNER EN ROND : MAXIT = INC-1
  101. C
  102.       IF (IIMPI.NE.0) THEN
  103.         WRITE(ioimp,*)' Debut des iterations (max =',MAXIT,')'
  104.       ENDIF
  105. C
  106.       DO 100 iter = 1, MAXIT
  107.  
  108.         IF (IERR.NE.0) then
  109.          if (nbthr.gt.1) call threadis
  110.          RETURN
  111.         endif
  112. *
  113. * Recalcul du critere et du residu toutes les 100 iterations :
  114.       IF (iter/100 * 100  - iter .EQ. 0) then
  115.        SEGACT MILIGN
  116.        CALL GRACO7( ilicre,MVECT1,MVECT6,inc,nbthr,lenb)
  117.        ENEN = DDOT2(INC,MVECT1.VECTBB(1),MVECT6.VECTBB(1))
  118.        IF (IIMPI.NE.0) THEN
  119.          WRITE(ioimp,*)' Nouvelle energie ',ENEN,crit,'(iter=',iter,')'
  120.        ENDIF
  121.        IF (ENEN.LE.0.D0) THEN
  122.          write(ioimp,*) 'ATTENTION : ENEN = ',ENEN,' !'
  123. *        INTERR(1) = 0
  124. *        CALL ERREUR(49)
  125. *        if (nbthr.gt.1) call threadis
  126. *        RETURN
  127.        ENDIF
  128.        ENEI = ENEN
  129.        IF (IIMPI.NE.0) THEN
  130.         if (iter/10 * 10  - iter .eq. 0) then
  131.          WRITE(IOIMP,FMT='('' ITERATION '',I6,'' CRITERE '',E12.5)')
  132.      *   ITER ,CRIT
  133.         endif
  134.        ENDIF
  135. *  recalcul residu reel
  136. *       IF (iter/1000 * 1000  - iter .EQ. 0) then
  137. *        DO   IB=1,INC
  138. *         MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT6.VECTBB(IB)
  139. *        ENDDO
  140. *       ENDIF
  141.       ENDIF
  142. *  fin recalcul
  143.       CALL GRACO7(Ilicre,MVECT4,MVECT6,inc,nbthr,lenb)
  144.       if (iter.EQ.1) then
  145.        R0RR0 = DDOT2(INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1))
  146.       endif
  147.       P0AP0 = DDOT2(INC,MVECT4.VECTBB(1),MVECT6.VECTBB(1))
  148.       IF (P0AP0.EQ.0.D0) P0AP0=1d-50
  149.       IF (ABS(P0AP0).lt.1d-45) write(ioimp,*) ' p0ap0 ',p0ap0
  150.       ALP0 = R0RR0 / P0AP0
  151.       DO 6 IB = 1, INC
  152.         MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)+ALP0*MVECT4.VECTBB(IB)
  153.   6   CONTINUE
  154.       DO 7 IB=1,INC
  155.         MVECT2.VECTBB(IB)=MVECT2.VECTBB(IB)-ALP0*MVECT6.VECTBB(IB)
  156.   7   CONTINUE
  157.       DO 61 IB=1,INC
  158.  
  159.         MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB)
  160.   61  CONTINUE
  161.       CALL GRACO8(ICHOLX,MVECT3,NOEN,ifacre,ifatra)
  162.       R1RR1 = DDOT2(INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1))
  163.       CRIT = R1RR1 / ENEI
  164.       IF (IIMPI.GT.1) WRITE(ioimp,*) ' critere ' , crit,iter
  165.       IF (ABS(CRIT) .LT. 1.d-20 .OR. ITER.GT.80000)  THEN
  166. * on a converge (mais en fait pas toujours si ITER est trop grand !)
  167.         GO TO 101
  168.       ENDIF
  169.       IF (R0RR0.EQ.0.D0) R0RR0=1d-30
  170.       IF (ABS(R0RR0).lt.1d-45) write(ioimp,*) ' r0rr0 ',r0rr0
  171. *pv   IF (R0RR0.EQ.0.D0) THEN
  172. *pv     CALL ERREUR(835)
  173. *pVC     WRITE(IOIMP,FMT='('' PROBLEME DE DIVISION PAR ZERO '')')
  174. *pv     RETURN
  175. *pv   ENDIF
  176.       BET1 = R1RR1 / R0RR0
  177.       DO 9 IB=1, INC
  178.         MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB)+BET1*MVECT4.VECTBB(IB)
  179.    9  CONTINUE
  180.       R0RR0 = R1RR1
  181.  100  CONTINUE
  182. * Fin des iterations !
  183.       if (nbthr.gt.1) call threadis
  184.       CALL ERREUR(460)
  185. C      WRITE( IOIMP,FMT='(''   PAS DE CONVERGENCE'')')
  186.       RETURN
  187.  101  MSOL = MVECT1
  188.       if (nbthr.gt.1) call threadis
  189.       if (iter.gt.80000) then
  190.         write(ioimp,*) 'pas de convergence en ',iter,' iterations',crit
  191.       else
  192.         write(ioimp,*) ' convergence en ',iter, ' iterations',crit
  193.       endif
  194.       DO 67 IB=1,INC
  195.         MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)*DNOR(IB)
  196.   67  CONTINUE
  197.       SEGSUP MVECT2,MVECT3,MVECT4,MVECT6
  198. *  suppression stockage morse matrice et transposee
  199.       SEGACT,MMATRI*mod
  200.       ilicre=jlicre
  201.       ligcre=ligcrp
  202.       segsup ilicre,ligcre
  203.       jlicre=0
  204.       ilicre=ifacre
  205.       ligcre=ligcrp
  206.       segsup ilicre,ligcre
  207.       ifacre=0
  208.       ilicre=ifatra
  209.       ligcre=ligcrp
  210.       ifatra=0
  211.       segsup ilicre,ligcre
  212.       SEGDES MDNOR,MMATRI,MILIGN
  213.  
  214.       RETURN
  215.       END
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  
  222.  
  223.  

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