program A_matrix_eigen

!**********************************************************************
! *              .      .    .
! * PROGRAM:    A_matrix_eigen
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-08-28
! *
! * ABSTRACT:   Eigenvalue decomposition (EVD) of the A matrix.
! * Perform an EVD of matrix A (NENS x NENS).
! * Save the matrix sqrt(I+A)**(-1/2) in full storage mode (simpler for subsequent use).
! * EVD performed using LAPACK and BLAS math libraries.
! * Employs the LAPACK program DSPEV  (double precision).
! *
! * PROGRAM LOG:
! *
! * 08/28/2003    ..... M. ZUPANSKI:
! * 10/08/2004    ..... M. ZUPANSKI: reunite IBMSP and NON_IBMSP pregrams
! *
! **********************************************************************
      real*8,parameter :: W_bound=1.d-7     !! Eigenvalue threshold - [IF W<W_bound THEN W=0.0]
      integer,parameter :: matrixa=21       !! Matrix A input file
      integer,parameter :: ia_eigen=51      !! Eigen (I+A)**(-1/2) file (for post-processing)
      integer,parameter :: ipa_out=52       !! (I+A)**(-1/2) output file
      integer,parameter :: idim=53
      integer,parameter :: i_evdA=744       !! Output eigenvalues 
      integer,parameter :: i_evdIPA=745     !! Output (I+A)**(-1/2) eigenvalues
!------------------------------------------------------------------------
      integer :: NENS,NENS_START,NALL
      integer :: NP,LDZ,mtot
!------------------------------
      real*8,dimension(:),allocatable    :: W_reverse
      real*8,dimension(:,:),allocatable  :: Z_reverse
      real*8,dimension(:),allocatable    :: W_scaled
      real,dimension(:,:),allocatable    :: Z_single
      real,dimension(:),allocatable      :: W_single
      real,dimension(:,:),allocatable    :: d_out
      real,dimension(:),allocatable      :: AP_single

      real :: trace,trace_inv,obsinf_cont,entropy

!
!      DECLARE NAMELIST
!
      NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS

!==============start calculation===================

!-----Read number of observations from res_cntrl--------------------
      CLOSE(idim)
      OPEN(idim,file='res_cntrl',status='unknown',form='unformatted')
      read(idim) mtot
      CLOSE(idim,status='KEEP')

!-- read ensemble size
       REWIND 15
       READ(15,ENSEMBLE_SIZE)

       write(*,*) "A_matrix_eigen: NENS_START, NENS=",NENS_START,NENS

!- Total number of degrees of freedom
       NALL=NENS-NENS_START+1

!====================================================================
!-- define dimensions -----------
!
       NP=(NALL+1)*NALL/2
       LDZ=NALL

!----------------------------------------------
     allocate(W_reverse(1:NALL))
     allocate(Z_reverse(1:LDZ,1:NALL))
     allocate(W_scaled(1:NALL))
!====================================================================
!
!   READ THE A-MATRIX ELEMENTS  (single precision) 
!
      allocate(AP_single(1:NP))
      CLOSE(matrixa)
      OPEN(matrixa,file='matrix_A',status='unknown',form='unformatted')
      read(matrixa) (AP_single(k),k=1,NP)
      CLOSE(matrixa,status='KEEP')

      CALL CALCULATE_EVD (NP,NALL,LDZ,AP_single,W_reverse,Z_reverse)

      deallocate(AP_single)

!-- Correct for close separation of eigenvalues (tail of sequence)

       do j=1,NALL
        if(W_reverse(j).lt.W_bound) then
         W_reverse(j)=0.0d0
          end if
            W_scaled(j)=W_reverse(j)*NALL/mtot
           end do
!         write(*,*) "------------ "
!         write(*,*) " EIGENVALUES "
!         write(*,*) "------------ "
!          DO K=1,NALL
!          write(*,*) " W_scaled(",K,")=",W_scaled(K)
!          END DO

!-- write out the corrected eigenvalues of A --------------
      CLOSE(i_evdA)
       OPEN(i_evdA,file='evd_A',status='unknown',form='formatted')
        do j=1,NALL
!      write(i_evdA,710) W_scaled(j)
       write(i_evdA,710) W_reverse(j)
      end do
     CLOSE(i_evdA,status='KEEP')
 710 format(E20.10)

