program svd
	character*50 FileName1,FileName2,FileName3
	FileName1='e:\data2\sst\sst1'
	FileName2='e:\data2\r160\r1601.dat'
	FileName3='e:\temp\svd.dat'
	call czqSVD(FileName1,FileName2,FileName3)
	end
	
	SUBROUTINE czqSVD(FileName1,FileName2,FileName3)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c              This program uses for SVD analysing                  c
c This  program is from zhangyizhi&chenyanshan
c Parameter interpetion
cccc  MR=MIN(M1,M2)
c     X:左场,M1个格点,序列长度为N
c     Y:右场,M2个格点,序列长度为N
!     RAB:典型场的相关系数
!     RXB,RYA:异性相关系数
!     RXA,RYB:同性相关系数
!     SIGMA:奇异值
!     UA,VA:为两变量的线性组合

C     UL:左奇异场,VR:右奇异场
c         for the first R( R<=min(NS,NZ)) diagonal elements(cekma)  c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      parameter(M1=286,N=46,M2=160,NUM=8)
      parameter(MR=min(m1,m2),KA=MAX(M1,M2)+1)
	character*50 FileName1,FileName2,FileName3
      dimension A(M1,M2),U(M1,M1),V(M2,M2),S(KA),E(KA),WORK(KA)
      dimension X(M1,N),Y(M2,N),UA(MR,N),VA(MR,N)
      dimension RAB(MR,MR),RXB(M1,MR),RYA(M2,MR)
      dimension RXA(M1,M2),RYB(M2,MR)
      dimension SIGMA(MR,MR),UL(M1,MR),VR(M2,MR)
      EPS=0.000001
c-----input data----------------------------------------
c-----input data---------------------------------------- 
      open(1,file=FileName1)
      read(1,*)((X(I,J),J=1,N),I=1,M1)
      close(1)

	open(1,file=FileName2)
      read(1,*)((Y(I,J),J=1,N),I=1,M2)
      close(1)
100   FORMAT(30F4.0)
101   FORMAT(30F6.1)
c--------------------------------------------------------
c--------------------------------------------------------
CCCC  A=X*Y X(M1,N),Y(M2,N),A(M1,M2)
!!!!!!!CYS子程序用来将原始资料进行标准化处理!!!!!!!
	CALL standard(X,Y,m1,m2,n)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	CALL RXY(A,X,Y,M1,N,M2)
	CALL UAV(A,M1,M2,U,V,L,EPS,KA,S,E,WORK)
	CALL LSR(A,M1,M2,U,V,UL,VR,SIGMA,MR)
	!A=UL*SIGMA*INVERSE(VR)
      CALL TC(UA,VA,UL,VR,M1,MR,M2,X,Y,N)
	CALL R(RAB,UA,VA,MR,MR,N)
      CALL R(RXB,X,VA,M1,MR,N)
      CALL R(RYA,Y,UA,M2,MR,N)
      CALL R(RXA,X,UA,M1,MR,N)
      CALL R(RYB,Y,VA,M2,MR,N)
	CALL OUTER(SIGMA,RAB,RXB,RYA,RXA,RYB,ul,
     * 	vr,M1,M2,MR,NUM,FileName3)
	END

      SUBROUTINE R(RAB,UA,VA,MR1,MR2,N)
      dimension RAB(MR1,MR2),UA(MR1,N),VA(MR2,N)
      DO 3 I=1,MR1
       DO 4 J=1,MR2
        SUM1=0.0;SUM2=0.0;SUM3=0.0
	  AVE1=0.0;AVE2=0.0
	  DO K=1,N
	   AVE1=AVE1+UA(I,K)/N
	   AVE2=AVE2+VA(J,K)/N
	  ENDDO
        DO 5 K=1,N
         SUM1=SUM1+(UA(I,K)-AVE1)*(VA(J,K)-AVE2)
         SUM2=SUM2+(UA(I,K)-AVE1)*(UA(I,K)-AVE1)
         SUM3=SUM3+(VA(J,K)-AVE2)*(VA(J,K)-AVE2)
