Real, Parameter :: Miss=-99.

      Real, Allocatable :: Xbar(:), Std(:) ! D(:), V(:)
      Real, Allocatable :: X(:,:), Y(:,:), YT(:,:), Xorg(:,:)
      ! SVD:
      Real, Allocatable :: A(:,:), U(:,:), V(:,:), WW(:,:), W(:), RV1(:), UU(:,:) 
      Real, Allocatable :: Uk(:,:), Vk(:,:), Wk(:), Percent(:), WWk(:,:)

      Real, Allocatable :: Xnew(:), Fnew(:)
      ! Jacobi:
      Real, Allocatable :: RR(:,:)
      
      Real, Allocatable :: Fx(:,:), Ax(:,:)   ! X = Fx*A'x
      Real, Allocatable :: EOF(:,:), PCs(:,:)
      !For testing:
      Real, Allocatable :: FF(:,:), AA(:,:)
      Real, Allocatable :: X1(:,:), EVtor2(:,:)
      Real, Allocatable :: YdX(:,:)

      Real, Allocatable :: EVtor(:,:), Lamda(:,:)
      Real, Allocatable :: EValues(:)

      Integer :: Ny, Nx
      Integer :: Mcount, Ncount, switch, Standard, Nmode, Nthre
      Real :: Threshold
      Character (Len=250) :: Fname
      Real :: Factor

      Namelist /Model_Info/ Fname, switch, standard, PC_Stard, Nmode 

      Integer :: i,j,k

      Open (1,file="Namelist_SVD_New.txt",Status="old")
      Read (1, Model_Info)

      Nx=10; Ny=15
      Factor = Sqrt(Real(Nx))

      Allocate (X(Nx,Ny), Y(Nx,Ny), YT(Ny,Nx))   ! Origin data, Centered or standardized data and its transpose
      Allocate (Xorg(Nx,Ny), Xnew(Ny),Fnew(Nmode))

      ! SVD: A(Nt,Nx) = U(Nt,Nx).W(Nx).VT(Nx,Nx) 
      Allocate (A(Nx,Ny), U(Nx,Ny), V(Ny,Ny), WW(Ny,Ny), W(Ny), RV1(Max(Nx,Ny)))
      Allocate (Uk(Nx,Nmode), Vk(Ny,Nmode), Wk(Nmode), Percent(Nmode), WWk(Nmode,Nmode))
      Allocate (FF(Nx,Nx), AA(Ny,Ny),UU(Nx,Nx))

      Open (1,file=fname)
      read (1,*)
      Do i=1, Nx
         read (1,*) (X(i,j),j=1,Ny)
      Enddo

      if (switch==1) then
         Mcount=Nx; Ncount=Ny
      else if (switch==0) then
         Mcount=Ny; Ncount=Nx
      endif
      !For testing:
      !Allocate (FF(Mcount,Mcount), AA(Mcount,Mcount))
      Allocate (X1(Nx,Ny), EVtor2(Mcount,Mcount))
      Allocate (YdX(Nx,Ny))

      Allocate (Xbar(Mcount), Std(Mcount))                      ! TB, DLC
      Allocate (EValues(Mcount))                                ! Tri rieng va Vecto rieng (1D)
      Allocate (Fx(Ncount,Mcount), Ax(Mcount,Mcount))           ! EOF va PCA
      Allocate (EVtor(Mcount,Mcount))                ! Working arr va matran Vecto rieng
      Allocate (Lamda(Mcount,Mcount))
      ! Jacobi:
      Allocate (RR(Mcount,Mcount))

      ! SVD:

     ! Test with SVD using centered matrix:
      Call Covar(Nx,Ny)
      A(:,:)=Y(:,:)/(Nx)

      ! Input for SVD: A(Nt,Nx) where Nt<Nx
      print*,"A-Call=",A
      Call SVD (A,Nx,Ny,Uk,Vk,Wk,Percent,Nmode)

         print*,"Uk=",Uk
         print*,"Vk=",Vk
         print*,"Wk=",Wk


      !Call svdcmpi (A,U,Nt,Nx)
      !Call Sort_svd (U,V,W,Nt,Nx)

      WW=0; !Do i=1,Nx;W(i)=W(i)**2; Enddo
      Do i=1,Nmode; WWk(i,i)=Wk(i); Enddo

      A=MatMul(MatMul(Uk,WWk),Transpose(Vk))
