C GRACO6 SOURCE PV090527 25/06/27 21:15:02 12303 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 PPARAM -INC CCOPTIO -INC CCASSIS -INC SMMATRI -INC SMVECTD segment ilicre * stockage matrice initiale en creux * ilideb position debut de ligne dans ligcre integer ilideb(nbinc+1) integer iliinc(ino+1) integer ligcrp endsegment segment ligcre * lmatr: longueur reelle ligne * posm: numero inconnue * valm: valeur terme integer posm(lmat) real*8 valm(lmat) endsegment POINTEUR MVECT4.MVECTD,MVECT5.MVECTD,MVECT6.MVECTD, $ MVECT7.MVECTD nbthr=nbthrs 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) * du fait des imprecisions nume©riques on peut mettre plus que INC iterations MAXIT=INC-1+1000 MILIGN=IILIGN SEGACT MILIGN SEGINI MVECT1, MVECT3,MVECT4 segini mvect2,mvect6 * conversion matrice en structure ligne complete creuse call graco11(mmatri,ilicre,0) * idem matrice factorisee call graco12(mmatri,ifacre,ifatra) * * on fait une premiere pseudo resolution pour avoir un ordre de * grandeur de l'energie * MVECTD=MSECO DO 31 IB=1,INC MVECT1.VECTBB(IB)=VECTBB(IB)*DNOR(IB) 31 CONTINUE CALL GRACO8 ( ICHOLX,MVECT1,NOEN,ifacre,ifatra) SEGACT MILIGN *p DO 38 IB=1,INC *p 38 CONTINUE CALL GRACO7(ilicre,MVECT1,MVECT2,inc,nbthr,lenb) ENEI = DDOT2(INC,MVECT1.VECTBB(1),MVECT2.VECTBB(1)) IF (IIMPI.NE.0) THEN WRITE(IOIMP,*) ' premiere estimation de l energie',ENEI ENDIF IF (ENEI.LE.0.D0) THEN write(ioimp,*) 'ATTENTION : ENEI = ',ENEI,' !' * INTERR(1) = 0 * CALL ERREUR(49) * if (nbthr.gt.1) call threadis * RETURN ENDIF 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 MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB) 59 CONTINUE CALL GRACO8(ICHOLX,MVECT3,NOEN,ifacre,ifatra) DO 60 IB=1,INC MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB) 60 CONTINUE C C DEBUT DE TOURNER EN ROND : MAXIT = INC-1 C IF (IIMPI.NE.0) THEN WRITE(ioimp,*)' Debut des iterations (max =',MAXIT,')' ENDIF C DO 100 iter = 1, MAXIT IF (IERR.NE.0) then if (nbthr.gt.1) call threadis RETURN endif * * Recalcul du critere et du residu toutes les 100 iterations : IF (iter/100 * 100 - iter .EQ. 0) then SEGACT MILIGN CALL GRACO7( ilicre,MVECT1,MVECT6,inc,nbthr,lenb) ENEN = DDOT2(INC,MVECT1.VECTBB(1),MVECT6.VECTBB(1)) IF (IIMPI.NE.0) THEN WRITE(ioimp,*)' Nouvelle energie ',ENEN,crit,'(iter=',iter,')' ENDIF IF (ENEN.LE.0.D0) THEN write(ioimp,*) 'ATTENTION : ENEN = ',ENEN,' !' * INTERR(1) = 0 * CALL ERREUR(49) * if (nbthr.gt.1) call threadis * RETURN ENDIF ENEI = ENEN IF (IIMPI.NE.0) THEN if (iter/10 * 10 - iter .eq. 0) then 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(Ilicre,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(ioimp,*) ' 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 MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB) 61 CONTINUE CALL GRACO8(ICHOLX,MVECT3,NOEN,ifacre,ifatra) R1RR1 = DDOT2(INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1)) CRIT = R1RR1 / ENEI IF (IIMPI.GT.1) WRITE(ioimp,*) ' critere ' , crit,iter IF (ABS(CRIT) .LT. 1.d-20 .OR. ITER.GT.80000) THEN * on a converge (mais en fait pas toujours si ITER est trop grand !) GO TO 101 ENDIF IF (R0RR0.EQ.0.D0) R0RR0=1d-30 IF (ABS(R0RR0).lt.1d-45) write(ioimp,*) ' r0rr0 ',r0rr0 *pv IF (R0RR0.EQ.0.D0) THEN *pv CALL ERREUR(835) *pVC WRITE(IOIMP,FMT='('' PROBLEME DE DIVISION PAR ZERO '')') *pv RETURN *pv ENDIF BET1 = R1RR1 / R0RR0 DO 9 IB=1, INC MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB)+BET1*MVECT4.VECTBB(IB) 9 CONTINUE R0RR0 = R1RR1 100 CONTINUE * Fin des iterations ! if (nbthr.gt.1) call threadis CALL ERREUR(460) C WRITE( IOIMP,FMT='('' PAS DE CONVERGENCE'')') RETURN 101 MSOL = MVECT1 if (nbthr.gt.1) call threadis if (iter.gt.80000) then write(ioimp,*) 'pas de convergence en ',iter,' iterations',crit else write(ioimp,*) ' convergence en ',iter, ' iterations',crit endif DO 67 IB=1,INC MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)*DNOR(IB) 67 CONTINUE SEGSUP MVECT2,MVECT3,MVECT4,MVECT6 * suppression stockage morse matrice et transposee SEGACT,MMATRI*mod ilicre=jlicre ligcre=ligcrp segsup ilicre,ligcre jlicre=0 ilicre=ifacre ligcre=ligcrp segsup ilicre,ligcre ifacre=0 ilicre=ifatra ligcre=ligcrp ifatra=0 segsup ilicre,ligcre SEGDES MDNOR,MMATRI,MILIGN RETURN END