!============================================================== ! CHUONG TRINH TIM K TRI RIENG SUBROUTINE TRACE(M,R,CON,K,D) INCLUDE 'FTNT.INC' REAL :: R(MMAX*MMAX), D(MMAX) Integer :: M, K, I, L ! WRITE(*,*)' BAT DAU C.T TRACE' FM=M L=0 DO I=1,M L=L+I D(I)=R(L) ENDDO K=0 DO I=1,M IF (D(I) < CON) Exit K=K+1 D(I)=D(I)/FM ENDDO DO I=2,K D(I)=D(I)+D(I-1) ENDDO END SUBROUTINE TRACE ! !----------------------------------------------------------------------- ! CHUONG TRINH TINH TAI TRONG NHAN TO SUBROUTINE LOAD(M,K,R,V) INCLUDE 'FTNT.INC' REAL :: R(MMAX*MMAX),V(MMAX*MMAX),SQ Integer :: M, K, I, J, L ! WRITE(*,*)' BAT DAU C.T LOAD' L=0 JJ=0 DO J=1,K JJ=JJ+J SQ=SQRT(R(JJ)) DO I=1,M L=L+1 V(L)=SQ*V(L) ENDDO ENDDO END SUBROUTINE LOAD !----------------------------------------------------------------------- ! CHUONG TRINH TINH TRI RIENG VA VEC TO RIENG ! SUBROUTINE EIGEN(A,R,N,MV) INCLUDE 'FTNT.INC' REAL :: A(MMAX*MMAX),R(MMAX*MMAX) REAL :: ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,COSX2,SINCS,RANGE Integer :: N, I, J, M ! WRITE(*,*)' BAT DAU EIGEN' 5 RANGE=1.E-6 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0. IF(I-J) 20,15,20 15 R(IJ)=1. 20 CONTINUE 25 ANORM=0. DO 35 I=1,N DO 35 J=1,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM)165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(ABS(A(LM))-THR)130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=.5*(A(LL)-A(MM)) 68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/SQRT(2.*(1.+(SQRT(1.-Y*Y)))) SINX2=SINX*SINX 78 COSX=SQRT(1.-SINX2) COSX2=COSX*COSX SINCS=SINX*COSX ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GOTO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GOTO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X 130 IF(M-N) 135,140,135 135 M=M+1 GOTO 60 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GOTO 55 150 IF(IND-1) 160,155,160 155 IND=0 GOTO 50 160 IF(THR-ANRMX) 165,165,45 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE WRITE(*,*)' KET THUC EIGEN' RETURN END SUBROUTINE EIGEN ! !----------------------------------------------------------------------- ! CHUONG TRINH QUAY NHAN TO SUBROUTINE VARMX(M,K,A,NC,TV,H,F,D) INCLUDE 'FTNT.INC' REAL :: A(MMAX*MMAX),TV(NMAX),H(MMAX),F(MMAX),D(MMAX) Integer :: M, K, I, J WRITE(*,*)' BAT DAU C.T VARMX' EPS=0.0000116 TVLT=0. LL=K-1 NV=1 NC=0 FN=M FFN=FN*FN CONS=.7071066 DO 110 I=1,M H(I)=0 DO 110 J=1,K L=M*(J-1)+I 110 H(I)=H(I)+A(L)*A(L) DO 120 I=1,M 115 H(I)=SQRT(H(I)) DO 120 J=1,K L=M*(J-1)+I 120 A(L)=A(L)/H(I) GOTO 132 130 NV=NV+1 TVLT=TV(NV-1) 132 TV(NV)=0. DO 150 J=1,K AA=0. BB=0. LB=M*(J-1) DO 140 I=1,M L=LB+I CC=A(L)*A(L) AA=AA+CC 140 BB=BB+CC*CC 150 TV(NV)=TV(NV)+(FN*BB-AA*AA)/FFN IF(NV-51) 160 ,430,430 160 IF((TV(NV)-TVLT)-(1.E-7)) 170,170,190 170 NC=NC+1 IF(NC-3)190,190 ,430 190 DO 420 J=1,LL L1=M*(J-1) II=J+1 DO 420 K1=II,K L2=M*(K1-1) AA=0 BB=0 CC=0 DD=0 DO 230 I=1,M L3=L1+I L4=L2+I U=(A(L3)+A(L4))*(A(L3)-A(L4)) T=A(L3)*A(L4) T=T+T CC=CC+(U+T)*(U-T) DD=DD+2*U*T AA=AA+U 230 BB=BB+T T=DD-2*AA*BB/FN B=CC-(AA*AA-BB*BB)/FN IF(T-B) 280,240,320 240 IF((T+B)-EPS)420,250,250 250 COS4T=CONS SIN4T=CONS GOTO 350 280 TAN4T=ABS(T)/ABS(B) IF(TAN4T-EPS)300,290,290 290 COS4T=1/SQRT(1+TAN4T*TAN4T) SIN4T=TAN4T*COS4T GOTO 350 300 IF(B) 310,420,420 310 SINP=CONS COSP=CONS GOTO 400 320 CTN4T=ABS(T/B) IF(CTN4T-EPS) 340,330,330 330 SIN4T=1/SQRT(1+CTN4T*CTN4T) COS4T=CTN4T*SIN4T GOTO 350 340 COS4T=0. SIN4T=1. 350 COS2T=SQRT((1+COS4T)/2) SIN2T=SIN4T/(2*COS2T) 355 COST=SQRT((1+COS2T)/2) SINT=SIN2T/(2*COST) IF(B) 370,370,360 360 COSP=COST SINP=SINT GOTO 380 370 COSP=CONS*COST+CONS*SINT 375 SINP=ABS(CONS*COST-CONS*SINT) 380 IF(T) 390,390,400 390 SINP=-SINP 400 DO 410 I=1,M L3=L1+I L4=L2+I AA=A(L3)*COSP+A(L4)*SINP A(L4)=-A(L3)*SINP+A(L4)*COSP 410 A(L3)=AA 420 CONTINUE GOTO 130 430 DO 440 I=1,M DO 440 J=1,K L=M*(J-1)+I 440 A(L)=A(L)*H(I) NC=NV-1 DO 450 I=1,M 450 H(I)=H(I)*H(I) DO 470 I=1,M F(I)=0. DO 460 J=1,K L=M*(J-1)+I 460 F(I)=F(I)+A(L)*A(L) 470 D(I)=H(I)-F(I) RETURN END SUBROUTINE VARMX ! !----------------------------------------------------------------------- ! SUBROUTINE CORRE(N,M,I0,X,XBAR,STD,RX,R,B,D,T) INCLUDE 'FTNT.INC' REAL :: X(NMAX*MMAX),XBAR(MMAX),STD(MMAX),RX(MMAX*MMAX) REAL :: R(MMAX*MMAX),B(MMAX),D(MMAX),T(MMAX) Integer :: N, M, L, J, K, IO, I IO = 1 DO 100 J=1,M B(J)=0. 100 T(J)=0. K=(M*M+M)/2 DO 102 I=1,K 102 R(I)=0. FN=N L=0 IF(I0)105,127,105 105 DO 108 J=1,M DO 107 I=1,N L=L+1 107 T(J)=T(J)+X(L) XBAR(J)=T(J) 108 T(J)=T(J)/FN DO 115 I=1,N JK=0 L=I-N DO 110 J=1,M L=L+N D(J)=X(L)-T(J) 110 B(J)=B(J)+D(J) DO 115 J=1,M DO 115 K=1,J JK=JK+1 115 R(JK)=R(JK)+D(J)*D(K) GOTO 205 127 IF(N-M)130,130,135 130 KK=N GOTO 137 135 KK=M 137 DO 140 I=1,KK ! CALL DATA(M,D) DO 140 J=1,M T(J)=T(J)+D(J) L=L+1 140 RX(L)=D(J) FRR=KK DO 150 J=1,M XBAR(J)=T(J) 150 T(J)=T(J)/FRR L=0 DO 180 I=1,KK JK=0 DO 170 J=1,M L=L+1 170 D(J)=RX(L)-T(J) DO 180 J=1,M B(J)=B(J)+D(J) DO 180 K=1,J JK=JK+1 180 R(JK)=R(JK)+D(J)*D(K) IF(N-KK) 205,205,185 185 KK=N-KK DO 200 I=1,KK JK=0 ! CALL DATA(M,D) DO 190 J=1,M XBAR(J)=XBAR(J)+D(J) D(J)=D(J)-T(J) 190 B(J)=B(J)+D(J) DO 200 J=1,M DO 200 K=1,M JK=JK+1 200 R(JK)=R(JK)+D(J)*D(K) 205 JK=0 DO 210 J=1,M XBAR(J)=XBAR(J)/FN DO 210 K=1,J JK=JK+1 210 R(JK)=R(JK)-B(J)*B(K)/FN JK=0 DO 220 J=1,M JK=JK+J 220 STD(J)=SQRT(ABS(R(JK))) DO 230 J=1,M DO 230 K=J,M JK=J+(K*K-K)/2 L=M*(J-1)+K RX(L)=R(JK) L=M*(K-1)+J RX(L)=R(JK) IF(STD(J)-STD(K)) 225,222,225 222 R(JK)=1.0 GOTO 230 225 STD2=STD(J)*STD(K) ! PRINT*,STD2 R(JK)=R(JK)/STD2 ! 225 R(JK)=R(JK)/(STD(J)*STD(K)) ! WRITE(3,500) R(JK) ! 500 FORMAT(5X,F10.4) 230 CONTINUE ! WRITE(*,*)' Dung luong mau=',fn FN=SQRT(FN-1.) DO 240 J=1,M 240 STD(J)=STD(J)/FN L=-M DO 250 I=1,M L=L+M+1 250 B(I)=RX(L) ! WRITE(*,*)' KET THUC CORRE' RETURN END SUBROUTINE CORRE !----------------------------------------------------------------------- ! SUBROUTINE DATA(M,D) RETURN END SUBROUTINE DATA !----------------------------------------------------------------------- ! SUBROUTINE FACTOR (X,F,A,N,M,K) ! !...TINH F(N,K) TU PHUOG TRINH X(N,M) = F(N,K)*AT(K,M) !...VOI AT LA CHUYEN VI CUA A(M,K) ! INCLUDE 'FTNT.INC' INTEGER N,M,K, J, I REAL :: X(NMAX,MMAX), F(NMAX,MMAX), A(MMAX,MMAX), XA(NMAX,MMAX) REAL :: AT(MMAX,MMAX), A1(MMAX,MMAX), AA(MMAX*MMAX),A2(MMAX,MMAX) !...CHUYEN VI DO I = 1,M DO J = 1,K AT(J,I) = A(I,J) ENDDO ENDDO !...TICH AT(K,M)*A(M,K) = A1(K,K) DO I = 1,K DO J = 1,K A1(I,J) = 0.0 DO L = 1,M A1(I,J) = A1(I,J) + AT(I,L)*A(L,J) ENDDO ENDDO ENDDO !...CHUYEN THANH MANG 1 CHIEU DE TINH NGHICH DAO L = 0 DO J = 1,K DO I = 1,K L = L + 1 AA(L) = A1(I,J) ENDDO ENDDO !...NGHICH DAO MA TRAN AA CALL MINV(AA,K,D) !...LUU AA MOT CHIEU VAO A1 HAI CHIEU L = 0 DO J = 1,K DO I = 1,K L = L + 1 A1(I,J) = AA(L) ENDDO ENDDO !...TICH X(N,M)*A(M,K)*A1(K,K) = F(N,K) = XA(N,K)*A1(K,K) ! DO I = 1,N ! DO J = 1,K ! XA(I,J) = 0.0 ! DO L = 1,M ! XA(I,J) = XA(I,J) + X(I,L)*A(L,J) ! ENDDO ! ENDDO ! ENDDO ! ! DO I = 1,N ! DO J = 1,K ! F(I,J) = 0.0 ! DO L = 1,M ! F(I,J) = F(I,J) + XA(I,L)*A1(L,J) ! ENDDO ! ENDDO ! ENDDO !! !...TICH X(N,M)*A(M,K)*A1(K,K) = F(N,K) = X(N,K)*A2(K,K) DO I = 1,M DO J = 1,K A2(I,J) = 0.0 DO L = 1,K A2(I,J) = A2(I,J) + A(I,L)*A1(L,J) ENDDO ENDDO ENDDO ! DO I = 1,N DO J = 1,K F(I,J) = 0.0 DO L = 1,M F(I,J) = F(I,J) + X(I,L)*A2(L,J) ENDDO ENDDO ENDDO RETURN END SUBROUTINE FACTOR !----------------------------------------------------------------------- SUBROUTINE MINV(A,N,D) ! ,L,M) INCLUDE 'FTNT.INC' REAL :: A(MMAX*MMAX),L(MMAX*MMAX),M(MMAX*MMAX) Integer :: N, k, I D=1. NK=-N DO 80 K=1,N NK=NK+N L(K)=K M(K)=K KK=NK+K BIGA=A(KK) DO 20 J=K,N IZ=N*(J-1) DO 20 I=K,N IJ=IZ+I 10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20 15 BIGA=A(IJ) L(K)=I M(K)=J 20 CONTINUE J=L(K) IF(J-K)35,35,25 25 KI=K-N DO 30 I=1,N KI=KI+N HOLD=-A(KI) JI=KI-K+J A(KI)=A(JI) 30 A(JI)=HOLD 35 I=M(K) IF(I-K)45,45,38 38 JP=N*(I-1) DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A(JK) A(JK)=A(JI) 40 A(JI)=HOLD 45 IF(BIGA)48,46,48 46 D=0. RETURN 48 DO 55 I=1,N IF(I-K)50,55,50 50 IK=NK+I A(IK)=A(IK)/(-BIGA) 55 CONTINUE DO 65 I=1,N IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N IF(I-K)60,65,60 60 IF(J-K)62,65,62 62 KJ=IJ-I+K A(IJ)=HOLD*A(KJ)+A(IJ) 65 CONTINUE KJ=K-N DO 75 J=1,N KJ=KJ+N IF(J-K)70,75,70 70 A(KJ)=A(KJ)/BIGA 75 CONTINUE D=D*BIGA A(KK)=1./BIGA 80 CONTINUE K=N 100 K=(K-1) IF(K)150,150,105 105 I=L(K) IF(I-K)120,120,108 108 JQ=N*(K-1) JR=N*(I-1) DO 110 J=1,N JK=JQ+J HOLD=A(JK) JI=JR+J A(JK)=-A(JI) 110 A(JI)=HOLD 120 J=M(K) IF(J-K)100,100,125 125 KI=K-N DO 130 I=1,N KI=KI+N HOLD=A(KI) JI=KI-K+J A(KI)=-A(JI) 130 A(JI)=HOLD GOTO 100 150 RETURN END SUBROUTINE MINV !---------------------------------- End of this file -------------------------------------