program min_test

implicit none
include 'mpif.h'
integer :: na,nblk,np_rows,np_cols,na_rows,na_cols,myid,nprocs,my_prow,my_pcol,mpierr,my_blacs_ctxt,sc_desc(9),info
complex, allocatable :: a(:,:),mat(:,:)

integer, external    :: numroc
!-------------------------------------------------------------------------------


   ! default parameters
   na = 8
   nblk = 2

   call mpi_init(mpierr)
   call mpi_comm_rank(mpi_comm_world,myid,mpierr)
   call mpi_comm_size(mpi_comm_world,nprocs,mpierr)

   do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
     if(mod(nprocs,np_cols) == 0 ) exit
   enddo
   ! at the end of the above loop, nprocs is always divisible by np_cols
   np_rows = nprocs/np_cols

   ! initialise BLACS
   my_blacs_ctxt = mpi_comm_world
   call BLACS_Gridinit(my_blacs_ctxt, 'C', np_rows, np_cols)
   call BLACS_Gridinfo(my_blacs_ctxt, np_rows, np_cols, my_prow, my_pcol)

   ! determine the neccessary size of the distributed matrices,
   ! we use the scalapack tools routine NUMROC

   na_rows = numroc(na, nblk, my_prow, 0, np_rows)
   na_cols = numroc(na, nblk, my_pcol, 0, np_cols)

   call descinit(sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info)

   if (myid==0) then
     print '(a)',' |- BLACS and DESC initialized '
   end if

! the total matrix mat(na,na) is allocated
   allocate(mat(na,na))
!  and generate (as a generic random matrix)
   call fill_matrix(mat,na)
   call mpi_barrier(mpi_comm_world,mpierr)
   call mpi_bcast(mat,na*na,mpi_complex,0,mpi_comm_world,mpierr)

! the matrix mat is distributed to the local matrices a(na_rows,na_cols) 
!  using the block-cyclic convention
   allocate(a(na_rows,na_cols))
   call distrib_mat(mat,na,a,na_rows,na_cols,np_cols,np_rows,my_pcol,my_prow,nblk)
                    
   print '(a,3i5)',' |- Block cyclic mat with dim (and proc) ',na_rows,na_cols,myid

! let's diagonalize mat sequentially
   if (myid==0) then
    call diago_seq(mat,na)
    print *
    print '(a)',' |- Reference matrix diagonalized sequentially, using lapack'
    print *
   endif
   call mpi_barrier(mpi_comm_world,mpierr)
! Now we do the diago using scalapack
   call diago_scalapack(a,na_rows,na_cols,sc_desc,na)
   if (myid==0) then
    print '(a)',' |- diagonalised using SCALAPACK ... '
    print *
   end if
   call blacs_gridexit(my_blacs_ctxt)
   call mpi_finalize(mpierr)
end program min_test

subroutine diago_scalapack(matrix,n,m,desca,na)
use inter
implicit none
integer :: n,m,na,info,ierr,rang,nprocs, lwork, liwork, idum,desca(9)
complex :: matrix(n,m), zdum(1,1)
complex, allocatable :: eigenvec(:,:),work(:)
real, allocatable :: wr(:),wi(:),rwork(:),scale(:)
real :: dum
include 'mpif.h'

 call mpi_comm_size(mpi_comm_world,nprocs,ierr)
 call mpi_comm_rank(mpi_comm_world,rang,ierr)

! allocate eigenvectors and eigenvalues (real and imaginary part)
 allocate(eigenvec(n,m))
 allocate(wr(na))
 allocate(wi(na))
! first call to get the ideal minimum work size 
 lwork = -1
 allocate(work(1))
 allocate(rwork(na))
 allocate(scale(na))
 call pcgeevx('N', 'N', 'V', 'N', na, matrix, desca, wr,wi, zdum,desca, eigenvec, desca, &
              idum, liwork,scale, dum, rwork, dum, work, lwork, info)
 lwork=int(abs(work(1)))
 write(*,*)'info ',info,rang
 call mpi_barrier(mpi_comm_world,ierr)
 deallocate(work)
! second call to get the eigenvectors and eigenvalues 
 allocate(work(lwork))
 call pcgeevx('B', 'N', 'V', 'N', na, matrix, desca, wr, wi, zdum,desca, eigenvec, desca, &
              idum, liwork, scale, dum, rwork, dum, work, lwork, info)
 write(*,*)'info ',info,rang
end subroutine diago_scalapack


subroutine fill_matrix(mat,na)
implicit none
integer :: na,i,j
real :: re,im
complex :: mat(na,na)
 
! this generates a random quasi-hermitian matrix
do i=1,na/2
 do j=1,na/2
  call random_number(re)
  call random_number(im)
  mat(i,j)=cmplx(re*2.-1.,im*2.-1.)
 enddo
enddo
do i=1,na/2
 do j=1,i-1
  mat(i,j)=conjg(mat(j,i))
 enddo
 mat(i,i)=cmplx(real(mat(i,i)),0.)
enddo
! second block (C = symmetric)
do i=1,na/2
 do j=na/2+i,na
  call random_number(re)
  call random_number(im)
  mat(i,j)=cmplx(re*2.-1.,im*2.-1.)
 enddo
enddo
do i=1,na/2
 do j=na/2+1,na/2+i-1
  mat(i,j)=mat(j-na/2,i+na/2)
 enddo
 mat(i,i+na/2)=cmplx(real(mat(i,i+na/2)),0.)
enddo
! third block (-C* = symmetric)
do i=na/2+1,na
 do j=1,na/2
  mat(i,j)=-conjg(mat(j,i))
 enddo
enddo
! fourth block (A=-R*)
do i=na/2+1,na
 do j=na/2+1,na
  mat(i,j)=-mat(j-na/2,i-na/2)
 enddo
enddo
end subroutine fill_matrix

subroutine distrib_mat(mat,na,a,na_rows,na_cols,np_cols,np_rows,my_pcol,my_prow,nblk)
implicit none
integer :: na,na_cols,na_rows,nblk,np_rows,np_cols,my_pcol,my_prow
integer :: i,j,il,jl
complex :: a(na_rows,na_cols),mat(na,na)
integer, external :: INDXL2G
do jl = 1 , na_cols
 j = INDXL2G(jl,nblk,my_pcol,0,np_cols)
 do il = 1 , na_rows
  i = INDXL2G(il,nblk,my_prow,0,np_rows)
  a(il,jl)= mat(i,j)
 end do
end do
end subroutine distrib_mat

subroutine diago_seq(mat,na)
implicit none
integer :: na,i,j,lwork,info
complex :: mat(na,na), vl(1,1)
complex, allocatable :: work(:),eigenvec(:,:), eigenval(:)
real, allocatable :: rwork(:)

 lwork = max(1,2*na)
 allocate(work(lwork))
 allocate(rwork(2*na))
 allocate(eigenval(na))
 allocate(eigenvec(na,na))

 call cgeev('N','V',na,mat,na,eigenval,vl,1,eigenvec,na,work,lwork,rwork,info)
 deallocate(work,rwork)
 
 do j=1,na
  write(50, *) 'eigenval ',j,'=',real(eigenval(j)),aimag(eigenval(j))
  do i=1,na
   write(50, fmt = 9999 ) 'z',i,j,real(eigenvec(i,j)),aimag(eigenvec(i,j))
  enddo
 enddo

 9999 format(a,'(',i6,',',i6,')=',e16.8, '+i*(',e16.8, ')')
end subroutine diago_seq

