module magma_opts

! DSYEVD Constants as defined in magmamic-1.1.0/include 10/07/2014

integer, parameter :: MagmaUpper = 121
integer, parameter :: MagmaLower = 122
integer, parameter :: MagmaNoVec = 301
integer, parameter :: MagmaVec = 302
integer, parameter :: MagmaErr = -100

contains

 function magma_jobz(jobz)
   ! Convert character to integer for Magmamic DSYEVD call
   integer :: magma_jobz
   character, intent(in) :: jobz

   select case(jobz)
     case('V')
       magma_jobz = MagmaVec
     case('v')
       magma_jobz = MagmaVec
     case('N')
       magma_jobz = MagmaNoVec
     case('n')
       magma_jobz = MagmaNoVec
     case default
       magma_jobz = MagmaErr
   end select

end function magma_jobz

 function magma_uplo(uplo)
   ! Convert character to integer for Magmamic DSYEVD call
   integer :: magma_uplo
   character, intent(in) :: uplo

   select case(uplo)
     case('L')
       magma_uplo = MagmaLower
     case('l')
       magma_uplo = MagmaLower
     case('U')
       magma_uplo = MagmaUpper
     case('u')
       magma_uplo = MagmaUpper
     case default
       magma_uplo = MagmaErr
   end select

end function magma_uplo

end module magma_opts                   

module dist_mat
! Form diestributed sector Hamiltonian matrix and diagonalize
! Time-stamp: "2003-09-23 17:08:23 cjn"
  use io_units, only: fo
  use precisn, only: wp
!  use blacs, only: ctxt, myrow, mycol, p, q, nblock, mynumrows,&
!       mynumcols, p_error, dlen_, nb_, csrc_
  use error_out, only: error_check
  use magma
  use magma_opts
  implicit none

  integer, save               :: ny           ! a dimension
  real(wp), allocatable, save :: a(:)         ! distributed matrix
  real(wp), allocatable, save :: z(:)         ! eigenvectors
  real(wp), allocatable, save :: evals(:)     ! eigenvalues
!  integer, save               :: desca(dlen_) ! Hamiltonian descriptor
!  integer, save               :: descz(dlen_) ! Eigenvector descriptor

  private
  public A_fill, sh_diag, kill_az
  public evals, z

contains

  subroutine A_fill (n)
! fill lower triangle of symmetric matrix, A
    use sec_ham, only: h_el
    integer, intent(in)   :: n       ! matrix dimension
    integer               :: j, jj, jloc, jglob, i, ii, iloc, iglob
    integer               :: status, info, numroc
    integer :: ind

    ny = n

! Compute local dimensions of matrix A
!    mynumrows = NUMROC(n, nblock, myrow, 0, p)
!    mynumcols = NUMROC(n, nblock, mycol, 0, q)

! Allocate local part of A
    allocate (a(ny*ny), z(ny*ny), stat=status)
    call error_check (status, ' Allocation Error, A_fill')

    a = 0.0_wp        ! initialize A

! Form array descriptor od destributed matrix A
!    call descinit (desca, n, n, nblock, nblock, 0, 0, ctxt, mynumrows,&
!         info)
!    descz = desca

!    do j = 1, mynumcols, nblock                   ! local column blocks
!       do jj = 1, MIN(nblock, mynumcols-j+1)      ! all columns
!          jloc = j - 1 + jj                       ! local column index
!          jglob = (j-1) * q + mycol * nblock + jj ! global column index

!          do i = 1, mynumrows, nblock
!             do ii = 1, MIN(nblock, mynumrows-i+1)
!                iloc = i - 1 + ii                       ! local index
!                iglob = (i-1) * p + myrow * nblock + ii ! global index
!                if (iglob > jglob) cycle
                
!                a(iloc+(jloc-1)*mynumrows) = h_el(iglob, jglob)

!             end do
!          end do
!       end do
!    end do

! sequential version - full matrix
  ind = 0
  do j = 1, ny
    do i = 1, j
      ind = ind + 1
      a(ind) = h_el(i,j)
    end do
  end do

  end subroutine A_fill

  subroutine sh_diag
! diagonalize the matrix a
    real(wp), allocatable, save :: r(:)
    real(wp), allocatable        :: work(:) 
    integer, allocatable         :: iwork(:)
    integer                      :: lwork, liwork
    real(wp)                     :: alpha, beta
    character(len=1)             :: jobz, uplo, MagmaRange
    logical                      :: ismycol
    real(wp)                     :: dnorm2
    integer                      :: status, ii, nstrt, trilwmin
    integer                      :: np, nq, lda, iproc, info, numroc
    integer                      :: indxg2p, lwork1, ip, iq, imyrow, imycol
    integer :: rows, cols
    real(wp)                     :: t0, t1
    integer                      :: c0, c1, cr
    real(wp), external :: dnrm2
