!  test_cluster_solver_v1.f90 
!
!  FUNCTIONS:
!  test_cluster_solver - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: test_cluster_solver
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program test_cluster_solver
    Use mkl_cluster_sparse_solver
!***--------------------------------------------------------------------
    implicit none
    include 'mpif.h'
!***--------------------------------------------------------------------
!*  Variables
    integer, parameter:: dp = kind(1.0D0)
    integer(4):: mpi_stat, rank, is_distr, ierr
    integer(4):: MKL_COMM
    Character(128):: fn_Ainfo,fn_AMat,fn_BVecs,fn_XVecs
    Integer:: lfn_Ainfo,lfn_AMat,lfn_BVecs,lfn_XVecs,fid
!***--------------------------------------------------------------------
!.. Internal solver memory pointer for 64-bit architectures    
    TYPE(MKL_CLUSTER_SPARSE_SOLVER_HANDLE):: pt(64) 
    Integer:: iparm(64)
    INTEGER:: maxfct, mnum, mtype, phase,error, msglvl
!    Integer, Save:: kset = 1
    Integer:: MatType = 6  !,n_refine = 9, ncores=1 
    Logical:: is_init = .FALSE., is_set_data=.FALSE., is_dbg=.FALSE.
    
    INTEGER,save:: num_equs, nrhs, num_nz
    Integer:: idum(1)
    Complex(kind=dp):: zdum(1)
    Integer, Allocatable :: ia(:)
    Integer, Allocatable :: ja(:)
    Complex(kind=dp), Allocatable :: aij(:)
    Complex(kind=dp), Allocatable :: b(:,:)
!..    Complex(kind=dp), allocatable :: bs(:,:)
    Complex(kind=dp), allocatable :: x(:,:)
!..    Complex(kind=dp), allocatable :: r(:,:)
!***--------------------------------------------------------------------
    Integer:: kr,kc
    Integer:: ip
    Real(kind=dp):: aij_r,aij_i
    Real(kind=dp), Allocatable:: xb_r(:),xb_i(:)
    
!***--------------------------------------------------------------------
!*  Body of test_cluster_solver
    fn_Ainfo = 'amat_1000_csr.txt'
    fn_AMat  = 'amat_1000_csr.txt'
    fn_BVecs = 'bvec_1000.txt'
    fn_XVecs = 'solu_1000.txt'
    
    lfn_Ainfo = len_trim(fn_Ainfo)
    lfn_AMat  = len_trim(fn_AMat)
    lfn_BVecs = len_trim(fn_BVecs)
    lfn_XVecs = len_trim(fn_XVecs)
!.. Initialize task data on master MPI process
    MKL_COMM = MPI_COMM_WORLD
    call mpi_init(mpi_stat)
    call mpi_comm_rank(MKL_COMM, rank, mpi_stat)
!    Write(*,'(/3x,A,I5,3x,A,I5)')'The Rank=',rank,' Call mpi_comm_rank(...), error code =', mpi_stat
!***--------------------------------------------------------------------
    is_distr = 0
    If (is_dbg) Write(*,'(/3x,A,I5)')'The Machine Rank=',rank
!*----------------------------------------------------------------------
!   Set up Cluster Sparse Solver control parameters
    iparm(1:64) = 0
    pt(1:64)%dummy = 0
    
    If (is_distr==0) Then   ! [A],{B} input from master node
       If (rank>=0)  Then   ! Master node 
          maxfct = 1
          mnum = 1
    
         iparm(1) = 1 ! =1 no solver default
         iparm(2) = 2 ! =2 fill-in reduce ordering from METIS
         iparm(3) = 0 ! reverved, set to 0
         iparm(4) = 0 ! reverved, set to 0
         iparm(5) = 0 ! user permutation: = 0 user permutation in perm(:) is ignored
         iparm(6) = 1 ! write solution to b: =1 solution on the first n compoments of b
         iparm(8) = 9 ! iterative refinement step: numbers of iterative refinement steps
         iparm(10) = 8 ! pivoting perturbation: perturbe the pivot elements with 1E-10
         iparm(11) = 0 ! scaling vectors: =1 scaling
         iparm(13) = 0 ! improved accuracy using (non-)symmetric weighted matchings.
         iparm(18) = 0 ! Output: number of nonzeros in the factor LU
         iparm(19) = 0 ! reverved, set to 0
         iparm(21) = 0 ! Pivoting for symmetric indefinte matrices.
                                ! Bunch and Kaufman pivoting is available for matrices: 
                                ! mtype=-2, mtype=-4, or mtype=6
         iparm(24) = 0 ! reverved, set to 0
         iparm(25) = 0 ! reverved, set to 0
         iparm(27) = 1 ! matrix checker: =0 default value, the PARDISO doesn't check the sparse matrix representation.
         iparm(28) = 0 ! single or double precision of PARDISO: =0 the input arrays (matrix a, vectors x and b) as well as all internal arrays are supposed to be in double precision
         iparm(31) = 0 ! Partial solve and computing selected components of the solution vectors
                       ! = 0 disables this option
         iparm(35) = 0 ! One- or zero-based indexing of columns and rows; = 0 ine-based indexing
         iparm(36) = 0 ! Schur complement matrix computation control;
                       ! = 0 Do not compute Schur complement
         iparm(37) = 0 ! Format for matrix storage; = 0 Use CSR format
         iparm(40) = 0 ! Matrix input format;  = 0 centralized input format, master process stores all data of [A] & {B}
         
         iparm(60) = 1 ! in-core PARDISO /OOC PARDISO is used otherwise.