5       CONTINUE
        RAB(I,J)=SUM1/SQRT(SUM2)/SQRT(SUM3)
4      CONTINUE
3     CONTINUE
      RETURN
      END

      SUBROUTINE OUTER(SIGMA,RAB,RXB,RYA,RXA,RYB,
     *	ul,vr,M1,M2,MR,NUM,FileName3)
	character*50 FileName3
      dimension SIGMA(MR,MR),RAB(MR,MR),RXB(M1,MR),RYA(M2,MR)
      dimension RXA(M1,MR),RYB(M2,MR),ul(m1,mr),vr(m2,mr)
      REAL SUM1,SUM2,PRE1,PRE2
      OPEN(1,FILE=FileName3)

      SUM1=0.0
      DO 10 I=1,MR
10    SUM1=SUM1+SIGMA(I,I)*SIGMA(I,I)
      write(1,100)
      SUM2=0.0
      PER2=0.0
      DO 20 I=1,NUM
       SUM2=SUM2+SIGMA(I,I)*SIGMA(I,I)
       PRE1=SIGMA(I,I)*SIGMA(I,I)/SUM1
       PRE2=PRE2+PRE1
       write(1,101) SIGMA(I,I)*SIGMA(I,I),SUM2,PRE1,PRE2       
20    CONTINUE
      write(1,*)
      write(1,102)
      DO 30 I=1,NUM
30    write(1,103)I,RAB(I,I)
      write(1,*)
      write(1,104)
104   FORMAT(4X,'  异性相关系数              ')
      write(1,105)
105   FORMAT(4X,'  X  对右特征场  ')
      write(1,106)((RXB(I,J),J=1,NUM),I=1,M1)
      write(1,107)
107   FORMAT(4X,'  Y  对左特征场 ')
      write(1,106)((RYA(I,J),J=1,NUM),I=1,M2)
      write(1,108)
108   FORMAT(4X,'  同性相关系数         ')
      write(1,109)
109   FORMAT(4X,'  X 对左特征场 ')
      write(1,106)((RXA(I,J),J=1,NUM),I=1,M1)
      write(1,110)
110   FORMAT(4X,'  Y 对右特征场  ')
      write(1,106)((RYB(I,J),J=1,NUM),I=1,M2)
      write(1,112)
      write(1,111)(I,(ul(i,j),j=1,num),i=1,M1)
      write(1,113)
      write(1,111)(I,(vr(i,j),j=1,num),i=1,m2)
111   format(i4,8f8.3)
112   format(1x,'左特征场')
113   format(1x,'右特征场')
CCC
106   FORMAT(8F8.3)
CCC
100   FORMAT(4X,'  特征值    ',3X,'累积特征值  ',
     $        3X,' 贡献百分率 ',3X,'累积贡献    ')
101   FORMAT(4X,F12.4,3X,F12.4,3X,F12.4,3X,F12.4)
102   FORMAT(4X,'  模态相关:')
103   FORMAT(1X,I3,F12.4)
      close(1)
      RETURN
      END

      SUBROUTINE TC(UA,VA,UL,VR,M1,MR,M2,X,Y,N)
      dimension UA(MR,N),VA(MR,N),UL(M1,MR),VR(M2,MR),X(M1,N),Y(M2,N)
      DO 2 I=1,MR
      DO 2 J=1,N
       UA(I,J)=0.0
       DO 1 L=1,M1
        UA(I,J)=UA(I,J)+UL(L,I)*X(L,J)
1     CONTINUE
2     CONTINUE
      DO 4 I=1,MR
      DO 4 J=1,N
       VA(I,J)=0.0
       DO 3 L=1,M2
        VA(I,J)=VA(I,J)+VR(L,I)*Y(L,J)
3     CONTINUE
4     CONTINUE
      RETURN
      END

      SUBROUTINE LSR(A,M1,M2,U,V,UL,VR,SIGMA,MR)
      dimension A(M1,M2),U(M1,M1),V(M2,M2),UL(M1,MR),VR(M2,MR)
      dimension SIGMA(MR,MR)
       DO 91 I=1,MR
       DO 91 J=1,MR
       SIGMA(I,J)=0.0
       IF(I .EQ. J) SIGMA(I,J)=A(I,J)