! Magma settings
    integer :: ijobz = 999 ! default setting
    integer :: iuplo = 999 ! default setting
    magma_devptr_t :: dA, dEvals, dWork, dIwork
    integer :: nGPU=8
    integer :: sizeof_double=8
    integer :: sizeof_int=4
    integer :: queue = -66
    real(wp):: vl, vu
    integer :: il, iu
    integer, pointer ::  mout(:)

! initialize eigenvector distributed array, z.
     alpha = 0.0_wp
     beta = 1.0_wp
!     call pdlaset ('all', ny, ny, alpha, beta, z, 1, 1, descz)
!     np = mynumrows
!     nq = mynumcols5

! perform diagonalization:
     jobz = 'V'
     uplo = 'U'
     liwork = 3 + 5*ny
     lwork = 1 + 6*ny + 2*ny*ny

     vl=0.0
     vu=0.0
     il=1
     iu=ny
     MagmaRange = 'A'
     !ijobz = magma_jobz(jobz)
     !iuplo = magma_uplo(uplo)

     !call cublas_alloc( lwork,    sizeof_int, dWork )
     !call cublas_alloc( liwork,   sizeof_int, dIwork )
     !call cublas_alloc( ny,       sizeof_double, dEvals )
     !call cublas_alloc( ny*ny,    sizeof_double, dA)
     !allocate (r(ny*ny) , stat=status)
     !call error_check (status, ' Allocation Error, r')
    
     !call CUBLAS_SET_VECTOR(lwork , sizeof_int    ,   work  , 1, dWork , 1)
     !call CUBLAS_SET_VECTOR(liwork, sizeof_int    ,   iwork , 1, dIwork, 1)
     !call CUBLAS_SET_VECTOR(ny    , sizeof_double ,   evals , 1, dEvals, 1)
     !call cublas_set_matrix(ny, ny, sizeof_double , a , ny, dA, ny)

     allocate (work(lwork), iwork(liwork), evals(ny),  stat=status)
     call error_check (status, 'sh_diag: allocation error')

     write (fo,'(/,a,/)') 'start of DSYEVD'
     call cpu_time (t0)
     call system_clock (count=c0)

    !call magmaf_dsyevd(jobz, uplo, ny, a, ny, evals, work, lwork, iwork,liwork, info)
    !call magmaf_dsyevd_gpu(jobz, uplo, ny, dA, ny, evals, r, ny,  work, lwork, iwork,liwork, info)
    !call magmaf_dsyevdx_m(nGPU, jobz, uplo, ny, a, ny, evals, work, lwork, iwork,liwork, info)
    call magmaf_dsyevdx( jobz, MagmaRange,  uplo, ny, a, ny, vl, vu, il, iu, mout, evals, work, lwork, iwork,liwork, info)
    !call dsyevd (jobz, uplo, ny, a, ny, evals, work, lwork, iwork, liwork, info)

     write (fo,'(/,a,/)') 'end of DSYEVD'
     call cpu_time (t1)
     write (fo,'(a,f16.4,a)') 'DSYEVD CPU time     = ', t1 - t0, ' secs'
     call system_clock (count=c1, count_rate=cr)
     write (fo,'(a,f16.4,a)') 'DSYEVD Elapsed time = ', REAL(c1-c0,wp) / &
             REAL(cr,wp), ' secs'

     call error_check (info, 'sh_diag: pdsyevd error')

     deallocate (work, iwork, stat=status)
     call error_check (status, 'sh_diag: deallocation error')

     !call cublas_get_vector (ny,  sizeof_double, dEvals, 1, evals, 1)
    
     print *, "DEBUG : Eigen value = ",evals(1)
     print *, "DEBUG : Eigen value = ",evals(ny)
  
     z = a
     evecs: do ii = 1, ny ! normalize eigenvectors
! indxg2p computes process coord which posseses entry of a
! distributed matrix specified by a global index INDXGLOB.
!        iproc = indxg2p(ii, descz(nb_), mycol, descz(csrc_), q)
!        ismycol = (iproc == mycol)
!        if (ismycol) then
!           call pdnrm2 (ny, dnorm2, z, 1, ii, descz, 1)
            
           dnorm2 = dnrm2(ny,z,1)    
           alpha = 1.0_wp / REAL(dnorm2,wp)
           call dscal(ny, alpha, z, 1)
!        end if
     end do evecs
     !call cublas_free( dA )
     !deallocate (r, stat=status)
     !call error_check (status, 'r : deallocation error')
   end subroutine sh_diag

   subroutine kill_az
     integer       :: status
     deallocate (a, z, evals, stat=status)
     call error_check (status, 'kill_az: deallocation error')
   end subroutine kill_az

end module dist_mat

