Subroutine MCA_SVD (A,B,C,Nx,Ny,Nt,Uk,Vk,Wk,Ak,Bk,Percent,Nmode,Miss) !Integer :: Nx, Ny, Nmode ! Nx must be less than or equal to Ny Real, Allocatable :: A(:,:), B(:,:) Real, Allocatable :: C(:,:), U(:,:), V(:,:), WW(:,:), W(:) Real, Allocatable :: Uk(:,:), Vk(:,:), Wk(:), Percent(:) !, WWk(:,:) Real, Allocatable :: Ak(:,:), Bk(:,:) Real :: Miss, tmp, tmp1, TotalVar Integer :: i, k if (.Not.Allocated(A)) Allocate (A(Nx,Nt)) if (.Not.Allocated(B)) Allocate (B(Ny,Nt)) if (.Not.Allocated(C)) Allocate (C(Nx,Ny)) if (.Not.Allocated(U)) Allocate (U(Nx,Ny)) if (.Not.Allocated(V)) Allocate (V(Ny,Ny)) if (.Not.Allocated(W)) Allocate (W(Ny)) if (.Not.Allocated(WW)) Allocate (WW(Ny,Ny)) if (.Not.Allocated(Uk)) Allocate (Uk(Nx,Nmode)) if (.Not.Allocated(Vk)) Allocate (Vk(Ny,Nmode)) if (.Not.Allocated(Wk)) Allocate (Wk(Nmode)) if (.Not.Allocated(Percent)) Allocate (Percent(Nmode)) if (.Not.Allocated(Ak)) Allocate (Ak(Nt,Nmode)) if (.Not.Allocated(Bk)) Allocate (Bk(Nt,Nmode)) !if (.Not.Allocated(WWk)) Allocate (WWk(Nmode,Nmode)) ! Input A(Nx,Nt), B(Ny,Nt), return to C(Nx,Ny) print*,"Calculating Covar matrx ... " Call Covar_Mtrx(A,B,C,Nt,Nx,Ny,Miss) print*,"Calculating Covar matrx ... Done" Call Write_C_Mtrx ("/work/projects/G09/SEACLID/MCA/C_matrx.nc") ! Input for SVD: C(Nx,Ny) where Nx<=Ny, return to U(Nx,Ny),V(Ny,Ny),W(Ny) print*,"Calculating SVD ..." Call svdcmp (C,U,V,W,Nx,Ny) print*,"Calculating SVD ... Done" Call Sort_svd (U,V,W,Nx,Ny) ! Extract leading Nmode: Do i=1,Nmode Uk(1:Nx,i)=U(1:Nx,i) Vk(1:Ny,i)=V(1:Ny,i) Wk(i) =W(i) Enddo TotalVar = Sum(W(1:Nx)) Do i=1,Nmode tmp=Sum(W(1:i)) Percent(i)=tmp/TotalVar Enddo !WWk=0 !Do i=1,Nmode; WWk(i,i)=Wk(i); Enddo ! Ak=AT.Uk Bk=BT.Vk Ak = MatMul(Transpose(A),Uk) Bk = MatMul(Transpose(B),Vk) Deallocate (C,U,V,W,WW) !print*,"Uk=",Uk !print*,"Vk=",Vk !print*,"Wk=",Wk !print*,"A=Uk*WWk*VTk=",MatMul(MatMul(Uk,WWk),Transpose(Vk)) End Subroutine MCA_SVD !---------------------------------------------------------------------------------------------- Subroutine Write_C_Mtrx (Fc) Character (Len=*) :: Fc Integer :: ix, iy Real :: x1d(Nx), y1d(Ny) Do ix=1,Nx; x1d(ix)=ix; Enddo Do iy=1,Ny; y1d(iy)=iy; Enddo Call wr_ncdf_2d (Fc,x1d,y1d,Nx,Ny,C,"cor","Covariance Matrix", "-","x","y", Miss) End Subroutine Write_C_Mtrx !--------------------------------------------------------------------------------- Subroutine Sort_svd (U,V,W,Nx,Ny) ! Sort svd values from hgihgest to lowest Integer :: q, i,j,n,m,k , Nx, Ny Real :: temp1,temp2(Nx),temp3(Ny) Real :: U(Nx,Ny), V(Ny,Ny), W(Ny) m=Nx; n=Ny q=min(m,n) do k=1,q do j=1,(q-k) if (w(j)<=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 u(i,j)=u(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 u(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 End Subroutine Sort_svd !--------------------------------------------------------------------------------- SUBROUTINE svdcmp (A,U,V,W,Nx,Ny) !(a,m,n,mp,np,w,v) INTEGER :: m,mp,n,np , Nx,Ny INTEGER :: i,its,j,jj,k,l,nm REAL :: anorm,c,f,g,h,s,xscale,x,y,z !, pythag !rv1(NMAX),pythag Real :: A(Nx,Ny), U(Nx,Ny), V(Ny,Ny), W(Ny), rv1(Ny) ! print*,"A-SVD=",A !M=Nt; N=Nx M=Nx; N=Ny g=0.0 xscale=0.0 anorm=0.0 do i=1,n !print*,"svdcmp-1 i=",i l=i+1 rv1(i)=xscale*g g=0.0 s=0.0 xscale=0.0 if(i.le.m)then do k=i,m xscale=xscale+abs(a(k,i)) enddo if(xscale.ne.0.0)then do k=i,m a(k,i)=a(k,i)/xscale s=s+a(k,i)*a(k,i) enddo f=a(i,i) g=-sign(sqrt(s),f) h=f*g-s a(i,i)=f-g do j=l,n s=0.0 do k=i,m s=s+a(k,i)*a(k,j) enddo f=s/h do k=i,m a(k,j)=a(k,j)+f*a(k,i) enddo enddo do k=i,m a(k,i)=xscale*a(k,i) enddo endif endif w(i)=xscale *g g=0.0 s=0.0 xscale=0.0 if((i.le.m).and.(i.ne.n))then do k=l,n xscale=xscale+abs(a(i,k)) enddo if(xscale.ne.0.0)then do k=l,n a(i,k)=a(i,k)/xscale s=s+a(i,k)*a(i,k) enddo f=a(i,l) g=-sign(sqrt(s),f) h=f*g-s a(i,l)=f-g do k=l,n rv1(k)=a(i,k)/h enddo do j=l,m s=0.0 do k=l,n s=s+a(j,k)*a(i,k) enddo do k=l,n a(j,k)=a(j,k)+s*rv1(k) enddo enddo do k=l,n a(i,k)=xscale*a(i,k) enddo endif endif anorm=max(anorm,(abs(w(i))+abs(rv1(i)))) enddo do i=n,1,-1 !print*,"svdcmp-2 i=",i if(i.lt.n)then if(g.ne.0.0)then do j=l,n v(j,i)=(a(i,j)/a(i,l))/g enddo do j=l,n s=0.0 do k=l,n s=s+a(i,k)*v(k,j) enddo do k=l,n v(k,j)=v(k,j)+s*v(k,i) enddo enddo endif do j=l,n v(i,j)=0.0 v(j,i)=0.0 enddo endif v(i,i)=1.0 g=rv1(i) l=i enddo do i=min(m,n),1,-1 !print*,"svdcmp-3 i=",i l=i+1 g=w(i) do j=l,n a(i,j)=0.0 enddo if(g.ne.0.0)then g=1.0/g do j=l,n s=0.0 do k=l,m s=s+a(k,i)*a(k,j) enddo f=(s/a(i,i))*g do k=i,m a(k,j)=a(k,j)+f*a(k,i) enddo enddo do j=i,m a(j,i)=a(j,i)*g enddo else do j= i,m a(j,i)=0.0 enddo endif a(i,i)=a(i,i)+1.0 enddo do k=n,1,-1 !print*,"svdcmp-4 i=",i do its=1,30 do 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 enddo 1 c=0.0 s=1.0 do 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 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) enddo enddo 2 z=w(k) if(l.eq.k)then if(z.lt.0.0)then w(k)=-z do j=1,n v(j,k)=-v(j,k) enddo endif goto 3 endif if(its.eq.50) print*, "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 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 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) enddo 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 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) enddo enddo rv1(l)=0.0 rv1(k)=f w(k)=x enddo 3 continue enddo U(:,:)=A(:,:) !print*,"U-SVD=",U END SUBROUTINE svdcmp !--------------------------------------------------------------------------------- Real 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 END FUNCTION pythag !--------------------------------------------------------------------------------- SUBROUTINE Covar_Mtrx(X,Y,C,Nt,Nx,Ny,Miss) Integer :: Nt, Nx, Ny, N1, L, J, K, IO, I Real :: X(Nx,Nt), Y(Ny,Nt), C(Nx,Ny), Xp(Nx,Nt), Yp(Ny,Nt) Real :: XBAR(Nx), Stdx(Nx), YBar(Ny), Stdy(Ny) Real :: Sx, Xtb, Sy, Ytb, Miss !--------------------------------------------------- ! Anomaly: Do j=1, Nx Sx=0; Xtb=0; N1=0 Do i=1, Nt if (X(j,i)/=Miss) then n1=n1+1; Xtb=Xtb+X(j,i); Sx=Sx+X(j,i)*X(j,i) Endif Enddo if (n1>nt/2) then XBAR(j)=Xtb/real(n1); STDx(j)=Sqrt(Sx/real(n1)-XBAR(j)*XBAR(j)) endif Do i=1, Nt if (X(j,i)==Miss) X(j,i)=XBAR(j) ! Fill-up Miss by mean values Xp(j,i) = (X(j,i)-XBAR(j)) !if (standard==1) Y(i,j) = Y(i,j)/STD(j) Enddo Enddo ! Do j=1, Ny Sx=0; Xtb=0; N1=0 Do i=1, Nt if (Y(j,i)/=Miss) then n1=n1+1; Xtb=Xtb+Y(j,i); Sx=Sx+Y(j,i)*Y(j,i) Endif Enddo if (n1>nt/2) then Ybar(j)=Xtb/real(n1); STDy(j)=Sqrt(Sx/real(n1)-Ybar(j)*Ybar(j)) endif Do i=1, Nt if (Y(j,i)==Miss) Y(j,i)=Ybar(j) ! Fill-up Miss by mean values Yp(j,i) = (Y(j,i)-Ybar(j)) !if (standard==1) Y(i,j) = Y(i,j)/STD(j) Enddo Enddo C = MatMul(Xp,Transpose(Yp))/Real(Nt) End SUBROUTINE Covar_Mtrx !---------------------------------------------------------------------------------