c  Fortran program to perform singular value decomposition on 
c  a space-time array with eigenvectors and PCs rotated according
c  to the varimax criterion.  Varimax subroutine originally from
c  EOF/PC code used at NASA GSFC.
c
c  **NOTE: DO NOT use singular values to compute explained variance
c  of EOFs if using varimax rotation!  Compute explained variance of
c  the rotated mode by calculating the variance of the raw PCs and dividing
c  by the sum of the variance of all raw PCs
c
c  Input file is the space-time array in binary format (AT array)
c
c  Output files are:
c     V array (spacedim x spacedim): Normalized vectors (var=1)
c     U array (timedim x timedim): Normalized vectors (var=1)
c     Singular values (sigma) (maxdim) 
c     Raw PCs (timedim x timedim): Unnormalized Rotated PCs
c
c

c  Specify space and time dimensions
c
c  maxdim is the larger of timedim or spacedim
c  mindim is the smaller of timedim or spacedim
c  rotations is the number of eigenvector rotations
c
      parameter(timedim=60)
      parameter(spacedim=1000)
      parameter(mindim=60)
      parameter(maxdim=1000)
      parameter(rotations=10)

c
c  Specify variables and arrays
c
c    Integer variables
c       m = timedim
c       n = spacedim
c       rot = rotations
c       q = smaller of timedim or spacedim (mindim)
c       i,j,l,pass = counters
c       
c    Arrays 
c
c       at = original data in format space-time
c       a = original data in format time-space
c       v = V array (**Note it's not VT) 
c       w = Sigma array: Singular values
c       rawPC = Unnormalized PCs in format time-time
c      
c       temp arrays used for data sorting


      integer m,n,q,i,j,l,pass
      integer rot
      real at(spacedim,timedim)
      real a(timedim,spacedim),v(spacedim,spacedim)
      real w(maxdim)
      real rawPC(timedim,timedim),sumPC
      real temp1,temp2(timedim),temp3(spacedim)


c Input and output files


      open (10,file='Atarray.bin',
     + form='unformatted', access='direct', 
     + recl=4*spacedim*timedim)

      open (11,file='Uarray_rot.bin',
     + form='unformatted', access='direct', 
     + recl=4*timedim*timedim)

      open (12,file='Varray_rot.bin',
     + form='unformatted', access='direct', 
     + recl=4*spacedim*spacedim)

      open (13, file='Sigmavalues.asci',
     + form='formatted')

      open (14,file='PCsraw_rot.bin',
     + form='unformatted', access='direct', 
     + recl=4*timedim*timedim)    



c Read in data (AT array)

      read(10,rec=1) ((at(j,i),j=1,spacedim),i=1,timedim)

      n=spacedim
      m=timedim


c Transpose data so it is in the form time-space

      do j=1,n
        do i=1,m
           a(i,j)=at(j,i)
        end do
      end do

c Call SVDcmp subroutine

      print*, 'To svdcmp'

      call svdcmp(a,m,n,m,n,w,v)
  
      print*, 'Exit svdcmp'


c  Sort the data from highest to lowest eigenvalue
c  Null space excluded

      q=mindim

      do pass=1,q
        do j=1,(q-pass)
          if(w(j).le.w(j+1)) then
            temp1=w(j)
            do i=1,m
              temp2(i)=a(i,j)
            end do
            do i=1,n
              temp3(i)=v(i,j)
            end do
            w(j)=w(j+1)
            do i=1,m
              a(i,j)=a(i,j+1)
            end do
            do i=1,n
              v(i,j)=v(i,j+1)
            end do
            w(j+1)=temp1
            do i=1,m
              a(i,j+1)=temp2(i)
            end do
            do i=1,n
              v(i,j+1)=temp3(i)
            end do
          end if
        end do
      end do   



c Call to do varimax rotation on eigenvectors

       print*, 'To varimax rotation'

       rot=int(rotations)

       call varimax(v,a,n,m,rot,rot)
     
       print*, 'Completed rotation on first',10 ,' eigenvectors'


c  Project original data onto EOFs to get unnormalized PCs
c  Need these to compute updated explained variance for REOFs

     
      do l=1,timedim
        do j=1,timedim
          sumPC=0
          do i=1,spacedim
            sumPC=v(i,l)*at(i,j)+sumPC
          end do
          rawPC(j,l)=sumPC
        end do
      end do      



c Output U and V arrays, Singular values, Raw PCs

      write(11,rec=1) ((a(i,j),i=1,timedim),j=1,timedim)
      write(14,rec=1) ((rawPC(j,l),j=1,timedim),l=1,timedim)
      write(12,rec=1) ((v(i,j),i=1,spacedim),j=1,spacedim)
      
      do i=1,q
        write(13,*) w(i)
      end do          

      print*, 'SVD Analysis with rotated eigenvectors complete'


      end


c  SVDcmp subroutine
c    Given a matrix (1:m,1:n) with physical dimensions mp by np,
c    this routine computes its singular value decomposition,
c    A = U W VT.  The matrix U replaces a on output.  The diagonal
c    matrix of singular values W is output as a vector w(1:n)
c    The matrix V (not the transpose VT) is the output as v(1:n,1:n) 
c
      SUBROUTINE svdcmp(a,m,n,mp,np,w,v)
      INTEGER m,mp,n,np,NMAX
      REAL a(mp,np),v(np,np),w(np)
      PARAMETER (NMAX=10000)
CU    USES pythag
      INTEGER i,its,j,jj,k,l,nm
      REAL anorm,c,f,g,h,s,scale,x,y,z,rv1(NMAX),pythag
      g=0.0
      scale=0.0
      anorm=0.0
      do 25 i=1,n
        l=i+1
        rv1(i)=scale*g
        g=0.0
        s=0.0
        scale=0.0
        if(i.le.m)then
          do 11 k=i,m
            scale=scale+abs(a(k,i))
11        continue
          if(scale.ne.0.0)then
            do 12 k=i,m
              a(k,i)=a(k,i)/scale
              s=s+a(k,i)*a(k,i)
12          continue
            f=a(i,i)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,i)=f-g
            do 15 j=l,n
              s=0.0
              do 13 k=i,m
                s=s+a(k,i)*a(k,j)
