svd.for

svd气象中用的fortran程序

Clarge!
      PROGRAM SVD
      PARAMETER(N=686,IP=160,IT=50)
      DIMENSION X(N,IP),U(N,N),V(IP,IP),E(N),WORK(N),S(N)
      dimension TC1(N,IT),TC2(IP,IT),R(16),AM1(N),AM2(IP)
      dimension SSTA(IP,IT),SLP(N,IT)
    
            open(3,file='d:\TANG\data\rj7.DAT',status='old')
            read(3,*) ((ssta(i,j),i=1,ip),j=1,it)
            open(1,file='d:\tang\data\qs1.dat',status='old')
            read(1,*) ((slp(i,j),i=1,n),j=1,it)
	      close(1)
            close(3)



	
      LDX=n
      LDU=n
      LDV=ip
      im=16
 
      CALL CACOR(SSTA,SLP,X,N,IP,IT)
 
      CALL SSVDC(X,LDX,N,IP,S,E,U,LDU,V,LDV,WORK,11,INFO)
 
c
      CALL TCFC(U,N,TC1,SLP,AM1,IT)
      CALL TCFC(V,IP,TC2,SSTA,AM2,IT)
      CALL TCHXG(R,TC1,TC2,iM,IP,N,IT)
      W=0.0
      DO 21 I=1,N
  21  W=W+S(I)*S(I)
      DO 22 I=1,10
  22  S(I)=S(I)*S(I)/W
      open(11,file='d:\tang\svdjg\1s7r.dat')
      write(11,*) 'u=left'
      write(11,101) ((u(i,j),j=1,5),i=1,n)
      write(11,*) 'v=right'
      write(11,101) ((v(i,j),j=1,5),i=1,ip)
      write(11,*) 'tc1=left'
      write(11,101) ((tc1(i,j),i=1,5),j=1,it)
      write(11,*) 'tc2=right'
      write(11,101) ((tc2(i,j),i=1,5),j=1,it)
      write(11,*) 's(i)='
      write(11,100) (s(i),i=1,10)
      write(11,*) 'am1='
      write(11,100) (am1(I),I=1,10)
      write(11,*) 'am2='
      write(11,100) (am2(I),I=1,10)
      write(11,*) 'r(i)='
      write(11,100) (r(I),I=1,10)
 101  format(5f10.3)
 100  format(10f7.3)
	close(11)
      STOP
      END
