- 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