13            continue
              f=s/h
              do 14 k=i,m
                a(k,j)=a(k,j)+f*a(k,i)
14            continue
15          continue
            do 16 k=i,m
              a(k,i)=scale*a(k,i)
16          continue
          endif
        endif
        w(i)=scale *g
        g=0.0
        s=0.0
        scale=0.0
        if((i.le.m).and.(i.ne.n))then
          do 17 k=l,n
            scale=scale+abs(a(i,k))
17        continue
          if(scale.ne.0.0)then
            do 18 k=l,n
              a(i,k)=a(i,k)/scale
              s=s+a(i,k)*a(i,k)
18          continue
            f=a(i,l)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,l)=f-g
            do 19 k=l,n
              rv1(k)=a(i,k)/h
19          continue
            do 23 j=l,m
              s=0.0
              do 21 k=l,n
                s=s+a(j,k)*a(i,k)
21            continue
              do 22 k=l,n
                a(j,k)=a(j,k)+s*rv1(k)
22            continue
23          continue
            do 24 k=l,n
              a(i,k)=scale*a(i,k)
24          continue
          endif
        endif
        anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
25    continue
      do 32 i=n,1,-1
        if(i.lt.n)then
          if(g.ne.0.0)then
            do 26 j=l,n
              v(j,i)=(a(i,j)/a(i,l))/g
26          continue
            do 29 j=l,n
              s=0.0
              do 27 k=l,n
                s=s+a(i,k)*v(k,j)
27            continue
              do 28 k=l,n
                v(k,j)=v(k,j)+s*v(k,i)
28            continue
29          continue
          endif
          do 31 j=l,n
            v(i,j)=0.0
            v(j,i)=0.0
