xxmr2d:out of memory when using pzheev

[wm_atkinson]

I just recently installed ScaLapack on a cluster of Opterons using the Portland Group compilers (pgf90). Everything seems to have compiled fine, but when I write a test code using pzheev to find the eigenvalues/vectors of a complex*16 distributed matrix I get an "xxmr2d:out of memory" error. The same code works fine with pdsyev routine (providing you change complex (8) to real (8) as appropriate). Does anybody have any idea what's going on? I've attached the test code below, if that makes a diference.

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


[Stan Tomov]

In the text of your message you say that the matrix type is complex*16 but in
the code you use complex*8. I changed that and the size of lrwork
to 158 (4*N-2), and the test completes without problems. I guess it
is a bug that rwork(1) didn't give the minimum size required for the
array rwork.
Stan

[wm_atkinson]

In the text of your message you say that the matrix type is complex*16 but in
the code you use complex*8. I changed that and the size of lrwork
to 158 (4*N-2), and the test completes without problems. I guess it
is a bug that rwork(1) didn't give the minimum size required for the
array rwork.
Stan

Hold on. I need to make sure I'm not doing something stupid here. If I check

Code: Select all

program testit
  implicit none
  complex (8) :: c1
  complex*16 :: c2

  write(6,*) kind(c1)
  write(6,*) kind(c2)
end program testit


then I get

Code: Select all

            8
            8


If I am not mistaken, then complex*16 is the same as complex (8).

I think the kind of error you describe with rwork is the same one I got when I spent a (very frustrating) day using the wrong precision subroutine. To summarise, I think the test code I submitted is fine, except that I forgot to define real (8) :: r1.

-Bill

[Stan Tomov]

You are right, the precision is O.K. I didn't see your definition of float
and thought this is single precision specific for your compiler.
I also had to define r1. With these changes I still need to put lrwork
at least 158, otherwise I get segmentation fault (pzheev returns 82).
I used the GNU compilers (g77/gfortran) to run the test.

[wm_atkinson]

Remarkable. It executes! Ok so there seems to be something wrong with the workspace query in pzheev. At least I can be somewhat confident that the problem doesn't lie in how I built ScaLapack.

Just out of curiosity, did you try running with pdsyev instead of pzheev? That works just fine for me.

Thanks for the help.

[Stan Tomov]

I had tried the test with pdsyev and
it worked fine for me too.
Stan