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 :: 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 :: Nt, 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)

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

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

      ! SVD: A(Nt,Nx) = U(Nt,Nx).W(Nx).VT(Nx,Nx) 
      Allocate (A(Nt,Nx), U(Nt,Nx), V(Nx,Nx), WW(Nx,Nx), W(Nx), RV1(Max(Nt,Nx)))
      Allocate (FF(Nt,Nt), AA(Nx,Nx),UU(Nt,Nt))

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

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

      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(Nt,Nx)
      A(:,:)=Y(:,:)/(Nt)

      ! Input for SVD: A(Nt,Nx) where Nt<Nx
      Call svdcmp
      Call Sort_svd
      WW=0; !Do i=1,Nx;W(i)=W(i)**2; Enddo
      Do i=1,Nx; WW(i,i)=W(i); Enddo
      A=MatMul(MatMul(U,WW),Transpose(V))
      Y=MatMul(U,WW)
      FF=MatMul(Y,Transpose(Y))
      UU=MatMul(U,Transpose(U))
      AA=MatMul(V,Transpose(V))
    
      print*,W 
      print*,"Sum of Cova=", Sum(W) 
      Call Print_test_SVD

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

      !Allocate (EOF(Nt,K), PCs(Nx,K))

      !Call Calc_AX

      !if (PC_Stard==1) then
      !   Call Calc_EOF_PCA_Standard
      !else
      !   Call Calc_EOF_PCA
      !endif
      !FF=MatMul(Transpose(EOF),EOF)
      !AA=MatMul(Transpose(PCs),PCs)
      !if (switch==0) then
      !   X1=MatMul(EOF,Transpose(PCs))*Factor
      !else if (switch==1) then
      !   X1=MatMul(EOF,Transpose(PCs))*Factor
         !Do j=1,Nx
         !   Do i=1,Nt
         !      YdX(i,j)=Y(i,j)/X1(i,j)
         !   Enddo
         !Enddo
      !endif
      ! Restore the origin values of X:
      
         !Do j=1,Nx
         !   Do i=1,Nt
         !      if (standard==1) then
         !         Xorg(i,j)=X1(i,j)*Std(j)+Xbar(j)
         !      else
         !         Xorg(i,j)=X1(i,j)+Xbar(j)
         !      endif 
         !   Enddo
         !Enddo
      
      !if (standard==1) then
      !   Do j=1, Mcount
         !   X1(:,j)=X1(:,j)*Factor
      !   Enddo
      !endif
      !EVtor2=MatMul(Transpose(EVtor),EVtor)

      !Call Write_test


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

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

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

      End Subroutine Print_test_SVD
      !---------------------------------------------------------------------------------
      Subroutine Sort_svd
         ! Sort svd values from hgihgest to lowest
         Integer :: q, i,j,n,m,pass
         Real :: temp1,temp2(Nt),temp3(Nx)

         m=Nt; n=Nx
         q=min(m,n)

         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
                     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,m,n,mp,np,w,v)
      INTEGER :: m,mp,n,np
      !REAL a(mp,np),v(np,np),w(np)
      !PARAMETER (NMAX=500)
