module globals
  use mpi ! uncomment for LP64 mpi library
implicit none

  !include 'mpif.h'   ! Uncomment for ILP64 intel mpi library

  integer, parameter :: I4B = selected_int_kind( 8) !4-byte integer
  integer, parameter :: I8B = selected_int_kind(12) !8-byte integer

  ! this is the integer type for mpi control values 
  !integer, parameter :: mpi_int_kind = I8B ! uncomment for ILP64 mpi library
  integer, parameter :: mpi_int_kind = I4B  ! uncomment for LP64 mpi library

  integer, parameter :: db = kind(1.0d0)
  integer, parameter :: stdout = 6

  contains

  !======================================================================
  ! Create a diagonal matrix with ones on the diagonal
  subroutine create_matrix(m,n,eye,mat)
    implicit none
    integer, intent(in) :: m, n !row,col length
    logical, intent(in) :: eye ! diagonal = 1?
    complex(db), allocatable, intent(out) :: mat(:,:) ! matrix
  
    integer :: i

    allocate(mat(m,n))
    ! set diagonal to 1, 0 otherwise
    mat = cmplx(0._db,0._db,db)
    if ( eye ) then
       do i = 1, min(m,n)
          mat(i,i) = cmplx(1.0_db,0.0_db,db)
       end do
    end if
    return
  
  
  end subroutine create_matrix
  
  !======================================================================

end module globals
!======================================================================
program test
  use globals
  implicit none

  integer(mpi_int_kind) :: rank_, ierr, nproc_, version, subversion
  integer(mpi_int_kind) :: mpi_int_size
  logical :: master
  integer :: rank, nproc
  integer(I4B) :: isize,v_,sv_

  real(db), allocatable :: ss(:)
  complex(db), allocatable :: uu(:,:), vvH(:,:), zz(:,:)


  integer :: lwork, m, n, info
  complex(db), allocatable :: work(:)
  real(db),allocatable :: rwork(:)
  complex(db) :: work_query(1)

  character(len=1) :: jbU, jbVH

  character(len=198) :: mversion

  ! Initialize the MPI environment
  call mpi_init(ierr)
  call mpi_comm_rank(MPI_COMM_WORLD, rank_, ierr)
  call mpi_comm_size(MPI_COMM_WORLD, nproc_, ierr)
  rank = rank_; nproc = nproc_
  master = rank .eq. 0

  if( master ) then
     call mpi_get_version(version,subversion,ierr)
     write(stdout,'(2(A,I1.1))') ' MPI Version : ', version, '.', subversion

     call mpi_get_version(v_,sv_,ierr)
     write(stdout,'(2(A,I1.1))') ' MPI Version : ', v_, '.', sv_
  
     call mpi_type_size(MPI_INTEGER,mpi_int_size,ierr)
     write(stdout,'(A,I)') ' MPI Integer Size : ', mpi_int_size

     call mpi_type_size(MPI_INTEGER,isize,ierr)
     write(stdout,'(A,I)') ' MPI Integer Size : ', isize
     
     call mkl_get_version_string(mversion)
     write(stdout,*) mversion
     write(stdout,*) 'Number of nodes = ', nproc
  endif
  
  !jbU = 'S'; jbVH = 'S'  ! OK
  jbU = 'N'; jbVH = 'S'  ! INCORRECT LWORK
  !jbU = 'S'; jbVH = 'N'  ! INCORRECT LWORK
  !jbU = 'N'; jbVH = 'N'  ! INCORRECT LWORK

  ! create diagonal matrices
  do m = 10,10,1
  do n = 1 ,10,1
     call create_matrix(m,n,.true.,zz)
     call create_matrix(m,min(m,n),.false.,uu)
     call create_matrix(min(m,n),n,.false.,vvH)
     
     allocate(ss(min(m,n)))
     
     write(stdout,*) 'Starting SVD', m, n

     allocate(rwork(5*min(m,n)))
     lwork=-1
     call zgesvd(jbU, jbVH, m,n, zz, max(1,m), ss, uu, &
                 max(1,m), vvH, max(1,min(m,n)), work_query, lwork, rwork, info)
     lwork = max(1,int(work_query(1)))
     !write(stdout,*) rank, lwork
     allocate(work(lwork)) 
     call zgesvd(jbU, jbVH, m,n, zz, max(1,m), ss, uu, &
                 max(1,m), vvH, max(1,min(m,n)), work, lwork, rwork, info)
     deallocate(rwork,work)

     if ( info .ne. 0 ) write(stdout,*) 'INFO = ', info
     !call mpi_barrier(MPI_COMM_WORLD,ierr)
     write(stdout,*) 'Done with SVD'


     deallocate(zz,uu,vvH,ss)
  end do
  end do

  call mpi_finalize(ierr)

end program test