!*---------------------------------------------------------------------------
         error =  0 ! initialize error flag
         msglvl = 1 ! print statistical information
         mtype =  MatType ! complex and symmetric matrix
         
         If (rank==0)  Write(*,'(/3x,A)')'Finish initializing iparm(:) ...'
         is_init = .TRUE.    
       Else
         is_init = .TRUE.    
       End If 
    Else     ! [A],{B} data distributed among the nodes
       Write(*,'(/3x,A)') '*** Warning:'
       Write(*,'( 7x,A)') 'The distributed setting is under developing ...'
       ierr = -100
       stop
    End If
    
    
!***--------------------------------------------------------------------
!*  Read in data on master
    If (is_init) Then
       If (is_distr==0) Then  ! Input on master node 
          If (rank==0)  Then  ! on Master node 
             lfn_Ainfo = len_trim(fn_Ainfo)
             lfn_AMat  = len_trim(fn_AMat)
             lfn_BVecs = len_trim(fn_BVecs)
             lfn_XVecs = len_trim(fn_XVecs)
             
             Open(newunit=fid, file=fn_AMat(1:lfn_AMat), Action='Read')
             Read(fid,*) num_equs, num_nz
             Write(*,'(/3x,A,I5,3x,A,I10)') 'Num. of Equ=',num_equs,'Num. of NZ=', num_nz
             If (Allocated(ia)) deallocate(ia, stat=ierr)
             If (Allocated(ja)) deallocate(ja, stat=ierr)
             If (Allocated(aij)) deallocate(aij, stat=ierr)
             Allocate(ia(num_equs+1),ja(num_nz),aij(num_nz), stat=ierr)
             If (ierr/=0) Then 
                 Write(*,'(/3x,A)') '*** Error ***'
                 Write(*,'( 7x,A)') 'Failed to allocate memory for ia(:),ja(:),aij(:), ...'
                 Close(fid)
                 stop
             End If
             Read(fid,*) (ia(kr),kr=1,num_equs+1)
             Do ip = 1, num_nz 
                 Read(fid,*) ja(ip), aij_r, aij_i
                 aij(ip) = cmplx(aij_r,aij_i)
             End Do
             Close(fid)
             
             Open(newunit=fid, file=fn_BVecs(1:lfn_BVecs), Action='READ')
             Read(fid,*) kr, nrhs
             Write(*,'(/3x,A,I5,3x,A,I10)') 'Num. of Equ=',kr,'Num. of rhs=', nrhs
             If (kr /= num_equs) Then 
                Write(*,'(/3x,A)')'*** Error ***'
                Write(*,'( 7x,A)')'The size of [A] and {B} is NOT matched.'
                Write(*,'( 7x,A)')'The input files are not correct. Return'
                ierr = -101
                Close(fid)
                stop
             End If
             If (Allocated(B)) Deallocate(B, stat=ierr)
             If (Allocated(X)) Deallocate(X, stat=ierr)
             If (Allocated(xb_r)) Deallocate(xb_r, stat=ierr)
             If (Allocated(xb_i)) Deallocate(xb_i, stat=ierr)
             Allocate(B(num_equs,nrhs), X(num_equs,nrhs), xb_r(nrhs),xb_i(nrhs),stat=ierr)
             If (ierr /= 0 ) Then 
                Write(*,'(/3x,A)')'*** Error ***'
                Write(*,'( 7x,A)')'Failed to allocate memory for B(:,:), X(:,:), ...'
                Close(fid)
                stop
             End If
             X(1:num_equs,1:nrhs) = cmplx(0.0d0,0.0d0)
             
             Do kr=1,num_equs 
                Read(fid,*) (xb_r(kc),xb_i(kc),kc=1,nrhs)
                Do kc=1,nrhs 
                   B(kr,kc) = cmplx(xb_r(kc), xb_i(kc))
                End Do
             End Do
             Close(fid)
             
             If (Allocated(xb_r)) Deallocate(xb_r, stat=ierr)
             If (Allocated(xb_i)) Deallocate(xb_i, stat=ierr)
             
             is_set_data = .TRUE.
             Write(*,'(/3x,A)')'Data have been read in. ...'
          Else
             is_set_data = .TRUE.
          End If
       Else  ! Distributed input 
          Write(*,'(/3x,A)') '*** Warning:'
          Write(*,'( 7x,A)') 'The distributed setting is under developing ...'
          ierr = -100
          stop 
       End If
    Else
       Write(*,'(/6x,A)')'PSOLVER Module is not initialized yet.' 
       Write(*,'( 6x,A)')'Do nothing, but return.' 
       Ierr = -100
       stop
    End If
    
