- Code: Select all
!
! TEST CODE TO COMPARE PDSYEV AND PZHEEV
!
program test
use mpi
implicit none
integer, parameter :: float=selected_real_kind(6,70)
integer :: i1,i2
integer, parameter :: N = 40
real (float) ::W(N)
!%%% Uncomment for PZHEEV
COMPLEX (float) :: A(N,N),Z(N,N)
COMPLEX (float), allocatable :: a0(:,:),z0(:,:)
COMPLEX (float), allocatable :: work(:)
!%%% Uncomment for PCHEEV
!REAL (float) :: A(N,N),Z(N,N)
!REAL (float), allocatable :: a0(:,:),z0(:,:)
!REAL (float), allocatable :: work(:)
real (float), allocatable :: w0(:)
integer :: lwork,lrwork
real (float), allocatable :: rwork(:)
! blacs
integer :: ctxt,MyRow,MyCol,nproc
integer, parameter :: Prow=2,Pcol=1
integer, parameter :: Brow=16,Bcol=16
integer, parameter :: dlen_ = 9
integer :: desca0(dlen_),descz0(dlen_),descw0(dlen_),LDA,LDB,tmp
integer :: numroc
! network topology
character (1) :: TOP = ' '
! mpi
integer :: ierr,myrank
external numroc
call SL_init(ctxt,Prow,Pcol)
call mpi_comm_size(mpi_comm_world,nproc,ierr)
if (nproc /= Prow*Pcol) then
write(6,*) "Error! Must have ",Prow*Pcol," processors"
stop
end if
call mpi_comm_rank(mpi_comm_world,myrank,ierr)
call blacs_gridinfo(ctxt,Prow,Pcol,MyRow,MyCol)
LDA = numroc(N,Brow,myrow,0,Prow)
LDB = numroc(N,Bcol,mycol,0,Pcol)
write(6,*) MyRow,MyCol,myrank,LDA,LDB
allocate(a0(lda,ldb),z0(lda,ldb),w0(N))
call descinit(desca0,N,N,Brow,Bcol,0,0,ctxt,LDA,ierr )
call descinit(descz0,N,N,Brow,Bcol,0,0,ctxt,LDA,ierr )
if (myrow+mycol==0) then
! generate a symmetric matrix with random elements
call random_seed()
do i1 = 1,N
do i2 = i1,N
call random_number(r1)
A(i2,i1) = r1
A(i1,i2)=A(i2,i1)
end do
end do
! distribute the matrix
!%%% Uncomment for PDSYEV:
!call DGEBS2D(ctxt,'A',TOP,N,N,A,N)
!%%% Uncomment for PZHEEV:
call ZGEBS2D(ctxt,'A',TOP,N,N,A,N)
else
! receive the matrix
!%%% Uncomment for PDSYEV:
!call DGEBR2D(ctxt,'A',TOP,N,N,A,N,0,0)
!%%% Uncomment for PZHEEV:
call ZGEBR2D(ctxt,'A',TOP,N,N,A,N,0,0)
end if
do i1 = 1,N
do i2 = 1,N
!%%% Uncomment for PDSYEV:
!call pdelset(a0,i1,i2,desca0,A(i1,i2))
!%%% Uncomment for PZHEEV:
call pzelset(a0,i1,i2,desca0,A(i1,i2))
end do
end do
lwork = -1
lrwork = -1
allocate(work(1),rwork(1))
!%%% Uncomment for PDSYEV:
!call pdsyev('V','U',N,a0,1,1,desca0,w0,z0,1,1,descz0,WORK,LWORK,&
! rwork,lrwork,ierr)
!%%% Uncomment for PZHEEV:
call pzheev('V','U',N,a0,1,1,desca0,w0,z0,1,1,descz0,WORK,LWORK,&
rwork,lrwork,ierr)
write(6,*) "done pzheev 1"
lwork = int(abs(work(1)))
lrwork = int(rwork(1)+1)
deallocate(work,rwork)
allocate(work(lwork),rwork(lrwork))
!call pdsyev('V','U',N,a0,1,1,desca0,w0,z0,1,1,descz0,WORK,LWORK,&
! rwork,lrwork,ierr)
!%%% Uncomment for PZHEEV:
call pzheev('V','U',N,a0,1,1,desca0,w0,z0,1,1,descz0,WORK,LWORK,&
rwork,lrwork,ierr)
write(6,*) "done pzheev 2"
deallocate(work,rwork)
deallocate(a0,z0,w0)
call blacs_gridexit(ctxt)
call blacs_exit(0)
end program test