31        continue
        endif
        v(i,i)=1.0
        g=rv1(i)
        l=i
32    continue
      do 39 i=min(m,n),1,-1
        l=i+1
        g=w(i)
        do 33 j=l,n
          a(i,j)=0.0
33      continue
        if(g.ne.0.0)then
          g=1.0/g
          do 36 j=l,n
            s=0.0
            do 34 k=l,m
              s=s+a(k,i)*a(k,j)
34          continue
            f=(s/a(i,i))*g
            do 35 k=i,m
              a(k,j)=a(k,j)+f*a(k,i)
35          continue
36        continue
          do 37 j=i,m
            a(j,i)=a(j,i)*g
37        continue
        else
          do 38 j= i,m
            a(j,i)=0.0
38        continue
        endif
        a(i,i)=a(i,i)+1.0
39    continue
      do 49 k=n,1,-1
        do 48 its=1,30
          do 41 l=k,1,-1
            nm=l-1
            if((abs(rv1(l))+anorm).eq.anorm)  goto 2
            if((abs(w(nm))+anorm).eq.anorm)  goto 1
41        continue
1         c=0.0
          s=1.0
          do 43 i=l,k
            f=s*rv1(i)
            rv1(i)=c*rv1(i)
            if((abs(f)+anorm).eq.anorm) goto 2
            g=w(i)
            h=pythag(f,g)
            w(i)=h
            h=1.0/h
            c= (g*h)
            s=-(f*h)
            do 42 j=1,m
              y=a(j,nm)
              z=a(j,i)
              a(j,nm)=(y*c)+(z*s)
              a(j,i)=-(y*s)+(z*c)
42          continue
43        continue
2         z=w(k)
          if(l.eq.k)then
            if(z.lt.0.0)then
              w(k)=-z
              do 44 j=1,n
                v(j,k)=-v(j,k)
44            continue
            endif
            goto 3
          endif
          if(its.eq.30) pause 'no convergence in svdcmp'
          x=w(l)
          nm=k-1
          y=w(nm)
          g=rv1(nm)
          h=rv1(k)
          f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
          g=pythag(f,1.0)
          f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
          c=1.0
          s=1.0
          do 47 j=l,nm
            i=j+1
            g=rv1(i)
            y=w(i)
            h=s*g
            g=c*g
            z=pythag(f,h)
            rv1(j)=z
            c=f/z
            s=h/z
            f= (x*c)+(g*s)
            g=-(x*s)+(g*c)
            h=y*s
            y=y*c
            do 45 jj=1,n
              x=v(jj,j)
              z=v(jj,i)
              v(jj,j)= (x*c)+(z*s)
              v(jj,i)=-(x*s)+(z*c)
45          continue
            z=pythag(f,h)
            w(j)=z
            if(z.ne.0.0)then
              z=1.0/z
              c=f*z
              s=h*z
            endif
            f= (c*g)+(s*y)
            x=-(s*g)+(c*y)
            do 46 jj=1,m
              y=a(jj,j)
              z=a(jj,i)
              a(jj,j)= (y*c)+(z*s)
              a(jj,i)=-(y*s)+(z*c)
46          continue
47        continue
          rv1(l)=0.0
          rv1(k)=f
          w(k)=x
48      continue
3       continue
49    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.

C  Varimax subroutine
C
C  THIS SUBROUTINE ROTATES THE VECTORS X USING THE VARIMAX CRITERION. 
C  ON INPUT,  X CONTAINS THE N EIGENVECTORS, EACH OF LENGTH M, TO BE ROTATED. 
C  ON OUTPUT, X CONTAINS THE ROTATED SYSTEM.

      subroutine varimax(X,Y,M,K,N,NM)   

      INTEGER M,K,N,NM
      REAL X(M,N),Y(K,N)
     
      REAL H(M),U(M),V(M)
      REAL T(K)
      REAL A,B,C,D,VVOLD,VVNEW
      REAL PHI,COSPHI,SINPHI
      REAL SUMTERM,SUMTERM2,SUMTERM3
      INTEGER P,Q,L,MAXITER

