C CHUONG TRINH FAN TICH NHAN TO $LARGE $DEBUG ! Sao y tu ban Tieng Nga INCLUDE 'FTNT.INC' real*8 B(MMax),D(MMax),S(MMax),XBAR(MMax) real*8 T(MMax),V(MMax*MMax),R(MMax*MMax),TV(MMax*MMax) real*8 xx(NMax,MMax), x(NMax*MMax), Y(NMax,MMax), TB(MMax) real*8 TriRieng(MMax),VTORieng(MMax,MMax) real*8 TaiTrong(MMax,MMax),NhanTo(NMax,MMax) !! dimension x(1) 1 FORMAT(//20X,'FAN TICH NHAN TO '//10X,'Dung luong mau = ', /I6/10X,'So bien = ',7X,I6) 2 FORMAT(10X,'Trung binh '/1X,15(F8.1,1X)) 3 FORMAT(10X,'Do lech chuan '/1x,15(F8.1,1X)) 4 FORMAT(10X,'Ma tran tuong quan ') 5 FORMAT(3X,I3,2X,15F8.3) 6 FORMAT(10X,'Cac gia tri rieng'/1X,10f12.3) 7 FORMAT(10X,'Ty so tich luy cua cac tri rieng'/1X,10F8.4) 8 FORMAT(/10X,'Cac vec to rieng') 9 FORMAT(2X,15F8.3) 10 FORMAT(/10X,'Ma tran cac nhan to (',I3,' Nhan to )') 11 FORMAT(2X,10F11.4) 12 FORMAT(/5X,'Vong lap ',5X,'Phuong sai ') 13 FORMAT(7X,I4,9X,E9.4) 14 FORMAT(/10X,'Ma tran cac nhan to da quay (',I3,' Nhan to)') 15 FORMAT(2X,10F11.4) 16 FORMAT(/15X,'Kiem tra tinh tuong thich'/10x,'Bien thu',6X, /'Dau',6X,'Cuoi',6X,'Hieu') 17 FORMAT(12X,I3,6X,3(E9.2 )) 19 FORMAT(/10X,'Chi co',I3,' nhan to duoc giu lai. /Khong thuc hien phep quay') WRITE(*,*)' BAT DAU C.T FACTOR' open(1,file='solieu.txt', status='old') open(3,file='ketqua.txt',status='unknown') ! 100 read(1,*) n,m,con read(1,*) n,m,con do i=1,n read(1,*) (xx(i,j),j=1,m) enddo k=0 do j=1,m do i=1,n k=k+1 x(k)=xx(i,j) print*,k,x(k) enddo enddo ! Thu voi M=1 ! M=2 WRITE(3,1)N,M IO=1 CALL CORRE(N,M,IO,x,xbar,s,V,R,d,b,T) ! Tinh Ma tran Do lech (Y) cua X do i=1,N do j=1,M Y(i,j) = (XX(i,j)-XBAR(j))/SQRT(REAL(N)) enddo enddo WRITE(3,2)(XBAR(J),J=1,M) WRITE(3,3)(S(J),J=1,M) WRITE(3,4) DO 120 I=1,M DO 110 J=1,M IF(I-J) 102,104,104 102 L=I+(J*J-J)/2 GOTO 110 104 L=J+(I*I-I)/2 110 D(J)=R(L) 120 WRITE(3,5)I,(D(J),J=1,M) MV=0 CALL EIGEN(R,V,M,MV) ! Luu cac tri rieng va Vecto rieng do I=1,M L=I+(I*I-I)/2 TriRieng(I) = R(L) enddo L=0 do J=1,M do I=1,M L=L+1 VTORieng(i,j) = V(L) enddo enddo WRITE(*,*)' TRI RIENG 1,2,3,4,5=',R(1),R(3),R(6),R(10),R(15) WRITE(*,'(A\)')' VAO LAI THAM SO (CON)=' READ(*,*)CON CALL TRACE(M,R,CON,K,D) DO 130 I=1,K L=I+(I*I-I)/2 c gtr(i)=r(l) 130 S(I)=R(L) WRITE(3,6)(S(J),J=1,K) WRITE(3,7)(D(J),J=1,K) WRITE(3,8) L=0 DO 150 J=1,K DO 140 I=1,M L=L+1 140 D(I)=V(L) write(3,'(5x,'' Vec to '',i4)') j 150 WRITE(3,9) (D(I),I=1,M) CALL LOAD(M,K,R,V) WRITE(3,10)K DO 180 I=1,M DO 170 J=1,K L=M*(J-1)+I c a1(i,j)=v(l) 170 D(J)=V(L) write(3,'(5x,'' Bien thu '',i4)') i ! Luu Ma tran Tai trong Nhan to (K nhan to) do j=1,K TaiTrong(I,J) = D(j) enddo 180 WRITE(3,11) (D(J),J=1,K) IF(K-1)185,185,188 185 WRITE(3,19)K GOTO 100 188 CALL VARMX(M,K,V,NC,TV,B,T,D) NV=NC+1 WRITE(3,12) DO 190 I=1,NV NC=I-1 190 WRITE(3,13)NC,TV(I) WRITE(3,14)K DO 220 I=1,M DO 210 J=1,K L=M*(J-1)+I c a2(i,j)=v(l) 210 S(J)=V(L) write(3,'(5x,'' Bien thu '',i4)') i 220 WRITE(3,15) (S(J),J=1,K) WRITE(3,16) DO 230 I=1,M 230 WRITE(3,17)I,B(I),T(I),D(I) ! Tinh Ma tran Nhan to (F) do i=1,N do j=1,K NhanTo(i,j)=0.0 do L=1,M NhanTo(i,j)=NhanTo(i,j)+Y(i,L)*TaiTrong(L,j) enddo NhanTo(i,j)=NhanTo(i,j)/SQRT(TriRieng(j)) enddo enddo ! In Ket qua write(3,*)' Tri rieng' write(3,'(15F10.4)') (TriRieng(i),i=1,M) write(3,*)' Vecto rieng' do j=1,M write(3,'(i3,15F10.4)') j,(VTORieng(i,j),i=1,M) enddo write(3,*)' Tai trong Nhan to' do i=1,M write(3,'(i3,15F10.4)') i, (TaiTrong(i,j),j=1,K) enddo write(3,*)' Ma tran Nhan to' do i=1,N write(3,'(i3,15F10.4)') i,(NhanTo(i,j),j=1,M) enddo write(3,*)' Ma tran Do lech (Y) ban dau' do i=1,N do j=1,M ! Y(i,j)=Y(i,j) +XBAR(J) enddo write(3,'(i3,15F10.4)') i,(Y(i,j),j=1,M) enddo write(3,*)' Ma tran Do lech (Y) Uoc luong' do i=1,N do j=1,M Y(i,j) = 0.0 do L=1,K Y(i,j)=Y(i,j)+NhanTo(i,L)*TaiTrong(j,L) enddo enddo do j=1,M ! Y(i,j)=Y(i,j)+XBAR(J) enddo write(3,'(i3,15F10.4)') i,(Y(i,j),j=1,M) enddo write(3,*)' Ma tran Do lech (X-Y)' do i=1,N do j=1,M Y(i,j)=Y(i,j)-XX(i,j) enddo ! write(3,'(i3,15F10.4)') i,(Y(i,j),j=1,M) enddo 100 STOP END C*********************************************** c To hop cac c.t eigen,load,trace,varmx,corre,data c ************************************************ C CHUONG TRINH TIM K TRI RIENG SUBROUTINE TRACE(M,R,CON,K,D) INCLUDE 'FTNT.INC' REAL*8 R(MMAX*MMAX), D(MMAX) WRITE(*,*)' BAT DAU C.T TRACE' FM=M L=0 DO 100 I=1,M L=L+I 100 D(I)=R(L) K=0 DO 110 I=1,M IF(D(I)-CON) 120,105,105 105 K=K+1 110 D(I)=D(I)/FM 120 DO 130 I=2,K 130 D(I)=D(I)+D(I-1) RETURN END c C CHUONG TRINH TINH TAI TRONG NHAN TO SUBROUTINE LOAD(M,K,R,V) INCLUDE 'FTNT.INC' REAL*8 R(MMAX*MMAX),V(MMAX*MMAX),SQ WRITE(*,*)' BAT DAU C.T LOAD' L=0 JJ=0 DO 160 J=1,K JJ=JJ+J 150 SQ=SQRT(R(JJ)) DO 160 I=1,M L=L+1 160 V(L)=SQ*V(L) RETURN END C CHUONG TRINH TINH TRI RIENG VA VEC TO RIENG C SUBROUTINE EIGEN(A,R,N,MV) INCLUDE 'FTNT.INC' REAL*8 A(MMAX*MMAX),R(MMAX*MMAX) REAL*8 ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,COSX2,SINCS,RANGE 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 c WRITE(*,*)' 2' 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) c WRITE(*,*) ' 3' IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 c WRITE(*,*)' 4' 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 c WRITE(*,*)' 5' 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 C 6 c WRITE(*,*)' G/DOAN 6,7' 130 IF(M-N) 135,140,135 135 M=M+1 GOTO 60 c WRITE(*,*)' 8' 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 c WRITE(*,*)' 9' 160 IF(THR-ANRMX) 165,165,45 c WRITE(*,*)' 10' 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 c c C CHUONG TRINH QUAY NHAN TO SUBROUTINE VARMX(M,K,A,NC,TV,H,F,D) INCLUDE 'FTNT.INC' REAL*8 A(MMAX*MMAX),TV(NMAX),H(MMAX),F(MMAX),D(MMAX) 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 c c SUBROUTINE CORRE(N,M,I0,X,XBAR,STD,RX,R,B,D,T) INCLUDE 'FTNT.INC' REAL*8 X(NMAX*MMAX),XBAR(MMAX),STD(MMAX),RX(MMAX*MMAX) REAL*8 R(MMAX*MMAX),B(MMAX),D(MMAX),T(MMAX) 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 180I=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 C 225 R(JK)=R(JK)/(STD(J)*STD(K)) C WRITE(3,500) R(JK) C500 FORMAT(5X,F10.4) 230 CONTINUE C 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) C WRITE(*,*)' KET THUC CORRE' RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE DATA(M,D) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE FACTOR (X,F,A,N,M,K) C C...TINH F(N,K) TU PHUOG TRINH X(N,M) = F(N,K)*AT(K,M) C...VOI AT LA CHUYEN VI CUA A(M,K) C INCLUDE 'FTNT.INC' INTEGER N,M,K REAL*8 X(NMAX,MMAX), F(NMAX,MMAX), A(MMAX,MMAX), XA(NMAX,MMAX) REAL*8 AT(MMAX,MMAX), A1(MMAX,MMAX), AA(MMAX*MMAX),A2(MMAX,MMAX) C...CHUYEN VI DO I = 1,M DO J = 1,K AT(J,I) = A(I,J) ENDDO ENDDO C...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 C...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 C...NGHICH DAO MA TRAN AA CALL MINV(AA,K,D) C...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 C...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 C ! 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 !! C...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 C 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 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE MINV(A,N,D) ! ,L,M) INCLUDE 'FTNT.INC' REAL*8 A(MMAX*MMAX),L(MMAX*MMAX),M(MMAX*MMAX) 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