Télécharger la plaquette Cast3M


fr en

Télécharger bcgstb.eso

Retour à la liste

Numérotation des lignes : 

bcgstb

  1. C BCGSTB    SOURCE    MB234859  26/09/22    21:15:05     12653          
  2.       SUBROUTINE BCGSTB(N,A,B,X,XCRIT,MAXIT,KERRE)
  3. C----------------------------------------------------------------------
  4. C            RESOLUTION DE AX=B PAR L'ALGORITHME BiCGSTAB
  5. C
  6. C Reference : "Templates for the Solution of Linear Systems :
  7. C              Building Blocks for Iterative Methods", R.Barrett et al.
  8. C
  9. C     Entrees :
  10. C     ---------
  11. C     N     : ENTIER donnant la taille de la matrice et des vecteurs
  12. C     A     : MLENTI stockant la matrice par lignes.
  13. C             Chaque ligne i a ses valeurs contenues dans un MLREEL
  14. C             dont le pointeur est stocke dans A.LECT(i).
  15. C     B     : MLREEL donnant le vecteur second membre
  16. C     MAXIT : ENTIER donnant le nombre d'iteration max souhaite
  17. C     XPREC : REEL   donnant la precision souhaitee
  18. C
  19. C     Sorties :
  20. C     ---------
  21. C     X     : MLREEL donnant le vecteur solution
  22. C     KERRE : ENTIER precisant si tout s'est bien passe (0) ou non.
  23. C----------------------------------------------------------------------
  24. C     DECLARATIONS
  25. C----------------------------------------------------------------------
  26.       IMPLICIT INTEGER(I-N)
  27.       IMPLICIT REAL*8 (A-H,O-Z)
  28. C
  29. -INC PPARAM
  30. -INC CCASSIS
  31. -INC CCOPTIO
  32. -INC CCREEL
  33. -INC SMLENTI
  34.       POINTEUR A.MLENTI
  35. -INC SMLREEL
  36.       POINTEUR X.MLREEL,B.MLREEL
  37.       POINTEUR R.MLREEL,R0.MLREEL,P.MLREEL,S.MLREEL,T.MLREEL,V.MLREEL
  38. C
  39.       REAL*8 ALPHA,BETA,OMEGA,RHOOLD,RHONEW,BNORM,SNORM,TNORM,VPSCA
  40. C
  41. C     Informations pour la parallelisation
  42.       LOGICAL BTHRD,BSGDES
  43.       COMMON/PDTSCA/NBTHRD,NVAL,NLIG,IMAT,IVEC,IRES
  44.       EXTERNAL BCGSi
  45. C----------------------------------------------------------------------
  46. C     INITIALISATIONS
  47. C----------------------------------------------------------------------
  48.       KERRE  = 0
  49.       ITER   = 0
  50.       ALPHA  = 0.0D0
  51.       BETA   = 0.0D0
  52.       OMEGA  = 1.0D0
  53.       RHOOLD = 1.0D0
  54. C
  55.       JG=N
  56.       SEGINI,R,R0,P,S,T,V
  57. C
  58.       NLA=0
  59.       BSGDES=.FALSE.
  60.       DO I=1,N
  61.         MLREE1 = A.LECT(I)
  62.         SEGACT /ERR=3/ MLREE1
  63.         NLA=NLA+1
  64.       ENDDO
  65.       GOTO 4
  66. C
  67.  3    CONTINUE
  68.       BSGDES=.TRUE.
  69.       MLREE1 = A.LECT(NLA)
  70.       SEGDES,MLREE1
  71.       NLA=NLA-1
  72.  4    CONTINUE
  73. C
  74. C     Norme du second membre b
  75.       BNORM = DDOTPV(N,B.PROG(1),B.PROG(1))
  76.       BNORM = SQRT(BNORM)
  77.       IF (BNORM.LE.XPETIT) BNORM = 1.0D0
  78. C
  79. C     Residu initial R0 = b - Ax0
  80.       DO I=1,NLA
  81.         MLREE1    = A.LECT(I)
  82.         R.PROG(I) = B.PROG(I) - DDOTPV(N,MLREE1.PROG(1),X.PROG(1))
  83.       ENDDO
  84.       IF (BSGDES) THEN
  85.         DO I=NLA+1,N
  86.           MLREE1    = A.LECT(I)
  87.           SEGACT,MLREE1
  88.           R.PROG(I) = B.PROG(I) - DDOTPV(N,MLREE1.PROG(1),X.PROG(1))
  89.           SEGDES,MLREE1
  90.         ENDDO
  91.       ENDIF
  92. C
  93. C     Convergence atteinte?
  94.       SNORM = DDOTPV(N,R.PROG(1),R.PROG(1))
  95.       IF (SQRT(SNORM)/BNORM .LT. XCRIT) GOTO 200
  96. C
  97.       DO I=1,N
  98.         R0.PROG(I) = R.PROG(I)
  99.         P .PROG(I) = R.PROG(I)
  100.       ENDDO
  101. C
  102. C     Parallelisation des produits matrice-vecteurs
  103.       NBTHR = MIN(MAX(N*NLA/15000000,1),NBTHRS)
  104.       BTHRD = .TRUE.
  105.       IF ((NBTHRS.EQ.1).OR.(OOTHRD.GT.0)) THEN
  106.         NBTHR = 1
  107.         BTHRD = .FALSE.
  108.       ENDIF
  109.       IF (BTHRD) CALL THREADII
  110.       NBTHRD=NBTHR
  111.       NVAL=N
  112.       NLIG=NLA
  113.       IMAT=A
  114. C----------------------------------------------------------------------
  115. C     BOUCLE PRINCIPALE DES ITERATIONS
  116. C----------------------------------------------------------------------
  117.  10   CONTINUE
  118.       ITER = ITER + 1
  119.       IF (IERR.NE.0) GOTO 100
  120. C
  121. C     Mise a jour de RHONEW = <R0, R(iter)>
  122.       RHONEW = DDOTPV(N,R0.PROG(1),R.PROG(1))
  123.       IF (ABS(RHONEW).LT.XPETIT) THEN
  124.         KERRE = 1
  125.         GOTO 100
  126.       ENDIF
  127. C
  128. C     Coefficient de direction + maj de la direction de recherche
  129.       IF (ITER.GT.1) THEN
  130.         BETA = (RHONEW / RHOOLD) * (ALPHA / OMEGA)
  131.         DO I=1,N
  132.           P.PROG(I) = R.PROG(I) + BETA * (P.PROG(I) - OMEGA * V.PROG(I))
  133.         ENDDO
  134.       ENDIF
  135. C
  136. C     Calcul de V=AP
  137.       IVEC=P
  138.       IRES=V
  139.       DO ITH=2,NBTHR
  140.         CALL THREADID(ITH,BCGSI)
  141.       ENDDO
  142.       CALL BCGSI(1)
  143.       DO ITH=2,NBTHR
  144.         CALL THREADIF(ITH)
  145.       ENDDO
  146.       IF (BSGDES) THEN
  147.         DO I=NLA+1,N
  148.           MLREE1    = A.LECT(I)
  149.           SEGACT,MLREE1
  150.           V.PROG(I) = DDOTPV(N,MLREE1.PROG(1),X.PROG(1))
  151.           SEGDES,MLREE1
  152.         ENDDO
  153.       ENDIF
  154. C
  155. C     Calcul du pas ALPHA = RHONEW / VPSCA  avec VPSCA = <R0, V>
  156.       VPSCA = DDOTPV(N,R0.PROG(1),V.PROG(1))
  157.       IF (ABS(VPSCA).LT.XPETIT) THEN
  158.         KERRE = 2
  159.         GOTO 100
  160.       ENDIF
  161.       ALPHA = RHONEW / VPSCA
  162. C
  163. C     Calcul du résidu intermédiaire S = R - ALPHA * V + test de convergence
  164.       DO I=1,N
  165.         S.PROG(I) = R.PROG(I) - ALPHA * V.PROG(I)
  166.       ENDDO
  167. C
  168. C     Premier test de convergence : ||S|| / ||B|| < XCRIT
  169.       SNORM = DDOTPV(N,S.PROG(1),S.PROG(1))
  170.       IF (SQRT(SNORM) / BNORM .LT. XCRIT) THEN
  171.         DO I=1,N
  172.           X.PROG(I) = X.PROG(I) + ALPHA * P.PROG(I)
  173.         ENDDO
  174.         GOTO 100
  175.       ENDIF
  176. C
  177. C     Calcul du second pas T = AS
  178.       IVEC=S
  179.       IRES=T
  180.       DO ITH=2,NBTHR
  181.         CALL THREADID(ITH,BCGSI)
  182.       ENDDO
  183.       CALL BCGSI(1)
  184.       DO ITH=2,NBTHR
  185.         CALL THREADIF(ITH)
  186.       ENDDO
  187.       IF (BSGDES) THEN
  188.         DO I=NLA+1,N
  189.           MLREE1    = A.LECT(I)
  190.           SEGACT,MLREE1
  191.           V.PROG(I) = DDOTPV(N,MLREE1.PROG(1),X.PROG(1))
  192.           SEGDES,MLREE1
  193.         ENDDO
  194.       ENDIF
  195. C
  196. C     Parametre de stabilisation OMEGA = <S, T> / <T, T>
  197.       TNORM = DDOTPV(N,T.PROG(1),T.PROG(1))
  198.       VPSCA = DDOTPV(N,S.PROG(1),T.PROG(1))
  199.       IF (TNORM.EQ.0.0D0) THEN
  200.         KERRE = 3
  201.         GOTO 100
  202.       ENDIF
  203.       OMEGA = VPSCA / TNORM
  204. C
  205. C     Mise a jour de la solution X et du residu R
  206.       DO I=1,N
  207.         X.PROG(I) = X.PROG(I) + ALPHA*P.PROG(I) + OMEGA*S.PROG(I)
  208.         R.PROG(I) = S.PROG(I) - OMEGA * T.PROG(I)
  209.       ENDDO
  210. C
  211. C     Second test de convergence : ||R|| / ||B|| < XCRIT
  212.       SNORM = DDOTPV(N,R.PROG(1),R.PROG(1))
  213.       IF (SQRT(SNORM) / BNORM .LT. XCRIT) GOTO 100
  214.  
  215.       IF (ABS(OMEGA).LT.XPETIT) THEN
  216.          KERRE = 4
  217.          GOTO 100
  218.       ENDIF
  219.  
  220.       RHOOLD = RHONEW
  221. C
  222.       IF (ITER.LT.MAXIT) GOTO 10
  223. C----------------------------------------------------------------------
  224. C     MAXIT atteint sans converger
  225.       KERRE = 5
  226. C
  227.  100  CONTINUE
  228. C
  229.       IF (BTHRD) CALL THREADIS
  230. C
  231.  200  CONTINUE
  232. C
  233. C     Impression
  234.       IF (IIMPI.GT.1) THEN
  235.         IF (KERRE.EQ.0) THEN
  236.           WRITE(IOIMP,*) 'Convergence en ',ITER, ' iterations.'
  237.         ELSE
  238.           WRITE(IOIMP,*) 'Erreur ',KERRE,' a l iteration ',ITER
  239.         ENDIF
  240.       ENDIF
  241. C
  242. C     Menage
  243.       SEGSUP,R,R0,P,S,T,V
  244. C
  245.       END
  246.  
  247.  

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