C  Normalize components

      do P=1,M
        SUMTERM=0
        do Q=1,N
          SUMTERM=X(P,Q)**2+SUMTERM
        end do
        H(P)=sqrt(SUMTERM)
      end do

      do Q=1,N
        do P=1,M
          if(H(P).ne.0) then
            X(P,Q)=X(P,Q)/H(P)
          end if
        end do
      end do

C Varimax iteration for all pairs of vectors



      SUMTERM=0
      do P=1,M
        do Q=1,N
          SUMTERM=X(P,Q)**4+SUMTERM
        end do
      end do

      SUMTERM2=0
      SUMTERM3=0
      do Q=1,N
        do P=1,M
          SUMTERM2=X(P,Q)**2+SUMTERM2
        end do
        SUMTERM3=SUMTERM2**2+SUMTERM3
        SUMTERM2=0
      end do

      MAXITER=40
      VVOLD= M*SUMTERM-SUMTERM3

      do L=1,MAXITER

        do P=N,N-NM+2,-1
          do Q=P-1,N-NM+1,-1
            
             do I=1,M
               U(I)=X(I,P)**2-X(I,Q)**2
               V(I)=2*X(I,P)*X(I,Q)
             end do
           
             SUMTERM=0
             SUMTERM2=0
             do I=1,M
               A=U(I)+SUMTERM
               B=V(I)+SUMTERM2
             end do

             SUMTERM=0
             SUMTERM2=0
             do I=1,M
               SUMTERM=U(I)*U(I)-V(I)*V(I)+SUMTERM
               SUMTERM2=U(I)*V(I)+SUMTERM2
             end do

             C=SUMTERM-(A*A-B*B)/FLOAT(M)
             D=2*(SUMTERM2-A*B/FLOAT(M))

             PHI=.25*dble((atan2(D,C)))
             COSPHI=dble(cos(PHI))
             SINPHI=dble(sin(PHI))
             
             do I=1,M
               U(I)=X(I,P)*COSPHI+X(I,Q)*SINPHI
               X(I,Q)=-X(I,P)*SINPHI+X(I,Q)*COSPHI
               X(I,P)=U(I)
             end do
  
             do I=1,K
               T(I)=Y(I,P)*COSPHI+Y(I,Q)*SINPHI
               Y(I,Q)=-Y(I,P)*SINPHI+Y(I,Q)*COSPHI
               Y(I,P)=T(I)
             end do
          end do
        end do

        SUMTERM=0
        do P=1,M
          do Q=1,N
            SUMTERM=X(P,Q)**4+SUMTERM
          end do
        end do

        SUMTERM2=0
        SUMTERM3=0
        do Q=1,N
          do P=1,M
            SUMTERM2=X(P,Q)**2+SUMTERM2
          end do
          SUMTERM3=SUMTERM2**2+SUMTERM3
          SUMTERM2=0
        end do

        VVNEW= M*SUMTERM-SUMTERM3    


        if((VVNEW-VVOLD)/VVNEW.lt.(0.0001)) goto 50
          
        if(L.eq.MAXITER) then
          print*, 'VARIMAX rotation iteration did not coverge'
          goto 50
        endif

        VVOLD=VVNEW
      
      end do

50    do Q=1,N
       do P=1,M
         X(P,Q)=X(P,Q)*H(P)
       end do
      end do
    
      return
      END
C
C

      FUNCTION pythag(a,b)
      REAL a,b,pythag
      REAL absa,absb
      absa=abs(a)
      absb=abs(b)
      if(absa.gt.absb)then
        pythag=absa*sqrt(1.+(absb/absa)**2)
      else
        if(absb.eq.0.)then
          pythag=0.
        else
          pythag=absb*sqrt(1.+(absa/absb)**2)
        endif
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.