!CU    USES pythag
      INTEGER :: i,its,j,jj,k,l,nm
      REAL :: anorm,c,f,g,h,s,scale,x,y,z !, pythag !rv1(NMAX),pythag

      M=Nt; N=Nx

      g=0.0
      scale=0.0
      anorm=0.0
      do 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 k=i,m
            scale=scale+abs(a(k,i))
          enddo
          if(scale.ne.0.0)then
            do k=i,m
              a(k,i)=a(k,i)/scale
              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)=scale*a(k,i)
            enddo
          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 k=l,n
            scale=scale+abs(a(i,k))
          enddo
          if(scale.ne.0.0)then
            do k=l,n
              a(i,k)=a(i,k)/scale
              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)=scale*a(i,k)
            enddo
          endif
        endif
        anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
      enddo
      do i=n,1,-1
        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
        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
        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.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 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(:,:)
      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 Eigen 
      INTEGER :: n,np,nrot
      INTEGER :: i,ip,iq,j
      REAL :: c,g,h,s,sm,t,tau,theta,tresh 
      Real, Allocatable :: A(:,:), D(:), V(:,:), B(:), Z(:)

      Allocate (A(Mcount,Mcount),D(Mcount),V(Mcount,Mcount),B(Mcount),Z(Mcount))
      A(:,:) = RR(:,:); N=Mcount
      do  ip=1,n
        do  iq=1,n
          v(ip,iq)=0.
        enddo
        v(ip,ip)=1.
      enddo
      do  ip=1,n
        b(ip)=a(ip,ip)
        d(ip)=b(ip)
        z(ip)=0.
      enddo
      nrot=0
      do  i=1,50
        sm=0.
        do  ip=1,n-1
          do iq=ip+1,n
            sm=sm+abs(a(ip,iq))
          enddo
        enddo
        if(sm.eq.0.) then  ! Normal terminal 
          EValues(:)=D(:)
          EVtor(:,:)=V(:,:)
          Deallocate (A, D, V, B, Z)
          return
        endif
        if(i.lt.4)then
          tresh=0.2*sm/n**2
        else
          tresh=0.
        endif
        do  ip=1,n-1
          do  iq=ip+1,n
            g=100.*abs(a(ip,iq))
            if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))).and. &
               (abs(d(iq))+g.eq.abs(d(iq))))then
              a(ip,iq)=0.
            else if(abs(a(ip,iq)).gt.tresh)then
              h=d(iq)-d(ip)
              if(abs(h)+g.eq.abs(h))then
                t=a(ip,iq)/h
              else
                theta=0.5*h/a(ip,iq)
                t=1./(abs(theta)+sqrt(1.+theta**2))
                if(theta.lt.0.)t=-t
              endif
              c=1./sqrt(1+t**2)
              s=t*c
              tau=s/(1.+c)
              h=t*a(ip,iq)
              z(ip)=z(ip)-h
              z(iq)=z(iq)+h
              d(ip)=d(ip)-h
              d(iq)=d(iq)+h
              a(ip,iq)=0.
              do  j=1,ip-1
                g=a(j,ip)
                h=a(j,iq)
                a(j,ip)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
              enddo
              do  j=ip+1,iq-1
                g=a(ip,j)
                h=a(j,iq)
                a(ip,j)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
              enddo
              do  j=iq+1,n
                g=a(ip,j)
                h=a(iq,j)
                a(ip,j)=g-s*(h+g*tau)
                a(iq,j)=h+s*(g-h*tau)
              enddo
              do  j=1,n
                g=v(j,ip)
                h=v(j,iq)
                v(j,ip)=g-s*(h+g*tau)
                v(j,iq)=h+s*(g-h*tau)
              enddo
              nrot=nrot+1
            endif
          enddo
        enddo
        do  ip=1,n
          b(ip)=b(ip)+z(ip)
          d(ip)=b(ip)
          z(ip)=0.
        enddo
      enddo
      Deallocate (A, D, V, B, Z)
      ! Not successful
      return
      END SUBROUTINE Eigen
      !---------------------------------------------------
      SUBROUTINE Eigsrt 
      INTEGER :: n,np
      Real, Allocatable :: D(:), V(:,:)
      INTEGER :: i,j,k
      REAL :: p

      Allocate (D(Mcount),V(Mcount,Mcount))
      N=Mcount
      D(:)=EValues(:); V(:,:)=EVtor(:,:)

      do  i=1,n-1
        k=i
        p=d(i)
        do  j=i+1,n
          if(d(j).ge.p)then
            k=j
            p=d(j)
          endif
        enddo
        if(k.ne.i)then
          d(k)=d(i)
          d(i)=p
          do  j=1,n
            p=v(j,i)
            v(j,i)=v(j,k)
            v(j,k)=p
          enddo
        endif
      enddo
      EValues(:)=D(:)
      EVtor(:,:)=V(:,:)
      Deallocate (D, V)
      return
      END SUBROUTINE Eigsrt
      !---------------------------------------------------
      Subroutine Calc_AX
      Integer :: i,j
      Real :: SQ
      DO j=1,K 
         SQ=Sqrt(EValues(j)) 
         DO i=1,Mcount 
            Ax(I,J)=EVtor(i,j)*SQ
         ENDDO
      ENDDO
      End Subroutine Calc_AX
      !---------------------------------------------------------------------------------
      Subroutine Calc_EOF_PCA_Standard
      Integer :: i,j
      Real :: Modul_Fx
        if (switch==0) then
           Ax(1:Nx,1:K)=EVtor(1:Nx,1:K)
           Fx(1:Nt,1:K) = MatMul(Y(1:Nt,1:Nx),Ax(1:Nx,1:K))
           Do j=1,K
              !Fx(:,j)=Fx(:,j)/(EValues(j))/Factor ! Fx Chuan hoa
              Fx(:,j)=Fx(:,j)/Factor               ! Voi PCs chuan hoa (Ax=EVtor)
           Enddo
           EOF(1:Nt,1:K) = Fx(1:Nt,1:K)
           PCs(1:Nx,1:K) = Ax(1:Nx,1:K)
        else if (switch==1) then
           Do j=1, K
              !Ax(:,j)=Ax(:,j)/sqrt(EValues(j))
              Lamda(j,j)=Sqrt(EValues(j))
           Enddo
           YT = Transpose(Y)
           !Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),Ax(1:Nt,1:K))
           Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),EVtor(1:Nt,1:K))
           Do j=1,K
              Modul_Fx=0
              Do i=1, Ncount
                 Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j)
              Enddo
              Modul_Fx = Sqrt(Modul_Fx)
              Fx(:,j)=Fx(:,j)/Modul_Fx
           Enddo
           !Fx=MatMul(Fx(1:Ncount,1:K),Lamda(1:K,1:K))
           EOF(1:Nt,1:K) = Ax(1:Nt,1:K)
           PCs(1:Nx,1:K) = Fx(1:Ncount,1:K)
        endif
      End Subroutine Calc_EOF_PCA_Standard
      !---------------------------------------------------
      Subroutine Calc_EOF_PCA
      Integer :: i,j
      Real :: Modul_Fx
        if (switch==0) then
           Fx(1:Nt,1:K) = MatMul(Y(1:Nt,1:Nx),Ax(1:Nx,1:K))
           Do j=1,K
              Fx(:,j)=Fx(:,j)/(EValues(j)) /Factor
           Enddo
           EOF(1:Nt,1:K) = Fx(1:Nt,1:K)
           PCs(1:Nx,1:K) = Ax(1:Nx,1:K)
        else if (switch==1) then
           Do j=1, K
              Ax(:,j)=Ax(:,j)/sqrt(EValues(j))
              Lamda(j,j)=Sqrt(EValues(j))
           Enddo
           YT = Transpose(Y)
           Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),Ax(1:Nt,1:K))
           Do j=1,K
              Modul_Fx=0
              Do i=1, Ncount
                 Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j)
              Enddo
              Modul_Fx = Sqrt(Modul_Fx)
              Fx(:,j)=Fx(:,j)/Modul_Fx
           Enddo
           Fx=MatMul(Fx(1:Ncount,1:K),Lamda(1:K,1:K))
           EOF(1:Nt,1:K) = Ax(1:Nt,1:K)
           PCs(1:Nx,1:K) = Fx(1:Ncount,1:K)
        endif
      End Subroutine Calc_EOF_PCA

      !---------------------------------------------------
      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(Nt,Nx)
      Integer :: Nt, Nx, N1, L, J, K, IO, I
      Real :: Sx, Xtb
      !---------------------------------------------------
      Do j=1, Nx
         Sx=0; Xtb=0; N1=0
         Do i=1, Nt
            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>nt/2) then
            XBAR(j)=Xtb/real(n1); STD(j)=Sqrt(Sx/real(n1)-XBAR(j)*XBAR(j))
         endif
         Do i=1, Nt
            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(Nt)
      End SUBROUTINE Covar

     !-----------------------------------------------------------------------

      END