91     CONTINUE
      DO 92 I=1,M1
      DO 92 J=1,MR
92    UL(I,J)=U(I,J)
      DO 93 I=1,MR
      DO 93 J=1,M2
	IF(V(I,J)  .GT. 1000.) THEN
	 WRITE(*,*) I,J
	ENDIF
93    VR(J,I)=V(I,J)

      RETURN
      END


      SUBROUTINE RXY(A,X,Y,M1,N,M2)
      dimension A(M1,M2),X(M1,N),Y(M2,N)
      DO 50 I=1,M1
      DO 50 J=1,M2
       A(I,J)=0.0
	 DO 10 L=1,N
        A(I,J)=A(I,J)+X(I,L)*Y(J,L)
10	 CONTINUE
50	CONTINUE
	RETURN
	END
      
	SUBROUTINE MUL(A,B,M,N,K,C)
	dimension A(M,N),B(N,K),C(M,K)
C      DOUBLE PRECISION A,B,C
	DO 50 I=1,M
	DO 50 J=1,K
	 C(I,J)=0.0
	 DO 10 L=1,N
	  C(I,J)=C(I,J)+A(I,L)*B(L,J)
10	 CONTINUE
50	CONTINUE
	RETURN
	END

      
	SUBROUTINE UAV(A,M,N,U,V,L,EPS,KA,S,E,WORK)
	dimension A(M,N),U(M,N),V(N,N),S(KA),E(KA),WORK(KA)
C      DOUBLE PRECISION A,U,V,S,E,D,WORK
      REAL DD,F,G,CS,SN,SHH,SK,EK,B,C,SM,SM1,EM1
	IT=60
	K=N
	IF(M-1 .LT. N) K=M-1
	L=M
	IF(N-2 .LT. M) L=N-2
	IF(L .LT. 0) L=0
	LL=K
	IF(L .GT. K) LL=L
	IF(LL .GE. 1) THEN
	 DO 150 KK=1,LL
	  IF(KK .LE. K) THEN
	   D=0.0
         DO 10 I=KK,M
10       D=D+A(I,KK)*A(I,KK)
         S(KK)=SQRT(D)
	   IF(S(KK) .NE. 0.0) THEN
	    IF(A(KK,KK) .NE. 0.0) S(KK)=SIGN(S(KK),A(KK,KK))
	    DO 20 I=KK,M
20        A(I,KK)=A(I,KK)/S(KK)
          A(KK,KK)=1.0+A(KK,KK)
         ENDIF
	   S(KK)=-S(KK)
	  ENDIF
	  IF(N .GE. KK+1) THEN
	   DO 50 J=KK+1,N
	    IF((KK .LE. K) .AND. (S(KK) .NE. 0.0)) THEN
	     D=0.0
		 DO 30 I=KK,M
30         D=D+A(I,KK)*A(I,J)
           D=-D/A(KK,KK)
		 DO 40 I=KK,M
40         A(I,J)=A(I,J)+D*A(I,KK)
          ENDIF
	    E(J)=A(KK,J)
50       CONTINUE
	  ENDIF
	  IF(KK .LE. K) THEN
	   DO 60 I=KK,M
60       U(I,KK)=A(I,KK)
        ENDIF
	  IF(KK .LE. L) THEN
	   D=0.0
	   DO 70 I=KK+1,N
70       D=D+E(I)*E(I)
         E(KK)=SQRT(D)
	   IF(E(KK) .NE. 0.0) THEN
	    IF(E(KK+1) .NE. 0.0) E(KK)=SIGN(E(KK),E(K+1))
	    DO 80 I=KK+1,N 
80        E(I)=E(I)/E(KK)
          E(KK+1)=1.0+E(KK+1)
         ENDIF
	   E(KK)=-E(KK)
	   IF((KK+1 .LE. M) .AND. (E(KK) .NE. 0.0)) THEN
	    DO 90 I=KK+1,M
90        WORK(I)=0.0
	    DO 110 J=KK+1,N
	     DO 100 I=KK+1,M
