第四章主分量分析.docVIP

PAGE 80 PAGE 10 第四章主分量分析 4．２主分量分析（ＰＣ）程序 4.2.1功能 计算要素场的主分量。 4.2.2方法说明 4.2.3子程序语句 SUBROUTINE ZFL(X,N,P,V,B) 4.2.4哑元说明 X——二维实型数组，大小为，输入参量，存放P个空间点的N次观测值。 N——整型变量，输入参数，序列的长度。 P——整型变量，输入参数，空间点的个数。 V——二维实型数组，大小为，输出参数，其列向量为特征向量。 B——一维实型数组，大小为P，输出参数，存放特征值，从大到小排列。 4.2.5子程序 SUBROUTINE ZFL(X,N,P,V,B,G) INTEGER::P REAL(8),DIMENSION(P,N)::X REAL(8),DIMENSION(P)::XV REAL(8),DIMENSION(P)::B,G REAL(8),DIMENSION(P,P)::V REAL(8),DIMENSION(P,P)::S ! 求每一点的均值 XV=0 DO J=1,P DO I=1,N XV(J)=XV(J)+X(I,J) END DO XV(J)=XV(J)/N END DO ! 求协方差阵,该协方差阵为对称方阵 S=0 DO I=1,P DO J=1,P DO K=1,N S(I,J)=S(I,J)+(X(I,K)-XV(I))*(X(J,K)-XV(J)) END DO S(I,J)=S(I,J)/N END DO END DO CALL JCB(S,P,1.0E-8,V,B,L) ! 调用雅可比法求解矩阵的特征值和特征向量 CALL MSORT(B,P,V) ! 将特征值按大小排列 GP=0 DO I=1,P GP=GP+B(I) END DO DO I=1,P G(I)=B(I)/GP END DO END SUBROUTINE JCB(A,N,EPS,V,B,L) ! A:调用时存放实对称矩阵 ! B:返回时存放矩阵的全部特征值 ! V:存放特征向量,其中第i列为与第i个特征值相对应的特征向量 ! EPS:存放精度要求 REAL(8),DIMENSION(N,N)::A,V REAL(8),DIMENSION(N)::B REAL(8)::FM,CN,SN,OMEGA,X,Y INTEGER::P,Q L=1 V=0 DO I=1,N V(I,I)=1 END DO 10 FM=0 DO I=1,N DO J=1,N IF(I/=J.AND.ABS(A(I,J))FM)THEN FM=ABS(A(I,J)) P=I Q=J END IF END DO END DO IF(FMEPS)THEN L=1 RETURN END IF IF(L100)THEN L=0 RETURN END IF L=L+1 X=-A(P,Q) Y=(A(Q,Q)-A(P,P))/2 OMEGA=X/SQRT(X*X+Y*Y) IF(Y0)OMEGA=-OMEGA SN=1+SQRT(1-OMEGA*OMEGA) SN=OMEGA/SQRT(2*SN) CN=SQRT(1-SN*SN) FM=A(P,P) A(P,P)=FM*CN*CN+A(Q,Q)*SN*SN+A(P,Q)*OMEGA A(Q,Q)=FM*SN*SN+A(Q,Q)*CN*CN-A(P,Q)*OMEGA A(P,Q)=0 A(Q,P)=0 DO J=1,N IF(J/=P.AND.J/=Q)THEN FM=A(P,J) A(P,J)=FM*CN+A(Q,J)*SN A(Q,J)=-FM*SN+A(Q,J)*CN END IF END DO DO I=1,N IF(I/=P.AND.I/=Q)THEN FM=A(I,P) A(I,P)=FM*CN+A(I,Q)*SN A(I,Q)=-FM*SN+A(I,Q)*CN END IF END DO DO I=1,N FM=V(I,P) V(I,P)=FM*CN+V(I,Q)*SN V(I,Q)=-FM*SN+V(I,Q)*CN END DO DO I=1,N B(I)=A(I,I) END DO GOTO 10 END SUBROUTINE MSORT(B,N,V) ! 将特征值按大小排列 REAL(8),DIMENS