C GRACO6 SOURCE PV 13/10/14 21:17:48 7838 SUBROUTINE GRACO6 ( ICHOLX, MSECO ,NOEN,MSOL,lenb) * * EXECUTE LES ITERATIONS DU GRADIENTS CONJUGUE. * IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC CCOPTIO -INC SMMATRI -INC SMVECTD -INC CCASSIS POINTEUR MVECT4.MVECTD,MVECT5.MVECTD,MVECT6.MVECTD, $ MVECT7.MVECTD nbthr=max(nbesc,1) nbthr=min(nbthr,64) if (nbthr.gt.1) call threadii MMATRI=ICHOLX * activation de la matrice une fois pour toute. SEGACT,MMATRI MILIGN=IASLIG SEGACT,MILIGN INO=ILIGN(/1) MDNOR=IDNORM SEGACT MDNOR DO I=1,INO LLIGN=ILIGN(I) SEGACT LLIGN enddo MVECTD = MSECO SEGACT MVECTD INC = VECTBB(/1) MAXIT=INC- 1 MILIGN=IILIGN SEGACT MILIGN SEGINI MVECT1, MVECT3,MVECT4 if (nbthr.ge.1) then INC=INC*nbthr segini mvect2,mvect6 INC=vectbb(/1) else segini mvect2,mvect6 endif C IX0=MVECT1 C IR0 =MVECT2 C IRR0=MVECT3 C IP0=MVECT4 C IAP0=MVECT6 * * on fait une premire pseudo resolution pour avoir un ordre de * grandeur de l'energie * MVECTD=MSECO DO 31 IB=1,INC 31 MVECT1.VECTBB(IB)=VECTBB(IB)*DNOR(IB) CALL GRACO8 ( ICHOLX,MVECT1,NOEN) SEGACT MILIGN *p DO 38 IB=1,INC *p 38 CONTINUE CALL GRACO7( ICHOLX,MVECT1,MVECT2,inc,nbthr,lenb) ENEI=DDOT2 (INC,MVECT1.VECTBB(1),MVECT2.VECTBB(1)) IF(IIMPI.NE.0)WRITE(IOIMP,*) * ' premiere estimation de l energie', ENEI MVECT7=MSECO DO 3 IB=1,INC MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT2.VECTBB(IB) 3 CONTINUE ** SEGDES MVECT7 DO 59 IB=1,INC 59 MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB) CALL GRACO8 ( ICHOLX,MVECT3,NOEN) DO 60 IB=1,INC 60 MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB) C C DEBUT DE TOURNER EN ROND C DO 100 ITER=1,MAXIT if (ierr.ne.0) return IF(iter/100 * 100 - iter . EQ . 0) then * * recalcul du critere et du residu * SEGACT MILIGN *p DO 138 IB=1,INC *p138 CONTINUE CALL GRACO7( ICHOLX,MVECT1,MVECT6,inc,nbthr,lenb) ENEN=DDOT2 (INC,MVECT1.VECTBB(1),MVECT6.VECTBB(1)) WRITE(6,*)' nouvelle energie',ENEN,crit ENEI = ENEN IF(IIMPI.NE.0) THEN if(iter/10 *10 - iter .eq.0) then IF(IIMPI.NE.0) * WRITE(IOIMP,FMT='('' ITERATION '',I6,'' CRITERE '',E12.5)') * ITER ,CRIT ENDIF ENDIF * recalcul residu reel * IF(iter/1000 * 1000 - iter . EQ . 0) then * DO IB=1,INC * MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT6.VECTBB(IB) * ENDDO * ENDIF ENDIF * fin recalcul CALL GRACO7(ICHOLX,MVECT4,MVECT6,inc,nbthr,lenb) IF(ITER.EQ.1) then R0RR0=DDOT2 (INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1)) endif P0AP0=DDOT2 (INC,MVECT4.VECTBB(1),MVECT6.VECTBB(1)) IF(P0AP0.EQ.0.D0) P0AP0=1d-50 IF(abs(P0AP0).lt.1d-45) write (6,*) ' p0ap0 ',p0ap0 ALP0=R0RR0 /P0AP0 DO 6 IB = 1, INC MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)+ALP0*MVECT4.VECTBB(IB) 6 CONTINUE DO 7 IB=1,INC MVECT2.VECTBB(IB)=MVECT2.VECTBB(IB)-ALP0*MVECT6.VECTBB(IB) 7 CONTINUE DO 61 IB=1,INC 61 MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB) CALL GRACO8 ( ICHOLX,MVECT3,NOEN) R1RR1=DDOT2 (INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1)) CRIT = R1RR1 / ENEI IF(IIMPI.GT.1) WRITE(6,*) ' critere ' , crit IF(ABS(CRIT) . LT . 1.d-20.OR.ITER.GT.80000) THEN * on a converge GO TO 101 ENDIF IF(R0RR0.EQ.0.D0) R0RR0=1d-30 IF(abs(R0RR0).lt.1d-45) write (6,*) ' r0rr0 ',r0rr0 *pv IF(R0RR0.EQ.0.D0) GO TO 1000 BET1 = R1RR1 / R0RR0 DO 9 IB=1, INC MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB)+BET1*MVECT4.VECTBB(IB) 9 CONTINUE R0RR0 = R1RR1 100 CONTINUE if (nbthr.gt.1) call threadis CALL ERREUR(462) C WRITE( IOIMP,FMT='('' PAS DE CONVERGENCE'')') RETURN 101 MSOL = MVECT1 if (nbthr.gt.1) call threadis if (iter.ge.80000) then write (6,*) 'pas de convergence en ',iter,' iterations',crit else write(6,*) ' convergence en ',iter, ' iterations',crit endif DO 67 IB=1,INC 67 MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)*DNOR(IB) SEGDES MDNOR,MMATRI,MILIGN SEGSUP MVECT2,MVECT3,MVECT4, MVECT6 RETURN 1000 CONTINUE CALL ERREUR(835) C WRITE(IOIMP,FMT='('' PROBLEME DE DIVISION PAR ZERO '')') RETURN END