100        WORK(I)=WORK(I)+E(J)*A(I,J)
110       CONTINUE
	    DO 130 J=KK+1,N
	     DO 120 I=KK+1,M
120        A(I,J)=A(I,J)-WORK(I)*E(J)/E(KK+1)
130       CONTINUE
	   ENDIF
	   DO 140 I=KK+1,N
140      V(I,KK)=E(I)
        ENDIF
150    CONTINUE
	ENDIF
	MM=N
	IF(M+1 .LT. N) MM=M+1
	IF(K .LT. N) S(K+1)=A(K+1,K+1)
	IF(M .LT. MM) S(MM)=0.0
	IF(L+1 .LT. MM) E(L+1)=A(L+1,MM)
	E(MM)=0.0
	NN=M
	IF(M .GT. N) NN=N
	IF(NN .GE. K+1) THEN
	 DO 190 J=K+1,NN
	  DO 180 I=1,M
180	  U(I,J)=0.0
	 U(J,J)=1.0
190	 CONTINUE
	ENDIF
	IF(K .GE. 1) THEN
	 DO 250 LL=1,K
	  KK=K-LL+1
	  IF(S(KK) .NE. 0.0) THEN
	   IF(NN .GE.KK+1) THEN
	    DO 220 J=KK+1,NN
	     D=0.0
	     DO 200 I=KK,M
200	     D=D+U(I,KK)*U(I,J)/U(KK,KK)
	     D=-D
	     DO 210 I=KK,M
210        U(I,J)=U(I,J)+D*U(I,KK)
220	    CONTINUE
	   ENDIF
	   DO 225 I=KK,M
225	   U(I,KK)=-U(I,KK)
	   U(KK,KK)=1.0+U(KK,KK)
	   IF(KK-1 .GE. 1) THEN
	    DO 230 I=1,KK-1
230	    U(I,KK)=0.0
	   ENDIF
	   ELSE
	    DO 240 I=1,M
240	    U(I,KK)=0.0
	    U(KK,KK)=1.0
	   ENDIF
250	  CONTINUE
	 ENDIF
	 DO 300 LL=1,N
	  KK=N-LL+1
	  IF((KK .LE. L) .AND. (E(KK) .NE. 0.0)) THEN
	   DO 280 J=KK+1,N
	    D=0.0
	    DO 260 I=KK+1,N
260	    D=D+V(I,KK)*V(I,J)/V(KK+1,KK)
	    D=-D
	    DO 270 I=KK+1,N
270       V(I,J)=V(I,J)+D*V(I,KK)
280	   CONTINUE
	  ENDIF
	  DO 290 I=1,N
290	  V(I,KK)=0.0
	  V(KK,KK)=1.0
300	 CONTINUE
	 DO 305 I=1,M
	 DO 305 J=1,N
305	 A(I,J)=0.0
	 M1=MM
	 IT=60
310	 IF(MM .EQ. 0) THEN
	  L=0
	  IF(M .GE. N) THEN
	  I=N
	 ELSE
	  I=M
	 ENDIF
	 DO 315 J=1,I-1
	  A(J,J)=S(J)
	  A(J,J+1)=E(J)
315	 CONTINUE
	 A(I,I)=S(I)
	 IF(M .LT. N) A(I,I+1)=E(I)
	 DO 314 I=1,N-1
	  DO 313 J=I+1,N
	   D=V(I,J)
	   V(I,J)=V(J,I)
	   V(J,I)=D
313     CONTINUE
314    CONTINUE
	 RETURN
	ENDIF
	IF(IT .EQ.0) THEN
	 L=MM
	 IF(M .GE. N) THEN
	  I=N
	  ELSE
	  I=M
	 ENDIF
	 DO 316 J=1,I-1
	  A(J,J)=S(J)
	  A(J,J+1)=E(J)
316    CONTINUE
	 A(I,I)=S(I)
	 IF(M .LT. N) A(I,I+1)=E(I)
	 DO 318 I=1,N-1
	  DO 317 J=I+1,N
	   D=V(I,J)
	   V(I,J)=V(J,I)
	   V(J,I)=D
