C BCGSTB SOURCE MB234859 26/09/22 21:15:05 12653 SUBROUTINE BCGSTB(N,A,B,X,XCRIT,MAXIT,KERRE) C---------------------------------------------------------------------- C RESOLUTION DE AX=B PAR L'ALGORITHME BiCGSTAB C C Reference : "Templates for the Solution of Linear Systems : C Building Blocks for Iterative Methods", R.Barrett et al. C C Entrees : C --------- C N : ENTIER donnant la taille de la matrice et des vecteurs C A : MLENTI stockant la matrice par lignes. C Chaque ligne i a ses valeurs contenues dans un MLREEL C dont le pointeur est stocke dans A.LECT(i). C B : MLREEL donnant le vecteur second membre C MAXIT : ENTIER donnant le nombre d'iteration max souhaite C XPREC : REEL donnant la precision souhaitee C C Sorties : C --------- C X : MLREEL donnant le vecteur solution C KERRE : ENTIER precisant si tout s'est bien passe (0) ou non. C---------------------------------------------------------------------- C DECLARATIONS C---------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C -INC PPARAM -INC CCASSIS -INC CCOPTIO -INC CCREEL -INC SMLENTI POINTEUR A.MLENTI -INC SMLREEL POINTEUR X.MLREEL,B.MLREEL POINTEUR R.MLREEL,R0.MLREEL,P.MLREEL,S.MLREEL,T.MLREEL,V.MLREEL C REAL*8 ALPHA,BETA,OMEGA,RHOOLD,RHONEW,BNORM,SNORM,TNORM,VPSCA C C Informations pour la parallelisation LOGICAL BTHRD,BSGDES COMMON/PDTSCA/NBTHRD,NVAL,NLIG,IMAT,IVEC,IRES EXTERNAL BCGSi C---------------------------------------------------------------------- C INITIALISATIONS C---------------------------------------------------------------------- KERRE = 0 ITER = 0 ALPHA = 0.0D0 BETA = 0.0D0 OMEGA = 1.0D0 RHOOLD = 1.0D0 C JG=N SEGINI,R,R0,P,S,T,V C NLA=0 BSGDES=.FALSE. DO I=1,N MLREE1 = A.LECT(I) SEGACT /ERR=3/ MLREE1 NLA=NLA+1 ENDDO GOTO 4 C 3 CONTINUE BSGDES=.TRUE. MLREE1 = A.LECT(NLA) SEGDES,MLREE1 NLA=NLA-1 4 CONTINUE C C Norme du second membre b BNORM = DDOTPV(N,B.PROG(1),B.PROG(1)) BNORM = SQRT(BNORM) IF (BNORM.LE.XPETIT) BNORM = 1.0D0 C C Residu initial R0 = b - Ax0 DO I=1,NLA MLREE1 = A.LECT(I) R.PROG(I) = B.PROG(I) - DDOTPV(N,MLREE1.PROG(1),X.PROG(1)) ENDDO IF (BSGDES) THEN DO I=NLA+1,N MLREE1 = A.LECT(I) SEGACT,MLREE1 R.PROG(I) = B.PROG(I) - DDOTPV(N,MLREE1.PROG(1),X.PROG(1)) SEGDES,MLREE1 ENDDO ENDIF C C Convergence atteinte? SNORM = DDOTPV(N,R.PROG(1),R.PROG(1)) IF (SQRT(SNORM)/BNORM .LT. XCRIT) GOTO 200 C DO I=1,N R0.PROG(I) = R.PROG(I) P .PROG(I) = R.PROG(I) ENDDO C C Parallelisation des produits matrice-vecteurs NBTHR = MIN(MAX(N*NLA/15000000,1),NBTHRS) BTHRD = .TRUE. IF ((NBTHRS.EQ.1).OR.(OOTHRD.GT.0)) THEN NBTHR = 1 BTHRD = .FALSE. ENDIF IF (BTHRD) CALL THREADII NBTHRD=NBTHR NVAL=N NLIG=NLA IMAT=A C---------------------------------------------------------------------- C BOUCLE PRINCIPALE DES ITERATIONS C---------------------------------------------------------------------- 10 CONTINUE ITER = ITER + 1 IF (IERR.NE.0) GOTO 100 C C Mise a jour de RHONEW = RHONEW = DDOTPV(N,R0.PROG(1),R.PROG(1)) IF (ABS(RHONEW).LT.XPETIT) THEN KERRE = 1 GOTO 100 ENDIF C C Coefficient de direction + maj de la direction de recherche IF (ITER.GT.1) THEN BETA = (RHONEW / RHOOLD) * (ALPHA / OMEGA) DO I=1,N P.PROG(I) = R.PROG(I) + BETA * (P.PROG(I) - OMEGA * V.PROG(I)) ENDDO ENDIF C C Calcul de V=AP IVEC=P IRES=V DO ITH=2,NBTHR CALL THREADID(ITH,BCGSI) ENDDO CALL BCGSI(1) DO ITH=2,NBTHR CALL THREADIF(ITH) ENDDO IF (BSGDES) THEN DO I=NLA+1,N MLREE1 = A.LECT(I) SEGACT,MLREE1 V.PROG(I) = DDOTPV(N,MLREE1.PROG(1),X.PROG(1)) SEGDES,MLREE1 ENDDO ENDIF C C Calcul du pas ALPHA = RHONEW / VPSCA avec VPSCA = VPSCA = DDOTPV(N,R0.PROG(1),V.PROG(1)) IF (ABS(VPSCA).LT.XPETIT) THEN KERRE = 2 GOTO 100 ENDIF ALPHA = RHONEW / VPSCA C C Calcul du résidu intermédiaire S = R - ALPHA * V + test de convergence DO I=1,N S.PROG(I) = R.PROG(I) - ALPHA * V.PROG(I) ENDDO C C Premier test de convergence : ||S|| / ||B|| < XCRIT SNORM = DDOTPV(N,S.PROG(1),S.PROG(1)) IF (SQRT(SNORM) / BNORM .LT. XCRIT) THEN DO I=1,N X.PROG(I) = X.PROG(I) + ALPHA * P.PROG(I) ENDDO GOTO 100 ENDIF C C Calcul du second pas T = AS IVEC=S IRES=T DO ITH=2,NBTHR CALL THREADID(ITH,BCGSI) ENDDO CALL BCGSI(1) DO ITH=2,NBTHR CALL THREADIF(ITH) ENDDO IF (BSGDES) THEN DO I=NLA+1,N MLREE1 = A.LECT(I) SEGACT,MLREE1 V.PROG(I) = DDOTPV(N,MLREE1.PROG(1),X.PROG(1)) SEGDES,MLREE1 ENDDO ENDIF C C Parametre de stabilisation OMEGA = / TNORM = DDOTPV(N,T.PROG(1),T.PROG(1)) VPSCA = DDOTPV(N,S.PROG(1),T.PROG(1)) IF (TNORM.EQ.0.0D0) THEN KERRE = 3 GOTO 100 ENDIF OMEGA = VPSCA / TNORM C C Mise a jour de la solution X et du residu R DO I=1,N X.PROG(I) = X.PROG(I) + ALPHA*P.PROG(I) + OMEGA*S.PROG(I) R.PROG(I) = S.PROG(I) - OMEGA * T.PROG(I) ENDDO C C Second test de convergence : ||R|| / ||B|| < XCRIT SNORM = DDOTPV(N,R.PROG(1),R.PROG(1)) IF (SQRT(SNORM) / BNORM .LT. XCRIT) GOTO 100 IF (ABS(OMEGA).LT.XPETIT) THEN KERRE = 4 GOTO 100 ENDIF RHOOLD = RHONEW C IF (ITER.LT.MAXIT) GOTO 10 C---------------------------------------------------------------------- C MAXIT atteint sans converger KERRE = 5 C 100 CONTINUE C IF (BTHRD) CALL THREADIS C 200 CONTINUE C C Impression IF (IIMPI.GT.1) THEN IF (KERRE.EQ.0) THEN WRITE(IOIMP,*) 'Convergence en ',ITER, ' iterations.' ELSE WRITE(IOIMP,*) 'Erreur ',KERRE,' a l iteration ',ITER ENDIF ENDIF C C Menage SEGSUP,R,R0,P,S,T,V C END