!***--------------------------------------------------------------------
!*  synchronize all processes
    call MPI_BARRIER(MKL_COMM,ierr)
    
!***--------------------------------------------------------------------
!*  Solve the problem
    If (is_init.and.is_set_data) Then 
!      Reordering and symbolic factorization
       If (is_dbg) Write(*,'(/3x,A,I5,3x,A)')'The Rank=',rank,' in re-ordering phase ...'
       phase = 11
       call cluster_sparse_solver(pt,maxfct,mnum,mtype,phase,num_equs,aij,ia,ja,idum,nrhs,iparm,msglvl,b,x,MKL_COMM,error)
       if (error.ne.0) then
           if (rank.eq.0) print *, 'ERROR during symbolic factorization: ', error
           ierr = error
           stop 
       End If
       If (rank.eq.0) write(*,'(/3x,A)') 'Reordering completed ... '

!      Numeric factorization
       If (is_dbg) Write(*,'(/3x,A,I5,3x,A)')'The Rank=',rank,' in factorization phase ...'    
       phase = 22
       call cluster_sparse_solver(pt,maxfct,mnum,mtype,phase,num_equs,aij,ia,ja,idum,nrhs,iparm,msglvl,b,x,MKL_COMM,error)
       if (error.ne.0) then
           if (rank.eq.0) print *, 'ERROR during numerical factorization: ', error
           ierr = error
           stop 
       End If
       If (rank.eq.0) write(*,'(/3x,A)') 'Factorization completed ... '
       
!      Forward & backward substitution
       If (is_dbg) Write(*,'(/3x,A,I5,3x,A)')'The Rank=',rank,' in solution phase ...'
       phase = 33
       call cluster_sparse_solver(pt,maxfct,mnum,mtype,phase,num_equs,aij,ia,ja,idum,nrhs,iparm,msglvl,b,x,MKL_COMM,error)
       if (error.ne.0) then
           if (rank.eq.0) print *, 'ERROR during during solution: ', error
           ierr = error
           stop 
       End If
       If (rank.eq.0) write(*,'(/3x,A)') 'Solution completed ... '
       
       If (rank==0) Then ! output solution to file
          lfn_XVecs = len_trim(fn_XVecs)
          Open(newunit=fid, file=fn_XVecs(1:lfn_XVecs), Action='Write')
          Do kr = 1,num_equs 
             Write(fid,'(*(G0,3x,G0))') (b(kr,kc),kc=1,nrhs)
          End Do
          Close(fid)
       End If
    Else 
       Write(*,'(/6x,A)')'PSOLVER Module either is not initialized, or' 
       Write(*,'( 6x,A)')'is not set data.'
       Write(*,'( 6x,A)')'Do nothing, but return.' 
       Ierr = -102
       stop
    End If
    
!***--------------------------------------------------------------------
!*  Close cluster_sparse_solver
    If (is_init) Then
!*---------------------------------------------------------------------------
       phase = -1 ! release internal memory
       CALL cluster_sparse_solver(pt,maxfct,mnum,mtype,phase,kr,zdum,idum,idum, &
                   &idum,kc,iparm,msglvl,zdum,zdum,MKL_COMM,error)
!*---------------------------------------------------------------------------
    Else
       Write(*,'(/6x,A)')'PSOLVER Module is not initialized yet.' 
       Write(*,'( 6x,A)')'Call PSOLVER_CLOSE is not needed.'
       Write(*,'( 6x,A)')'Do nothing, but return.' 
       Ierr = -100
       stop
    End If
!*---------------------------------------------------------------------------
!*  Release allocated data 
    If (Allocated(B)) Deallocate(B, stat=ierr)
    If (Allocated(X)) Deallocate(X, stat=ierr)
    If (Allocated(ia)) deallocate(ia, stat=ierr)
    If (Allocated(ja)) deallocate(ja, stat=ierr)
    If (Allocated(aij)) deallocate(aij, stat=ierr)
    
!*---------------------------------------------------------------------------
    is_init=.FALSE.
    
!***--------------------------------------------------------------------
!*  Close MPI process    
    call mpi_finalize(mpi_stat)
    If (rank==0) Write(*,'(/3x,A,I5,3x,A,I5)')'The Rank=',rank,' Call mpi_finalize(...), error code =', mpi_stat
    
!***--------------------------------------------------------------------
    end program test_cluster_solver
    
    
    
!    integer function mkl_progress( thread_process, step, stage )
!    integer*4 thread_process, step, lstr
!    character*(*) stage
!    lstr = len_trim(stage)
!    Write(*,'(3x,A,I7,3x,A,3x,A,I7)')'Thread:',thread_process,', stage: '//stage(1:lstr), ', step:',step
!    mkl_progress = 0
!    return
!    end function