CCCCCC
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE SSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
      INTEGER LDX,N,P,LDU,LDV,JOB,INFO
      REAL X(LDX,LDV),S(LDU),E(LDU),U(LDU,LDU),V(LDV,LDU),WORK(LDU)
      INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
     *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
      REAL SDOT,T
      REAL B,C,CS,EL,EMM1,F,G,SNRM2,SCALE,SHIFT,SL,SM,SN,SMM1,T1,TEST,
     *     ZTEST
      LOGICAL WANTU,WANTV
      MAXIT=30
      WANTU=.FALSE.
      WANTV=.FALSE.
      JOBU=MOD(JOB,100)/10
      NCU=N
      IF(JOBU.GT.1) NCU=MIN0(N,P)
      IF(JOBU.NE.0) WANTU=.TRUE.
      IF(MOD(JOB,10).NE.0) WANTV=.TRUE.
      INFO=0
      NCT=MIN0(N-1,P)
      NRT=MAX0(0,MIN0(P-2,N))
      LU=MAX0(NCT,NRT)
      IF(LU.LT.1) GO TO 170
      DO 160 L=1,LU
      LP1=L+1
      IF(L.GT.NCT) GO TO 20
      S(L)=SNRM2(N-L+1,X(L,L),1)
      IF(S(L).EQ.0.0E0) GO TO 10
      IF(X(L,L).NE.0.0E0) S(L)=SIGN(S(L),X(L,L))
      CALL SSCAL(N-L+1,1.0E0/S(L),X(L,L),1)
      X(L,L)=1.0E0+X(L,L)
   10 CONTINUE
      S(L)=-S(L)
   20 CONTINUE
      IF(P.LT.LP1) GO TO 50
      DO 40 J=LP1,P
      IF(L.GT.NCT) GO TO 30
      IF(S(L).EQ.0.0E0) GO TO 30
      T=-SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
      CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
   30 CONTINUE
      E(J)=X(L,J)
   40 CONTINUE
   50 CONTINUE
      IF(.NOT.WANTU .OR. L .GT.NCT) GO TO 70
      DO 60 I=L,N
      U(I,L)=X(I,L)
   60 CONTINUE
   70 CONTINUE
      IF(L.GT.NRT) GO TO 150
      E(L)=SNRM2(P-L,E(LP1),1)
      IF(E(L).EQ.0.0E0) GO TO 80
      IF(E(LP1).NE.0.0E0) E(L)=SIGN(E(L),E(LP1))
      CALL SSCAL(P-L,1.0E0/E(L),E(LP1),1)
      E(LP1)=1.0E0+E(LP1)
   80 CONTINUE
      E(L)=-E(L)
      IF(LP1.GT.N.OR.E(L).EQ.0.0E0) GO TO 120
      DO 90 I=LP1,N
      WORK(I)=0.0E0
   90 CONTINUE
      DO 100 J=LP1,P
      CALL SAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
  100 CONTINUE
      DO 110 J=LP1,P
      CALL SAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
  110 CONTINUE
  120 CONTINUE
      IF(.NOT.WANTV) GO TO 140
      DO 130 I=LP1,P
      V(I,L)=E(I)
  130 CONTINUE
  140 CONTINUE
  150 CONTINUE
  160 CONTINUE
  170 CONTINUE
      M=MIN0(P,N+1)
      NCTP1=NCT+1
      NRTP1=NRT+1
      IF(NCT.LT.P) S(NCTP1)=X(NCTP1,NCTP1)
      IF(N.LT.M) S(M)=0.0E0
      IF(NRTP1.LT.M) E(NRTP1)=X(NRTP1,M)
      E(M)=0.0E0
      IF(.NOT.WANTU) GO TO 300
      IF(NCU.LT.NCTP1) GO TO 200
      DO 190 J=NCTP1,NCU
      DO 180 I=1,N
      U(I,J)=0.0E0
  180 CONTINUE
      U(J,J)=1.0E0
  190 CONTINUE
  200 CONTINUE
      IF(NCT.LT.1) GO TO 290
      DO 280 LL=1,NCT
      L=NCT-LL+1
      IF(S(L).EQ.0.0E0) GO TO 250
      LP1=L+1
      IF(NCU.LT.LP1) GO TO 220
      DO 210 J=LP1,NCU
      T=-SDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
      CALL SAXPY(N-L+1,T,U(L,L),1,U(L,J),1)
  210 CONTINUE
  220 CONTINUE
      CALL SSCAL(N-L+1,-1.0E0,U(L,L),1)
      U(L,L)=1.0E0+U(L,L)
      LM1=L-1
      IF(LM1.LT.1) GO TO 240
      DO 230 I=1,LM1
      U(I,L)=0.0E0
  230 CONTINUE
  240 CONTINUE
      GO TO 270
  250 CONTINUE
      DO 260 I=1,N
      U(I,L)=0.0E0
  260 CONTINUE
      U(L,L)=1.0E0
  270 CONTINUE
  280 CONTINUE
  290 CONTINUE
  300 CONTINUE
      IF(.NOT.WANTV) GO TO 350
      DO 340 LL=1,P
      L=P-LL+1
      LP1=L+1
      IF(L.GT.NRT) GO TO 320
      IF(E(L).EQ.0.0E0) GO TO 320
      DO 310 J=LP1,P
      T=-SDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
      CALL SAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)
  310 CONTINUE
  320 CONTINUE
      DO 330 I=1,P
      V(I,L)=0.0E0
  330 CONTINUE
      V(L,L)=1.0E0
  340 CONTINUE
  350 CONTINUE
      MM=M
      ITER=0
  360 CONTINUE
      IF(M.EQ.0)  GO TO 620
      IF(ITER.LT.MAXIT) GO TO 370
      INFO=M
      GO TO 620
  370 CONTINUE
      DO 390 LL=1,M
      L=M-LL
      IF(L.EQ.0) GO TO 400
      TEST=ABS(S(L))+ABS(S(L+1))
      ZTEST=TEST+ABS(E(L))
      IF(ZTEST.NE.TEST) GO TO 380
      E(L)=0.0E0
      GO TO 400
  380 CONTINUE
  390 CONTINUE
  400 CONTINUE
      IF(L.NE.M-1) GO TO 410
      KASE=4
      GO TO 480
  410 CONTINUE
      LP1=L+1
      MP1=M+1
      DO 430 LLS=LP1,MP1
      LS=M-LLS+LP1
      IF(LS.EQ.L) GO TO 440
      TEST=0.0E0
      IF(LS.NE.M) TEST=TEST+ABS(E(LS))
      IF(LS.NE.L+1) TEST=TEST+ABS(E(LS-1))
      ZTEST=TEST+ABS(S(LS))
      IF(ZTEST.NE.TEST) GO TO 420
      S(LS)=0.0E0
      GO TO 440
  420 CONTINUE
  430 CONTINUE
  440 CONTINUE
      IF(LS.NE.L) GO TO 450
      KASE=3
      GO TO 470
  450 CONTINUE
      IF(LS.NE.M) GO TO 460
      KASE=1
      GO TO 470
  460 CONTINUE
      KASE=2
      L=LS
  470 CONTINUE
  480 CONTINUE
      L=L+1
      GO TO (490,520,540,570), KASE
  490 CONTINUE
      MM1=M-1
      F=E(M-1)
      E(M-1)=0.0E0
      DO 510 KK=L,MM1
      K=MM1-KK+L
      T1=S(K)
      CALL SROTG(T1,F,CS,SN)
      S(K)=T1
      IF(K.EQ.L) GO TO 500
      F=-SN*E(K-1)
      E(K-1)=CS*E(K-1)
  500 CONTINUE
      IF(WANTV) CALL SROT(P,V(1,K),1,V(1,M),1,CS,SN)
  510 CONTINUE
      GO TO 610
  520 CONTINUE
      F=E(L-1)
      E(L-1)=0.0E0
      DO 530 K=L,M
      T1=S(K)
      CALL SROTG(T1,F,CS,SN)
      S(K)=T1
      F=-SN*E(K)
      E(K)=CS*E(K)
      IF(WANTU) CALL SROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
  530 CONTINUE
      GO TO 610
  540 CONTINUE
      SCALE=AMAX1(ABS(S(M)),ABS(S(M-1)),ABS(E(M-1)),ABS(S(L)),
     *            ABS(E(L)))
      SM=S(M)/SCALE
      SMM1=S(M-1)/SCALE
      EMM1=E(M-1)/SCALE
      SL=S(L)/SCALE
      EL=E(L)/SCALE
      B=((SMM1+SM)*(SMM1-SM)+EMM1**2)/2.0E0
      C=(SM*EMM1)**2
      SHIFT=0.0E0
      IF(B.EQ.0.0E0.AND.C.EQ.0.0E0) GO TO 550
      SHIFT=SQRT(B**2+C)
      IF(B.LT.0.0E0) SHIFT=-SHIFT
      SHIFT=C/(B+SHIFT)
  550 CONTINUE
      F=(SL+SM)*(SL-SM)-SHIFT
      G=SL*EL
      MM1=M-1
      DO 560 K=L,MM1
      CALL SROTG(F,G,CS,SN)
      IF(K.NE.L) E(K-1)=F
      F=CS*S(K)+SN*E(K)
      E(K)=CS*E(K)-SN*S(K)
      G=SN*S(K+1)
      S(K+1)=CS*S(K+1)
      IF(WANTV) CALL SROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
      CALL SROTG(F,G,CS,SN)
      S(K)=F
      F=CS*E(K)+SN*S(K+1)
      S(K+1)=-SN*E(K)+CS*S(K+1)
      G=SN*E(K+1)
      E(K+1)=CS*E(K+1)
      IF(WANTU.AND.K.LT.N)
     *CALL SROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
  560 CONTINUE
      E(M-1)=F
      ITER=ITER+1
      GO TO 610
  570 CONTINUE
      IF(S(L).GE.0.0E0) GO TO 580
      S(L)=-S(L)
      IF(WANTV) CALL SSCAL(P,-1.0E0,V(1,L),1)
  580 CONTINUE
  590 IF(L.EQ.MM) GO TO 600
      IF(S(L).GE.S(L+1)) GO TO 600
      T=S(L)
      S(L)=S(L+1)
      S(L+1)=T
      IF(WANTV.AND.L.LT.P)
     *CALL SSWAP(P,V(1,L),1,V(1,L+1),1)
      IF(WANTU.AND.L.LT.N)
     *CALL SSWAP(N,U(1,L),1,U(1,L+1),1)
      L=L+1
      GO TO 590
  600 CONTINUE
      ITER=0
      M=M-1
  610 CONTINUE
      GO TO 360
  620 CONTINUE
      RETURN
      END