print*,"A=Uk*WWk*VTk=",A
      Y=MatMul(Uk,WWk)
      FF=MatMul(Y,Transpose(Y))
      UU=MatMul(Uk,Transpose(Uk))
      AA=MatMul(Vk,Transpose(Vk))
    
      print*,Wk 
      print*,"Sum of Cova=", Sum(Wk) 
      Call Print_test_SVD

      !Call Covar(Nt,Nx)
      !Call Eigen
      !Call Eigsrt
      K = Nmode



      Contains
      !====================================================
      Subroutine Print_test_SVD
        print*,"Xorg"
        Do i=1, Nx
           Do j=1, Ny; X(i,j)=X(i,j)-Xbar(j); Enddo
           print "(20F10.2)", (X(i,j),j=1, Ny)
        Enddo
        print*,"X-estimated"
        Do i=1, Nx
           print "(20F10.2)", (A(i,j)*(Nx),j=1, Ny)
        Enddo

        print*,"U"
        Do i=1, Nx !Mcount
           print "(20F10.2)", (U(i,j),j=1,Ny) ! Mcount)
        Enddo
        print*,"W"
        Do i=1, Ny !Mcount
           print "(20F10.2)", (WW(i,j),j=1,Ny) ! Mcount)
        Enddo

        print*,"FF"
        Do i=1, Nx !Mcount
           print "(20F10.2)", (FF(i,j),j=1,Nx) ! Mcount)
        Enddo
        print*,"AA"
        Do i=1, Ny
           print "(20F10.2)", (AA(i,j),j=1,Ny)! Mcount)
        Enddo
        print*,"UU"
        Do i=1, Nx
           print "(20F10.2)", (UU(i,j),j=1,Nx)! Mcount)
        Enddo

      End Subroutine Print_test_SVD
      !---------------------------------------------------------------------------------
      Subroutine Write_test
      print*,"PCA"
      Do i=1,Nx
         print "(20F10.2)", (PCs(i,j),j=1,K) ! Mcount)
      Enddo
      print*,"EOF"
      Do i=1, Nt
         print "(20F10.2)", (EOF(i,j),j=1,K) ! Mcount)
      Enddo
      print*,"AA"
      Do i=1, K !Mcount
         print "(20F10.2)", (AA(i,j),j=1,K) ! Mcount)
      Enddo
      print*,"FF"
      Do i=1, K
         print "(20F10.2)", (FF(i,j),j=1,K)! Mcount)
      Enddo
      print*,"Y"
      Do i=1, Ncount
         print "(20F10.2)", (Y(i,j),j=1, Mcount)
      Enddo
      print*,"X1"
      Do i=1, Ncount
         print "(20F10.2)", (X1(i,j),j=1, Mcount)
      Enddo
      print*,"Xorg"
      Do i=1, Ncount
         print "(20F10.2)", (X(i,j),j=1, Mcount)
      Enddo
      print*,"X-estimated"
      Do i=1, Ncount
         print "(20F10.2)", (Xorg(i,j),j=1, Mcount)
      Enddo

      print*,"EVtor2"
      Do i=1, Mcount
         print "(20F10.2)", (EVtor2(i,j),j=1, Mcount)
      Enddo
      End Subroutine Write_test
      !---------------------------------------------------

      SUBROUTINE Covar(Nx,Ny)
      Integer :: Ny, Nx, N1, L, J, K, IO, I
      Real :: Sx, Xtb
      !---------------------------------------------------
      Do j=1, Ny
         Sx=0; Xtb=0; N1=0
         Do i=1, Nx
            if (X(i,j)/=Miss) then
               n1=n1+1; Xtb=Xtb+X(i,j); Sx=Sx+X(i,j)*X(i,j)
            Endif
         Enddo
         if (n1>nx/2) then
            XBAR(j)=Xtb/real(n1); STD(j)=Sqrt(Sx/real(n1)-XBAR(j)*XBAR(j))
         endif
         Do i=1, Nx
            if (X(i,j)==Miss) X(i,j)=XBAR(j) ! Fill-up Miss by mean values
            Y(i,j) = (X(i,j)-XBAR(j)) 
            if (standard==1) Y(i,j) = Y(i,j)/STD(j)
         Enddo
      Enddo
      if (switch==1) then
         RR = MatMul(Y,Transpose(Y))
      else if (switch==0) then
         RR = MatMul(Transpose(Y),Y)
      endif
      RR(:,:)=RR(:,:)/real(Nx)
      End SUBROUTINE Covar

     !-----------------------------------------------------------------------
     Include "../SST/Sub_SVD.f90"

      END