IMPLICIT REAL*8 (A-H,O-Z) REAL*4 CPUT PARAMETER(MJ=100,M=MJ*10,M2=MJ*11-1,N=M*M2,EPS=1.0D-7) DIMENSION AA(N),AB(-M:N),AC(-M:N),F(N), WK(M) +, AB1(-M:N+M),AC1(-M:N+M) +, FQ(N+M) +, B1(N),BW(-M:N+M) +, DI(-M:N),P(-M:N+M),W(-M:N+M) +, X(-M:N+M),AP(N), UB(-M:N),AE(-M:N) +, AF(-M:N),AF1(-M:N+M), AE1(-M:N+M), UB1(-M:N+M) +, V(-M:N+M) +, RRA(-M:N+M),PA(-M:N+M),RR(-M:N+M), R0(N) +, R(N) C DATA EPSN/1.0D-5/,NIT/128/ EPSICW=EPSN/16.0 USHIRB=1.01D0 UU=0.95D0 * 0.98(Best)*** ISTB=3 ISISU=1 C J1=M*(5*MJ-1) +3*MJ J2=M*(6*MJ-1) +3*MJ J3=M*(5*MJ-1) +8*MJ J4=M*(6*MJ-1) +8*MJ C READ(5,*) DF WRITE(6,*) ' % DF= ',DF,' ISISU= ',ISISU WRITE(6,*) ' % M= ',M,' M2= ',M2,' N= ',N,' : ISTB= ',ISTB FMJ= DFLOAT(MJ) DHX=1.0D0/FMJ DHY=1.0D0/FMJ CC C0=10.0D0 WRITE(6,*) ' % C0= ',C0 DH=1.0D0/FMJ C0H=C0*DH C0H2=0.5D0*C0H CCCCC IF(ISISU.EQ.1) THEN CCCCC IF(DF .EQ. 0.0D0) THEN C0H=0.0D0 ENDIF DO 15 I=1,M-1 AA(I)= 0.5D0*(2.0D0+2.0D0*DF +DF*BER4(C0H)+ BER4(C0H) + + BER4(-C0H)+ DF*BER4(-C0H)) AB(I)=-0.5D0*(BER4(C0H)+ DF*BER4(C0H)) AC(I)=-0.5D0*2.0D0 AB1(I)=-0.5D0*(DF*BER4(-C0H)+ BER4(-C0H)) AC1(I)=0.0D0 15 CONTINUE AB1(1)=0.0D0 * AA(M)= 0.5D0*(DF*BER4(C0H)+ BER4(C0H)+1.0D0 +DF) AB(M)=0.0D0 AC(M)=-0.5D0*1.0D0 AB1(M)=-0.5D0*(DF*BER4(-C0H)+ BER4(-C0H)) AC1(M)=0.0D0 C C0H=C0*DH DO 35 J=1,M2-2 DO 25 K=1,M-1 AA(M*J+K)=0.5D0*(2.0D0*BER4(C0H)+ 2.0D0*BER4(-C0H)+4.0D0) AB(M*J+K)=-0.5D0*(2.0D0*BER4(C0H)) AC(M*J+K)=-0.5D0*2.0D0 AB1(M*J+K)=-0.5D0*(2.0D0*BER4(-C0H)) AC1(M*J+K)=-0.5D0*2.0D0 25 CONTINUE AB1(M*J+1)=0.0D0 * AA(M*J+M)= 0.5D0*(2.0D0 +2.0D0*BER4(C0H)) AB(M*J+M)=0.0D0 AC(M*J+M)=-0.5D0*1.0D0 AB1(M*J+M)=-0.5D0*(2.0D0*BER4(-C0H)) AC1(M*J+M)=-0.5D0*1.0D0 35 CONTINUE C IF(DF .EQ. 0.0D0) THEN C0H=0.0D0 ENDIF DO 45 I=N-M+1, N-1 AA(I)= 0.5D0*(2.0D0+2.0D0*DF +BER4(C0H)+ DF*BER4(C0H) + + DF*BER4(-C0H)+ BER4(-C0H)) AB(I)=-0.5D0*(DF*BER4(C0H)+ BER4(C0H)) AC(I)=0.0D0 AB1(I)=-0.5D0*(BER4(-C0H)+ DF*BER4(-C0H)) AC1(I)=-0.5D0*2.0D0 45 CONTINUE AB1(N-M+1)=0.0D0 * AA(N)= 0.5D0*(BER4(C0H)+ DF*BER4(C0H)+ DF+ 1.0D0) AB(N)=0.0D0 AC(N)=0.0D0 AB1(N)=-0.5D0*(BER4(-C0H)+ DF*BER4(-C0H)) AC1(N)=-0.5D0*1.0D0 CCCCC ELSE CCCCC IF(DF .EQ. 0.0D0) THEN C0H2=0.0D0 ENDIF DO 10 I=1,M-1 AA(I)= 0.5D0*(DF*4.0D0+4.0D0) AB(I)=-0.5D0*((1.0D0-C0H2)+DF*(1.0D0-C0H2)) AC(I)=-0.5D0*2.0D0 AB1(I)=-0.5D0*(DF*(1.0D0+C0H2)+(1.0D0+C0H2)) AC1(I)=0.0D0 10 CONTINUE AB1(1)=0.0D0 * AA(M)= 0.5D0*(DF*(2.0D0-C0H2)+(2.0D0-C0H2)) AB(M)=0.0D0 AC(M)=-0.5D0*1.0D0 AB1(M)=-0.5D0*(DF*(1.0D0+C0H2)+(1.0D0+C0H2)) AC1(M)=0.0D0 C C0H2=0.5D0*C0H DO 30 J=1,M2-2 DO 20 K=1,M-1 AA(M*J+K)=4.0D0 AB(M*J+K)=-0.5D0*(2.0D0-C0H) AC(M*J+K)=-0.5D0*2.0D0 AB1(M*J+K)=-0.5D0*(2.0D0+C0H) AC1(M*J+K)=-0.5D0*2.0D0 20 CONTINUE AB1(M*J+1)=0.0D0 * AA(M*J+M)= 0.5D0*(4.0D0-C0H) AB(M*J+M)=0.0D0 AC(M*J+M)=-0.5D0*1.0D0 AB1(M*J+M)=-0.5D0*(2.0D0+C0H) AC1(M*J+M)=-0.5D0*1.0D0 30 CONTINUE C IF(DF .EQ. 0.0D0) THEN C0H2=0.0D0 ENDIF DO 40 I=N-M+1, N-1 AA(I)= 0.5D0*(DF*4.0D0+4.0D0) AB(I)=-0.5D0*(DF*(1.0D0-C0H2)+(1.0D0-C0H2)) AC(I)=0.0D0 AB1(I)=-0.5D0*((1.0D0+C0H2)+DF*(1.0D0+C0H2)) AC1(I)=-0.5D0*2.0D0 40 CONTINUE AB1(N-M+1)=0.0D0 * AA(N)= 0.5D0*(DF*(2.0D0-C0H2)+(2.0D0-C0H2)) AB(N)=0.0D0 AC(N)=0.0D0 AB1(N)=-0.5D0*((1.0D0+C0H2)+DF*(1.0D0+C0H2)) AC1(N)=-0.5D0*1.0D0 CCCCC ENDIF CCCCC *C Uhen *** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 65 J= 5*MJ, 6*MJ-1 C DO 75 I=1,2*MJ C FQ(M*J +2*MJ+I)= 0.4D0 C FQ(M*J +6*MJ+I)=-0.4D0 C 75 CONTINUE C 65 CONTINUE CCC C DO 85 J= 5*MJ-1, 6*MJ-1 C DO 95 I=0,2*MJ C KP=M*J +2*MJ+I C F(KP)=0.25D0*(DHX*DHY)*(FQ(KP)+FQ(KP+M)+FQ(KP+1)+FQ(KP+1+M)) C KM=M*J +6*MJ+I C F(KM)=0.25D0*(DHX*DHY)*(FQ(KM)+FQ(KM+M)+FQ(KM+1)+FQ(KM+1+M)) C 95 CONTINUE C 85 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 60 J= 5*MJ-1, 6*MJ-1 DO 70 I=0,2*MJ F(M*J +2*MJ+I)= 0.2D0*(DHX*DHY) F(M*J +6*MJ+I)=-0.2D0*(DHX*DHY) 70 CONTINUE 60 CONTINUE C CALL CLOCK0 C *C BW(Ui) no Syokiti *** DO 111 K=1,N 111 BW(K)=0.0D0 C *C A*U=F wo Toku; Kotae --> F *** BNRM=0.0 DO 110 I=1,N 110 BNRM=BNRM+F(I)*F(I) BNRM=DSQRT(BNRM) * IF(ISTB.EQ.3) THEN CALL UPILU4(N,M,AA,AB,AC,AE,AF,UB,DI, USHIRB + ,AB1,AC1,AE1,AF1,UB1) CC CALL ILUCM3(N,M,AA,AB,AC,AE,AF,UB,DI, UU CC + ,AB1,AC1,AE1,AF1,UB1) ELSE IF(ISTB.EQ.2) THEN CALL UPILU2(N,M,AA,AB,AC,AE,UB,DI, USHIRB + ,AB1,AC1,AE1,UB1) CC CALL ILUCM2(N,M,AA,AB,AC,AE,UB,DI, UU CC + ,AB1,AC1,AE1,UB1) ENDIF C DO 226 I=1,N 226 X(I)=BW(I) C KCT=0 DO 230 IN=1,NIT IF(ISTB.EQ.3) THEN CALL BCGSTB3(AA,AB,AC,AE,AF,UB,DI,AB1,AC1,AE1,AF1,UB1,P,V,W,F + ,R0 ,X,AP,EPSICW,ERR,M,N,KCOUNT,BNRM,RRA,PA,RR) ELSE IF(ISTB.EQ.2) THEN CALL BCGSTB2(AA,AB,AC,AE, UB,DI,AB1,AC1,AE1, UB1,P,V,W,F + ,R0 ,X,AP,EPSICW,ERR,M,N,KCOUNT,BNRM,RRA,PA,RR) CC CALL BCG12( AA,AB,AC,AE, UB,DI,AB1,AC1,AE1, UB1,P,R,W,F CC + ,X,AP,EPSICW,ERR,M,N,KCOUNT,BNRM,RRA,PA,RR) ENDIF C KCT=KCT+KCOUNT CNR=0.0D0 DO 116 I=1,N B1(I)=F(I)-(AC1(I)*X(I-M)+AB1(I)*X(I-1) + +AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M)) CC B1(I)=F(I)-(AC(I-M)*X(I-M)+AB(I-1)*X(I-1) CC + +AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M)) 116 CNR=CNR+B1(I)*B1(I) CNR=DSQRT(CNR)/BNRM * WRITE(6,5800) IN,X(J1),X(J2),X(J3),X(J4),CNR,KCOUNT 5800 FORMAT(3H %,X,I4,F10.5,F10.5,F10.5,F10.5,D14.6,I5) IF(CNR.LT.EPSN) GO TO 250 EPSICW=EPSICW/2.0 C 230 CONTINUE 250 CONTINUE C DO 220 I=1,N *C BW(I)= X(I) *** 220 F(I) = X(I) C CC CALL BGLU( N, M, AR, EPS, IR, WK ) CC CALL BGSLV( N, M, AR, F ) C CALL CLOCK(CPUT) UMIN=F(1) JUMIN=1 KUMIN=1 UMAX=F(1) JUMAX=1 KUMAX=1 DO 130 J=1,M2 DO 120 K=1,M I=M*(J-1)+K IF(F(I).LT.UMIN) THEN UMIN=F(I) JUMIN=J KUMIN=K ENDIF IF(F(I).GT.UMAX) THEN UMAX=F(I) JUMAX=J KUMAX=K ENDIF 120 CONTINUE 130 CONTINUE C WRITE(6,*) ' % MJ= ',MJ WRITE(6,*) ' % CPUT= ',CPUT WRITE(6,*) ' % UMIN,(x,y)= ',UMIN,'(',JUMIN,',',KUMIN,')' WRITE(6,*) ' % UMAX,(x,y)= ',UMAX,'(',JUMAX,',',KUMAX,')' C WRITE(6,*) ' z=[ ' C CC DO 80 K= M,MJ,-MJ CC WRITE(6,2000) (F(M*J+K),J=MJ-1,MJ*10-1,MJ),' ; ' CC 80 CONTINUE *C MATLAB You *** DO 80 K= MJ,M,MJ WRITE(6,2000) (F(M*J+K),J=MJ-1,MJ*10-1,MJ),' ; ' 80 CONTINUE 2000 FORMAT(1H ,10F6.2,A3) C WRITE(6,*) ' ]; ' WRITE(6,*) ' contour(0:1:9, 0:1:9, z,10);' C STOP END CCC CCCCC* ** C *** FUJITSU REAL*8 FUNCTION BER4(BEX) REAL*8 BEX, HX C *** HITACHI C REAL FUNCTION BER4*8(BEX) C REAL*8 BER4, BEX, HX C** HX=1.000D-4 HX=1.000D-2 * IF(BEX.LT.-HX) THEN BER4=BEX/(DEXP(BEX)-1.0D0) ELSE IF(BEX.GE.-HX .AND. BEX.LE.HX) THEN BER4=1.0D0-(0.5D0-0.0833333*BEX)*BEX * CCCC BER4=1.0D0-0.5D0*BEX ELSE IF(BEX.GT.HX) THEN BER4=BEX*DEXP(-BEX)/(1.0D0-DEXP(-BEX)) ENDIF * END * CCC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * Aoki Aakusei(1,3)0.98 *** SUBROUTINE ILUCM3(N,M,AA,AB,AC,AE,AF,UB,DI,U,AB1,AC1,AE1,AF1,UB1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N),AF(-M:N),UB(-M:N), + AB1(-M:N+M), AC1(-M:N+M),AE1(-M:N+M),AF1(-M:N+M),UB1(-M:N+M) *** DO 100 I = 1, N *** CC W= ABS(AC1(I))+ABS(AB1(I))+ABS(AB(I))+ABS(AC(I)) CC UUI = MAX(W,AA(I)) *** AE1(I)= -AC1(I)*UB(I-M)*DI(I-M) AF1(I)= -AE1(I)*UB(I-M+1)*DI(I-M+1) UB1(I)= AB1(I)-AC1(I)*AE(I-M)*DI(I-M)-AE1(I)*AF(I-M+1)*DI(I-M+1) * C/ U: UTOP DIV=AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -AE1(I)*AE(I-M+1)*DI(I-M+1)-AF1(I)*AF(I-M+2)*DI(I-M+2) + -U*(AF1(I)*AC(I-M+2)*DI(I-M+2) + +AF(I-M)*AC1(I)*DI(I-M) + +UB1(I)*AF(I-1)*DI(I-1)+UB(I-M+2)*AF1(I)*DI(I-M+2)) C DIV= AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) C + -AE1(I)*AE(I-M+1)*DI(I-M+1)-AF1(I)*AF(I-M+2)*DI(I-M+2) DI(I)=1.0D0/DIV * UB(I)= AB(I)-AE1(I)*AC(I-M+1)*DI(I-M+1) + -AF1(I)*AE(I-M+2)*DI(I-M+2) AF(I)= -UB1(I)*AE(I-1)*DI(I-1) AE(I)= -UB1(I)*AC(I-1)*DI(I-1) * 100 CONTINUE RETURN END CCC * Aoki Aakusei(1,2) *** SUBROUTINE UPILU2(N,M,AA,AB,AC,AE,UB,DI,U,AB1,AC1,AE1,UB1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N),UB(-M:N), + AB1(-M:N+M), AC1(-M:N+M),AE1(-M:N+M),UB1(-M:N+M) *** DO 100 I = 1, N *** CC W= ABS(AC1(I))+ABS(AB1(I))+ABS(AB(I))+ABS(AC(I)) CC UUI = MAX(W,AA(I)) *** AE1(I)= -AC1(I)*UB(I-M)*DI(I-M) UB1(I)= AB1(I)-AC1(I)*AE(I-M)*DI(I-M) * C/ U: UTOP DIV=U*AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -AE1(I)*AE(I-M+1)*DI(I-M+1) C DIV= AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) C + -AE1(I)*AE(I-M+1)*DI(I-M+1) DI(I)=1.0D0/DIV * UB(I)= AB(I)-AE1(I)*AC(I-M+1)*DI(I-M+1) AE(I)= -UB1(I)*AC(I-1)*DI(I-1) * 100 CONTINUE RETURN END CCC SUBROUTINE UPILU4(N,M,AA,AB,AC,AE,AF,UB,DI,U,AB1,AC1,AE1,AF1,UB1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N),AF(-M:N),UB(-M:N), + AB1(-M:N+M), AC1(-M:N+M),AE1(-M:N+M),AF1(-M:N+M),UB1(-M:N+M) *** DO 100 I = 1, N *** CC W= ABS(AC1(I))+ABS(AB1(I))+ABS(AB(I))+ABS(AC(I)) CC UUI = MAX(W,AA(I)) *** AE1(I)= -AC1(I)*UB(I-M)*DI(I-M) AF1(I)= -AE1(I)*UB(I-M+1)*DI(I-M+1) UB1(I)= AB1(I)-AC1(I)*AE(I-M)*DI(I-M)-AE1(I)*AF(I-M+1)*DI(I-M+1) * C/ U: UTOP DIV=U*AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -AE1(I)*AE(I-M+1)*DI(I-M+1)-AF1(I)*AF(I-M+2)*DI(I-M+2) C DIV= AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) C + -AE1(I)*AE(I-M+1)*DI(I-M+1)-AF1(I)*AF(I-M+2)*DI(I-M+2) DI(I)=1.0D0/DIV * UB(I)= AB(I)-AE1(I)*AC(I-M+1)*DI(I-M+1) + -AF1(I)*AE(I-M+2)*DI(I-M+2) AF(I)= -UB1(I)*AE(I-1)*DI(I-1) AE(I)= -UB1(I)*AC(I-1)*DI(I-1) * 100 CONTINUE RETURN END CCC SUBROUTINE BCGSTB3(AA,AB,AC,AE,AF,UB,DI,AB1,AC1,AE1,AF1,UB1,P,V,W + ,B0 ,R0 ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM,AEE,EE,RR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),V(-M1:N+M1),UB(-M1:N),AE(-M1:N),AF(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, AF1(-M1:N+M1),AEE(-M1:N+M1),EE(-M1:N+M1),RR(-M1:N+M1),R0(N) * DO 10 I = 1, N 10 V(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB1(I)*X(I-1)+AC1(I)*X(I-M1)) *** KOKO DE WA V = B -AX DE ARU *** * DO 30 I = 1, N 30 W(I)=(V(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1) + -AF1(I)*W(I-M1+2))*DI(I) * DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1) + +AF(I)*W(I+M1-2)) *** KOKO DE WA W = INV(LU)*(B-AX) DE ARU ****** CR= 0.0D0 DO 60 I = 1, N P(I) = W(I) RR(I)= W(I) R0(I)= W(I) 60 CR = CR + W(I)*W(I) *** ITERATION; KONO ATO, V O V,W O Q TO SHITE TSUKAU *** DO 1000 K = 1, 400 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB1(I)*P(I-1)+AC1(I)*P(I-M1) 100 CONTINUE *** DO 110 I = 1, N 110 W(I)=(AP(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1) + -AF1(I)*W(I-M1+2))*DI(I) * DO 120 I = N, 1, -1 120 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1) + +AF(I)*W(I+M1-2)) *** W = INV(LU)*AP DE ARU *** SIG= 0.0D0 DO 130 I=1, N 130 SIG = SIG+R0(I)*W(I) ALFA = CR/SIG * DO 140 I=1,N 140 EE(I)= RR(I)-ALFA*W(I) * DO 150 I = 1, N 150 AEE(I)= AA(I)*EE(I)+AB(I)*EE(I+1)+AC(I)*EE(I+M1) + +AB1(I)*EE(I-1)+AC1(I)*EE(I-M1) *** DO 160 I = 1, N 160 V(I)=(AEE(I)-AC1(I)*V(I-M1)-UB1(I)*V(I-1)-AE1(I)*V(I-M1+1) + -AF1(I)*V(I-M1+2))*DI(I) * DO 170 I = N, 1, -1 170 V(I)=V(I)-DI(I)*(UB(I)*V(I+1)+AC(I)*V(I+M1)+AE(I)*V(I+M1-1) + +AF(I)*V(I+M1-2)) *** V = INV(LU)*A*E DE ARU *** EEV= 0.0D0 VV= 0.0D0 DO 180 I= 1,N EEV= EEV+EE(I)*V(I) 180 VV = VV +V(I)*V(I) OMEG =EEV/VV * DO 190 I= 1, N X(I)= X(I)+ALFA*P(I)+OMEG*EE(I) 190 RR(I)=EE(I)-OMEG*V(I) * CR=0.0D0 ERR = 0.0D0 DO 200 I = 1, N CR = CR+ R0(I)*RR(I) 200 ERR = ERR+RR(I)*RR(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 * BETA = CR/(OMEG*SIG) DO 300 I = 1, N 300 P(I)= RR(I)+BETA*(P(I)-OMEG*W(I)) * 1000 CONTINUE 2000 KCOUNT = K RETURN END * * * SUBROUTINE BCGSTB2(AA,AB,AC,AE,UB,DI,AB1,AC1,AE1,UB1,P,V,W,B0 + ,R0 ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM,AEE,EE,RR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),V(-M1:N+M1),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, AEE(-M1:N+M1),EE(-M1:N+M1),RR(-M1:N+M1),R0(N) * DO 10 I = 1, N 10 V(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB1(I)*X(I-1)+AC1(I)*X(I-M1)) *** KOKO DE WA V = B -AX DE ARU *** DO 30 I = 1, N 30 W(I)=(V(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) *** KOKO DE WA W = INV(LU)*(B-AX) DE ARU ****** CR= 0.0D0 DO 60 I = 1, N P(I) = W(I) RR(I)= W(I) R0(I)= W(I) 60 CR = CR + W(I)*W(I) *** ITERATION; KONO ATO, V O V,W O Q TO SHITE TSUKAU *** DO 1000 K = 1, 400 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB1(I)*P(I-1)+AC1(I)*P(I-M1) 100 CONTINUE DO 110 I = 1, N 110 W(I)=(AP(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 120 I = N, 1, -1 120 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) *** W = INV(LU)*AP DE ARU *** SIG= 0.0D0 DO 130 I=1, N 130 SIG = SIG+R0(I)*W(I) ALFA = CR/SIG * DO 140 I=1,N 140 EE(I)= RR(I)-ALFA*W(I) * * DO 150 I = 1, N 150 AEE(I)= AA(I)*EE(I)+AB(I)*EE(I+1)+AC(I)*EE(I+M1) + +AB1(I)*EE(I-1)+AC1(I)*EE(I-M1) DO 160 I = 1, N 160 V(I)=(AEE(I)-AC1(I)*V(I-M1)-UB1(I)*V(I-1)-AE1(I)*V(I-M1+1))*DI(I) DO 170 I = N, 1, -1 170 V(I)=V(I)-DI(I)*(UB(I)*V(I+1)+AC(I)*V(I+M1)+AE(I)*V(I+M1-1)) *** V = INV(LU)*A*E DE ARU *** EEV= 0.0D0 VV= 0.0D0 DO 180 I= 1,N EEV= EEV+EE(I)*V(I) 180 VV = VV +V(I)*V(I) OMEG =EEV/VV * DO 190 I= 1, N X(I)= X(I)+ALFA*P(I)+OMEG*EE(I) 190 RR(I)=EE(I)-OMEG*V(I) * CR=0.0D0 ERR = 0.0D0 DO 200 I = 1, N CR = CR+ R0(I)*RR(I) 200 ERR = ERR+RR(I)*RR(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 **************************************** KOKOMADE ******************* BETA = CR/(OMEG*SIG) DO 300 I = 1, N 300 P(I)= RR(I)+BETA*(P(I)-OMEG*W(I)) * 1000 CONTINUE 2000 KCOUNT = K RETURN END CCC* SUBROUTINE ILUCM2( N,M ,AA,AB,AC,AE,UB,DI,U,AB1,AC1,AE1,UB1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N),UB(-M:N), + AB1(-M:N+M), AC1(-M:N+M),AE1(-M:N+M),UB1(-M:N+M) DO 100 I = 1, N * AE1(I)= -AC1(I)*UB(I-M)*DI(I-M) UB1(I)= AB1(I)-AC1(I)*AE(I-M)*DI(I-M) * DIV=AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -AE1(I)*AE(I-M+1)*DI(I-M+1) + -(UB1(I)*AE(I-1)*DI(I-1)+AE1(I)*UB(I-M+1)*DI(I-M+1))*U DI(I)=1.0D0/DIV * UB(I)= AB(I)-AE1(I)*AC(I-M+1)*DI(I-M+1) AE(I)= -UB1(I)*AC(I-1)*DI(I-1) * 100 CONTINUE RETURN END *CCC SUBROUTINE BCG12(AA,AB,AC,AE,UB,DI,AB1,AC1,AE1,UB1,P,R,W,B0 + ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM,RRA,PA,RR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),R(N),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1) * DO 10 I = 1, N 10 R(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB1(I)*X(I-1)+AC1(I)*X(I-M1)) *** DO 30 I = 1, N 30 W(I)=(R(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) RR1 = 0.0D0 DO 60 I = 1, N P(I) = W(I) PA(I)= W(I) RR(I)= W(I) RRA(I)= W(I) 60 RR1 = RR1+RR(I)*RRA(I) *** ITERATION *** DO 1000 K = 1, 200 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB1(I)*P(I-1)+AC1(I)*P(I-M1) 100 CONTINUE DO 110 I = 1, N 110 W(I)=(AP(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 120 I = N, 1, -1 120 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) RAP= 0.0D0 DO 130 I=1, N 130 RAP = RAP+W(I)*PA(I) ALPHA = RR1/RAP ERR = 0.0D0 DO 200 I = 1, N X(I) = X(I)+ALPHA*P(I) RR(I) = RR(I)-ALPHA*W(I) 200 ERR = ERR+RR(I)*RR(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 *** DO 500 I = 1, N 500 W(I)=(PA(I)-AC(I-M1)*W(I-M1)-UB(I-1)*W(I-1) + -AE(I-M1+1)*W(I-M1+1))*DI(I) DO 510 I = N, 1, -1 510 W(I)=W(I)-DI(I)*(UB1(I+1)*W(I+1)+AC1(I+M1)*W(I+M1) + +AE1(I+M1-1)*W(I+M1-1)) DO 520 I= 1, N QA =AA(I)*W(I)+AB1(I+1)*W(I+1)+AC1(I+M1)*W(I+M1) + +AB(I-1)*W(I-1)+AC(I-M1)*W(I-M1) RRA(I)= RRA(I)-ALPHA*QA 520 CONTINUE RR2=RR1 RR1 = 0.0D0 DO 300 I = 1, N 300 RR1 = RR1+RR(I)*RRA(I) BETA = RR1/RR2 DO 400 I = 1, N P(I)= RR(I)+BETA*P(I) 400 PA(I) = RRA(I)+BETA*PA(I) 1000 CONTINUE 2000 KCOUNT = K RETURN END * SUBROUTINE CLOCK0 DIMENSION IA(3) COMMON /CLO/ IB(3) C CALL ITIME(IA) DO 200 I=1,3 200 IB(I)=IA(I) CC WRITE(6,*) (IB(I),I=1,3) RETURN END CCC SUBROUTINE CLOCK(CPUT) DIMENSION IA(3) COMMON /CLO/ IB(3) C CALL ITIME(IA) IJI= IA(1)-IB(1) IFUN=IA(2)-IB(2) IBYO=IA(3)-IB(3) CPUT=IJI*3600.0 +IFUN*60.0 +IBYO*1.0 CC WRITE(6,*) (IB(I),I=1,3),(IA(I),I=1,3),CPUT C CALL ITIME(IA) DO 200 I=1,3 200 IB(I)=IA(I) RETURN END