CCCCC ISAMAX
      INTEGER FUNCTION ISAMAX(N,SX,INCX)
      REAL SX(1),SMAX
      INTEGER I,INCX,IX,N
      ISAMAX=0
      IF(N.LT.1) RETURN
      ISAMAX=1
      IF(N.EQ.1) RETURN
      IF(INCX.EQ.1) GO TO 20
      IX=1
      SMAX=ABS(SX(1))
      IX=IX+INCX
      DO 10 I=2,N
      IF(ABS(SX(IX)).LE.SMAX) GO TO 5
      ISAMAX=I
      SMAX=ABS(SX(IX))
   5  IX=IX+INCX
   10 CONTINUE
      RETURN
   20 SMAX=ABS(SX(1))
      DO 30 I=2,N
      IF(ABS(SX(I)).LE.SMAX) GO TO 30
      ISAMAX=I
      SMAX=ABS(SX(I))
   30 CONTINUE
      RETURN
      END
CCCCC SASUM
      REAL FUNCTION SASUM(N,SX,INCX)
      REAL SX(1),STEMP
      INTEGER I,INCX,M,MP1,N,NINCX
      SASUM=0.0E0
      STEMP=0.0E0
      IF(N.LE.0) RETURN
      IF(INCX.EQ.1) GO TO 20
      NINCX=N*INCX
      DO 10 I=1,NINCX,INCX
      STEMP=STEMP+ABS(SX(I))
   10 CONTINUE
      SASUM=STEMP
      RETURN
   20 M=MOD(N,6)
      IF(M.EQ.0) GO TO 40
      DO 30 I=1,M
      STEMP=STEMP+ABS(SX(I))
   30 CONTINUE
      IF(N.LT.6) GO TO 60
   40 MP1=M+1
      DO 50 I=MP1,N,6
      STEMP=STEMP+ABS(SX(I))+ABS(SX(I+1))+ABS(SX(I+2))