1 : C GRACO6 SOURCE PV 13/10/14 21:17:48 7838 2 : SUBROUTINE GRACO6 ( ICHOLX, MSECO ,NOEN,MSOL,lenb) 3 : * 4 : * EXECUTE LES ITERATIONS DU GRADIENTS CONJUGUE. 5 : * 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8(A-H,O-Z) 8 : -INC CCOPTIO 9 : -INC SMMATRI 10 : -INC SMVECTD 11 : -INC CCASSIS 12 : POINTEUR MVECT4.MVECTD,MVECT5.MVECTD,MVECT6.MVECTD, 13 : $ MVECT7.MVECTD 14 : nbthr=max(nbesc,1) 15 : nbthr=min(nbthr,64) 15 : if (nbthr.gt.1) call threadii 16 : MMATRI=ICHOLX 17 : * activation de la matrice une fois pour toute. 18 : SEGACT,MMATRI 19 : MILIGN=IASLIG 20 : SEGACT,MILIGN 21 : INO=ILIGN(/1) 22 : MDNOR=IDNORM 23 : SEGACT MDNOR 24 : DO I=1,INO 25 : LLIGN=ILIGN(I) 26 : SEGACT LLIGN 27 : enddo 28 : MVECTD = MSECO 29 : SEGACT MVECTD 30 : INC = VECTBB(/1) 31 : MAXIT=INC- 1 32 : MILIGN=IILIGN 33 : SEGACT MILIGN 34 : SEGINI MVECT1, MVECT3,MVECT4 35 : if (nbthr.ge.1) then 36 : INC=INC*nbthr 37 : segini mvect2,mvect6 38 : INC=vectbb(/1) 39 : else 40 : segini mvect2,mvect6 41 : endif 42 : C IX0=MVECT1 43 : C IR0 =MVECT2 44 : C IRR0=MVECT3 45 : C IP0=MVECT4 46 : C IAP0=MVECT6 47 : * 48 : * on fait une premire pseudo resolution pour avoir un ordre de 49 : * grandeur de l'energie 50 : * 51 : MVECTD=MSECO 52 : DO 31 IB=1,INC 53 : 31 MVECT1.VECTBB(IB)=VECTBB(IB)*DNOR(IB) 53 : CALL GRACO8 ( ICHOLX,MVECT1,NOEN) 54 : SEGACT MILIGN 55 : *p DO 38 IB=1,INC 56 : *p 38 CONTINUE 56 : CALL GRACO7( ICHOLX,MVECT1,MVECT2,inc,nbthr,lenb) 57 : ENEI=DDOT2 (INC,MVECT1.VECTBB(1),MVECT2.VECTBB(1)) 58 : IF(IIMPI.NE.0)WRITE(IOIMP,*) 59 : * ' premiere estimation de l energie', ENEI 60 : MVECT7=MSECO 61 : DO 3 IB=1,INC 62 : MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT2.VECTBB(IB) 63 : 3 CONTINUE 64 : ** SEGDES MVECT7 65 : DO 59 IB=1,INC 66 : 59 MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB) 66 : CALL GRACO8 ( ICHOLX,MVECT3,NOEN) 67 : DO 60 IB=1,INC 68 : 60 MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB) 69 : C 70 : C DEBUT DE TOURNER EN ROND 71 : C 72 : DO 100 ITER=1,MAXIT 73 : if (ierr.ne.0) return 74 : IF(iter/100 * 100 - iter . EQ . 0) then 75 : * 76 : * recalcul du critere et du residu 77 : * 78 : SEGACT MILIGN 79 : *p DO 138 IB=1,INC 80 : *p138 CONTINUE 80 : CALL GRACO7( ICHOLX,MVECT1,MVECT6,inc,nbthr,lenb) 81 : ENEN=DDOT2 (INC,MVECT1.VECTBB(1),MVECT6.VECTBB(1)) 82 : WRITE(6,*)' nouvelle energie',ENEN,crit 83 : ENEI = ENEN 84 : IF(IIMPI.NE.0) THEN 85 : if(iter/10 *10 - iter .eq.0) then 86 : IF(IIMPI.NE.0) 87 : * WRITE(IOIMP,FMT='('' ITERATION '',I6,'' CRITERE '',E12.5)') 88 : * ITER ,CRIT 89 : ENDIF 90 : ENDIF 91 : * recalcul residu reel 92 : * IF(iter/1000 * 1000 - iter . EQ . 0) then 93 : * DO IB=1,INC 94 : * MVECT2.VECTBB(IB)=MVECT7.VECTBB(IB)*DNOR(IB)-MVECT6.VECTBB(IB) 95 : * ENDDO 96 : * ENDIF 97 : ENDIF 98 : * fin recalcul 98 : CALL GRACO7(ICHOLX,MVECT4,MVECT6,inc,nbthr,lenb) 99 : IF(ITER.EQ.1) then 100 : R0RR0=DDOT2 (INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1)) 101 : endif 102 : P0AP0=DDOT2 (INC,MVECT4.VECTBB(1),MVECT6.VECTBB(1)) 103 : IF(P0AP0.EQ.0.D0) P0AP0=1d-50 104 : IF(abs(P0AP0).lt.1d-45) write (6,*) ' p0ap0 ',p0ap0 105 : ALP0=R0RR0 /P0AP0 106 : DO 6 IB = 1, INC 107 : MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)+ALP0*MVECT4.VECTBB(IB) 108 : 6 CONTINUE 109 : DO 7 IB=1,INC 110 : MVECT2.VECTBB(IB)=MVECT2.VECTBB(IB)-ALP0*MVECT6.VECTBB(IB) 111 : 7 CONTINUE 112 : DO 61 IB=1,INC 113 : 61 MVECT3.VECTBB(IB)=MVECT2.VECTBB(IB) 113 : CALL GRACO8 ( ICHOLX,MVECT3,NOEN) 114 : R1RR1=DDOT2 (INC,MVECT2.VECTBB(1),MVECT3.VECTBB(1)) 115 : CRIT = R1RR1 / ENEI 116 : IF(IIMPI.GT.1) WRITE(6,*) ' critere ' , crit 117 : IF(ABS(CRIT) . LT . 1.d-20.OR.ITER.GT.80000) THEN 118 : * on a converge 119 : GO TO 101 120 : ENDIF 121 : IF(R0RR0.EQ.0.D0) R0RR0=1d-30 122 : IF(abs(R0RR0).lt.1d-45) write (6,*) ' r0rr0 ',r0rr0 123 : *pv IF(R0RR0.EQ.0.D0) GO TO 1000 124 : BET1 = R1RR1 / R0RR0 125 : DO 9 IB=1, INC 126 : MVECT4.VECTBB(IB)=MVECT3.VECTBB(IB)+BET1*MVECT4.VECTBB(IB) 127 : 9 CONTINUE 128 : R0RR0 = R1RR1 129 : 100 CONTINUE 129 : if (nbthr.gt.1) call threadis 129 : CALL ERREUR(462) 130 : C WRITE( IOIMP,FMT='('' PAS DE CONVERGENCE'')') 131 : RETURN 132 : 101 MSOL = MVECT1 132 : if (nbthr.gt.1) call threadis 133 : if (iter.ge.80000) then 134 : write (6,*) 'pas de convergence en ',iter,' iterations',crit 135 : else 136 : write(6,*) ' convergence en ',iter, ' iterations',crit 137 : endif 138 : DO 67 IB=1,INC 139 : 67 MVECT1.VECTBB(IB)=MVECT1.VECTBB(IB)*DNOR(IB) 140 : SEGDES MDNOR,MMATRI,MILIGN 141 : SEGSUP MVECT2,MVECT3,MVECT4, MVECT6 142 : RETURN 143 : 1000 CONTINUE 143 : CALL ERREUR(835) 144 : C WRITE(IOIMP,FMT='('' PROBLEME DE DIVISION PAR ZERO '')') 145 : RETURN 146 : END 147 : 148 : 149 : 150 : 151 : 152 : 153 : 154 : 155 : 156 : 157 : 158 : 159 : 160 :
© Cast3M 2003 - All rights reserved.
Disclaimer