!--- calculate the TRace of A ----------------
          trace=0.0
          Rodgers_inf_cont_scaled=0.0
          Rodgers_inf_cont=0.0
          entropy=0.0
           DO K=1,NALL
            trace=trace+W_scaled(k)
    Rodgers_inf_cont_scaled=Rodgers_inf_cont_scaled+W_scaled(k)/(1.+W_scaled(k))
    Rodgers_inf_cont=Rodgers_inf_cont+W_reverse(k)/(1.+W_reverse(k))
    entropy=entropy+0.5*log(1.+W_reverse(k))
             END DO
              write(*,*) "------------ "
             write(*,*) " trace=",trace
!            write(*,*) " Rodgers_inf_cont_scaled=",Rodgers_inf_cont_scaled
             write(*,*) " Rodgers_inf_cont=",Rodgers_inf_cont
             write(*,*) " entropy=",entropy
            write(*,*) "------------ "

!--- generate eigenvalues for the inverse ------------
     do j=1,NALL
      W_reverse(j)=1./dsqrt(1.0D0+W_reverse(j)) 
      W_scaled(j)=1./dsqrt(1.0D0+W_scaled(j))
       end do


          write(*,*) "------------ "
          write(*,*) "INVERSE EIGENVALUES "
          write(*,*) "------------ "
           DO K=1,NALL
!          write(*,*) "INVERSE W(",K,")=",W_scaled(K)
           write(*,*) "INVERSE W(",K,")=",W_reverse(K)
           END DO

!-- calculate the TRace of (I+A) inverse -----------
      trace_inv=0.0
       do j=1,NALL
         trace_inv=trace_inv+W_reverse(j)
!        trace_inv=trace_inv+W_scaled(j)
           end do
         write(*,*) "------------ "
       write(*,*) " trace_inv=",trace_inv
      write(*,*) "------------ "

      obsinf_cont=trace_inv/NENS
      write(*,*) "obsinf_cont=",obsinf_cont
      rewind 40
       if(obsinf_cont.gt.0.99) then
      write(40,*) "skip"
       else
      write(40,*) "run"
       endif

      CLOSE(54)
       OPEN(54,file='trace',status='unknown',form='formatted')
       write(54,*) Rodgers_inf_cont,trace_inv
      CLOSE(54)
!-----------------------------------------------------
     allocate(W_single(1:NALL))
     allocate(Z_single(1:LDZ,1:NALL))

     W_single(:)=W_reverse(:)
     Z_single(:,:)=Z_reverse(:,:)

!--- write the eigenvalues of (I+A)**(-1/2) --------------
       CLOSE(i_evdIPA)
        OPEN(i_evdIPA,file='evd_IPA',status='unknown',form='formatted')
         do j=1,NALL
!       write(i_evdIPA,710) W_scaled(j)
        write(i_evdIPA,710) W_reverse(j)
       end do
      CLOSE(i_evdIPA,status='KEEP')

     deallocate(W_reverse)
    deallocate(Z_reverse)

!-- write the file 'A_eigenfile' (needed for post-processing)
      OPEN(ia_eigen,file='A_eigenfile',status='unknown',form='unformatted')
        WRITE(ia_eigen) NALL
         WRITE(ia_eigen) (W_single(j),j=1,NALL)
           WRITE(ia_eigen) ((Z_single(i,j),i=1,NALL),j=1,NALL)
      CLOSE(ia_eigen,status='KEEP')

!--- calculate the matrix (I+A)**(-1/2) ----
!---- Linux PC Cluster (LAPACK, BLAS) ---------------------

      allocate(d_out(1:LDZ,1:NALL))

      CALL CALCULATE_MATMUL (NALL,LDZ,W_single,Z_single,d_out)

      deallocate(W_single)
      deallocate(Z_single)
      deallocate(W_scaled)

!--- write out the matrix (I+A)**(-1/2) ----
!
!  (since it has small dimensions, 
!   more convenient for subsequent use in full mode (non-packed))

        write(*,*) "d_out min,max =",minval(d_out),maxval(d_out)

        OPEN(ipa_out,file='sqrtIPA_inv',status='unknown',form='unformatted')
         WRITE(ipa_out) ((d_out(i,j),i=1,LDZ),j=1,NALL)
          CLOSE(ipa_out,status='KEEP')

      deallocate(d_out)

      STOP
      END program A_matrix_eigen