317	  CONTINUE
318	 CONTINUE
	 RETURN
	ENDIF
	KK=MM	
320	KK=KK-1
	IF(KK .NE. 0) THEN
	 D=ABS(S(KK))+ABS(S(KK+1))
	 DD=ABS(E(KK))
	 IF(DD .GT. EPS*D) GOTO 320
	 E(KK)=0.0
	ENDIF
	IF(KK .EQ. MM-1) THEN
	 KK=KK+1
	 IF(S(KK) .LT. 0.0) THEN
	  S(KK)=-S(KK)
	  DO 330 I=1,N
330	  V(I,KK)=-V(I,KK)
	 ENDIF
335	 IF(KK .NE. M1) THEN
	  IF(S(KK) .LT. S(KK+1)) THEN
	   D=S(KK)
	   S(KK)=S(KK+1)
	   S(KK+1)=D
	   IF(KK .LT. N) THEN
	    DO 340 I=1,N
	     D=V(I,KK)
	     V(I,KK)=V(I,KK+1)
	     V(I,KK+1)=D
340	    CONTINUE
	   ENDIF
	   IF(KK .LT. N) THEN
	    DO 350 I=1,M
	     D=U(I,KK)
	     U(I,KK)=U(I,KK+1)
	     U(I,KK+1)=D
350	    CONTINUE
	   ENDIF
	   KK=KK+1
	   GOTO 335
	  ENDIF
	 ENDIF
	 IT=60
	 MM=MM-1
	 GOTO 310
	ENDIF
	KS=MM+1
360	KS=KS-1
	IF(KS .GT. KK) THEN
	 D=0.0
	 IF(KS .NE. MM) D=D+ABS(E(KS))
	 IF(KS .NE. KK+1) D=D+ABS(E(KS-1))
	 DD=ABS(S(KS))
	 IF(DD .GT. EPS*D) GOTO 360
	 S(KS)=0.0
	ENDIF
	IF(KS .EQ. KK) THEN
	 KK=KK+1
	 D=ABS(S(MM))
	 IF(ABS(S(MM-1)) .GT. D) D=ABS(S(MM-1))
	 IF(ABS(E(MM-1)) .GT. D) D=ABS(E(MM-1))
	 IF(ABS(S(KK)) .GT. D) D=ABS(S(KK))
	 IF(ABS(E(KK)) .GT. D) D=ABS(E(KK))
	 SM=S(MM)/D
	 SM1=S(MM-1)/D
	 EM1=E(MM-1)/D
	 SK=S(KK)/D
	 EK=E(KK)/D
	 B=((SM1+SM)*(SM1-SM)+EM1*EM1)/2.0
	 C=SM*EM1
	 C=C*C
	 SHH=0.0
	 IF((B .NE. 0.0) .OR. (C .NE. 0.0)) THEN
	  SHH=SQRT(B*B+C)
	  IF(B .LT. 0.0) SHH=-SHH
	  SHH=C/(B+SHH)
	 ENDIF
	 F=(SK+SM)*(SK-SM)-SHH
	 G=SK*EK
	 DO 400 I=KK,MM-1
	  CALL SSS(F,G,CS,SN)
	  IF(I .NE. KK) E(I-1)=F
	  F=CS*S(I)+SN*E(I)
	  E(I)=CS*E(I)-SN*S(I)
	  G=SN*S(I+1)
	  S(I+1)=CS*S(I+1)
	  IF((CS .NE. 1.0) .OR. (SN .NE. 0.0)) THEN
	   DO 370 J=1,N
	    D=CS*V(J,I)+SN*V(J,I+1)
	    V(J,I+1)=-SN*V(J,I)+CS*V(J,I+1)
	    V(J,I)=D
370	   CONTINUE
	  ENDIF
	  CALL SSS(F,G,CS,SN)
	  S(I)=F
	  F=CS*E(I)+SN*S(I+1)
	  S(I+1)=-SN*E(I)+CS*S(I+1)
	  G=SN*E(I+1)
	  E(I+1)=CS*E(I+1)
	  IF(I .LT. M) THEN
	   IF((CS .NE. 1.0) .OR. (SN .NE. 0.0)) THEN
	    DO 380 J=1,M
	     D=CS*U(J,I)+SN*U(J,I+1)
	     U(J,I+1)=-SN*U(J,I)+CS*U(J,I+1)
	     U(J,I)=D
