Real, Parameter :: Miss=-99.

      Real, Allocatable :: Xbar(:), Std(:) ! D(:), V(:)
      Real, Allocatable :: X(:,:), Y(:,:), YT(:,:), Xorg(:,:)
      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(:)
      Real*8, Allocatable :: E_lapack(:),V_lapack(:,:)
      ! New:
      double precision, Allocatable :: Cor(:,:), Evect(:,:)
      double precision :: abserr

      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_PCA_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))

      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))

      Call Covar(Nt,Nx)
!      Call Eigen

!          EValues(:)=D(:)
!          EVtor(:,:)=V(:,:)

      Allocate (E_lapack(Mcount),V_lapack(Mcount,Mcount))
      Call Eigen_Lapack(RR, Mcount, E_lapack,V_lapack)
      EValues(:)=E_lapack(:); EVtor(:,:)=V_lapack(:,:)

      ! New:
      !Allocate (Cor(Mcount,Mcount), Evect(Mcount,Mcount))
      !Cor(:,:)=RR(:,:)
      !abserr=1E-4
      !Call Eigen_New (Cor,Evect,abserr,Mcount)
      !Do i=1,Mcount
      !   EValues(i)=Cor(i,i)
      !Enddo
      !EVtor(:,:)=Evect(:,:)
      !-----------------------

      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
!      Call Calc_EOF_PCA
      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)

      print*,EValues
      print*,"Sum of Cova=",Sum(EValues)
      Call Write_test


      Contains
      !====================================================
      Subroutine Eigen_New (a,x,abserr,n)
      !===========================================================
      ! Evaluate eigenvalues and eigenvectors
      ! of a real symmetric matrix a(n,n): a*x = lambda*x 
      ! method: Jacoby method for symmetric matrices 
      ! Alex G. (December 2009)
      !-----------------------------------------------------------
      ! input ...
      ! a(n,n) - array of coefficients for matrix A
      ! n      - number of equations
      ! abserr - abs tolerance [sum of (off-diagonal elements)^2]
      ! output ...
      ! a(i,i) - eigenvalues
      ! x(i,j) - eigenvectors
      ! comments ...
      !===========================================================
         implicit none
         integer i, j, k, n
         double precision a(n,n),x(n,n)
         double precision abserr, b2, bar
         double precision beta, coeff, c, s, cs, sc

         ! initialize x(i,j)=0, x(i,i)=1
         ! *** the array operation x=0.0 is specific for Fortran 90/95
         x = 0.0
         do i=1,n
            x(i,i) = 1.0
         end do

         ! find the sum of all off-diagonal elements (squared)
         b2 = 0.0
         do i=1,n
            do j=1,n
               if (i.ne.j) b2 = b2 + a(i,j)*a(i,j) !**2
            end do
         end do

         if (b2 <= abserr) return
            ! average for off-diagonal elements /2
         bar = 0.5*b2/float(n*n)

         do while (b2.gt.abserr)
            print*,b2,abserr
            do i=1,n-1
               do j=i+1,n
                  if (a(j,i)**2 <= bar) cycle  ! do not touch small elements
                  b2 = b2 - 2.0*a(j,i)**2
                  bar = 0.5*b2/float(n*n)
                  ! calculate coefficient c and s for Givens matrix
                  beta = (a(j,j)-a(i,i))/(2.0*a(j,i))
                  coeff = 0.5*beta/sqrt(1.0+beta**2)
                  s = sqrt(max(0.5+coeff,0.0))
                  c = sqrt(max(0.5-coeff,0.0))
                  ! recalculate rows i and j
                  do k=1,n
                     cs =  c*a(i,k)+s*a(j,k)
                     sc = -s*a(i,k)+c*a(j,k)
                     a(i,k) = cs
                     a(j,k) = sc
                  end do
                  ! new matrix a_{k+1} from a_{k}, and eigenvectors 
                  do k=1,n
                     cs =  c*a(k,i)+s*a(k,j)
                     sc = -s*a(k,i)+c*a(k,j)
                     a(k,i) = cs
                     a(k,j) = sc
                     cs =  c*x(k,i)+s*x(k,j)
                     sc = -s*x(k,i)+c*x(k,j)
                     x(k,i) = cs
                     x(k,j) = sc
                  end do
               end do
            end do
         end do
      End Subroutine Eigen_New
      !----------------------------------------------------
      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

     !-----------------------------------------------------------------------
subroutine Eigen_Lapack(A, M, E, V)
  implicit none
  integer, intent(in) :: M
  real,  intent(in)  :: A(M,M)
  real*8,  intent(out) :: E(M)
  real*8,  intent(out) :: V(M,M)

  integer :: info, lwork, liwork, istat
  real*8  :: wkopt
  integer :: iwkopt
  real*8,  allocatable :: work(:)
  integer, allocatable :: iwork(:)
  integer :: i, j

  interface
     subroutine dsyevd(jobz,uplo,n,a,lda,w,work,lwork,iwork,liwork,info)
       character(len=1), intent(in) :: jobz, uplo
       integer, intent(in) :: n, lda, lwork, liwork
       double precision, intent(inout) :: a(lda,*)
       double precision, intent(out) :: w(*)
       double precision, intent(inout) :: work(*)
       integer, intent(inout) :: iwork(*)
       integer, intent(out) :: info
     end subroutine dsyevd
  end interface

  ! Copy input into V, enforce symmetry for safety
  V(:,:) = A(:,:)
  do j = 1, M
     do i = j+1, M
        V(i,j) = 0.5d0*(V(i,j) + V(j,i))
        V(j,i) = V(i,j)
     end do
  end do

  ! Workspace query
  lwork  = -1
  liwork = -1

  allocate(work(1),  stat=istat); if (istat/=0) stop "Eigen: alloc work(1) failed"
  allocate(iwork(1), stat=istat); if (istat/=0) stop "Eigen: alloc iwork(1) failed"

  call dsyevd('V','U', M, V, M, E, work, lwork, iwork, liwork, info)

  wkopt  = work(1)
  iwkopt = iwork(1)

  deallocate(work,  stat=istat)
  deallocate(iwork, stat=istat)

  if (info /= 0) then
     write(*,*) "Eigen: DSYEVD workspace query failed, info=", info
     stop
  end if

  ! Allocate work arrays (ceiling + small buffer)
  lwork  = max(1, int(wkopt + 0.999999d0) + 64)
  liwork = max(1, iwkopt + 64)

  allocate(work(lwork),  stat=istat); if (istat/=0) stop "Eigen: alloc work failed"
  allocate(iwork(liwork),stat=istat); if (istat/=0) stop "Eigen: alloc iwork failed"

  call dsyevd('V','U', M, V, M, E, work, lwork, iwork, liwork, info)

  deallocate(work,  stat=istat)
  deallocate(iwork, stat=istat)

  if (info /= 0) then
     write(*,*) "Eigen: DSYEVD failed, info=", info
     stop
  end if

end subroutine Eigen_Lapack

      END