380	    CONTINUE
	   ENDIF
	  ENDIF
400	 CONTINUE
	 E(MM-1)=F
	 IT=IT-1
	 GOTO 310
	ENDIF
	IF(KS .EQ. MM) THEN
	 KK=KK+1
	 F=E(MM-1)
	 E(MM-1)=0.0
	 DO 420 LL=KK,MM-1
	  I=MM+KK-LL-1
	  G=S(I)
	  CALL SSS(G,F,CS,SN)
	  S(I)=G
	  IF(I .NE. KK) THEN
	   F=-SN*E(I-1)
	   E(I-1)=CS*E(I-1)
	  ENDIF
	  IF((CS .NE. 1.0) .OR. (SN .NE. 0.0)) THEN
	   DO 410 J=1,N
	    D=CS*V(J,I)+SN*V(J,MM)
	    V(J,MM)=-SN*V(J,I)+CS*V(J,MM)
	    V(J,I)=D
410	   CONTINUE
	  ENDIF
420	 CONTINUE
	 GOTO 310
	ENDIF
	KK=KS+1
	F=E(KK-1)
	E(KK-1)=0.0
	DO 450 I=KK,MM
	 G=S(I)
	 CALL SSS(G,F,CS,SN)
	 S(I)=G
	 F=-SN*E(I)
	 E(I)=CS*E(I)
	 IF((CS .NE. 1.0) .OR. (SN .NE. 0.0)) THEN
	  DO 430 J=1,M
	   D=CS*U(J,I)+SN*U(J,KK-1)
	   U(J,KK-1)=-SN*U(J,I)+CS*U(J,KK-1)
	   U(J,I)=D
430	  CONTINUE
	 ENDIF
450	CONTINUE
	GOTO 310
	END

	SUBROUTINE SSS(F,G,CS,SN)
      REAL D,R
	IF((ABS(F)+ABS(G)) .EQ. 0.0) THEN
	 CS=1.0
	 SN=0.0
	 D=0.0
	ELSE
	 D=SQRT(F*F+G*G)
	 IF(ABS(F) .GT. ABS(G)) D=SIGN(D,F)
	 IF(ABS(G) .GE. ABS(F)) D=SIGN(D,G)
	 CS=F/D
	 SN=G/D
	ENDIF
	R=1.0
	IF(ABS(F) .GT. ABS(G)) THEN
	 R=SN
	ELSE
	 IF(CS.NE.0.0) R=1.0/CS
	ENDIF
	F=D
	G=R
	RETURN
	END
	
      SUBROUTINE standard(x,y,m1,m2,n)
	dimension x(m1,n),y(m2,n)
	dimension xm(m1),xd(m1),ym(m2),yd(m2)
      fnx=real(m1)
      fny=real(m2)
      fn=real(n)
      do 10 i=1,m1
      xm(i)=0.0
      do 12  k=1,n
 12   xm(i)=xm(i)+x(i,k)
      xm(i)=xm(i)/fn
      xd(i)=0.0
      do 13  k=1,n
 13   xd(i)=xd(i)+(x(i,k)-xm(I))**2
      xd(i)=sqrt(xd(i)/fn)
 10   continue
      do 20 i=1,m2
      ym(i)=0.0
      do 22  k=1,n
 22   ym(i)=ym(i)+y(i,k)
      ym(i)=ym(i)/fn
      yd(i)=0.0
      do 23  k=1,n
 23   yd(i)=yd(i)+(y(i,k)-ym(I))**2
      yd(i)=sqrt(yd(i)/fn)
 20   continue
      do  30  i=1,m1
      do  30  k=1,n
 30   x(i,k)=(x(i,k)-xm(I))/xd(I)
      do  32  i=1,m2
      do  32  k=1,n
 32   y(i,k)=(y(i,k)-ym(I))/yd(I)
	RETURN
	END