PCHEEVX and PCGEMR2d

[sipahi]

Hi,
I have been trying to use pcheevx and pcgemr2d together.
I have a code that needs to diagonalize a series of matrices and use their eigenvalues together to perform extensive multiplications.
The time spent in diagonalizations compares to the multiplications one.
In previous codes, I used to distribute a set of result eigenvectors in different processors and gather them selectively as i needed.
As the matrices sizes have grown, I started to use scalapack to find the vectors. The next step would be redistribute the vectors in order to store each vector in one processor to facilitate the gather procedure.
I have been trying to use pcheevx and pcgemr2d in the same routine, but I cannot succeed. Each time I start using the two routines alltogether I get a SIGSEGV message on execution.

I have been playing with some examples that I found in the net and I produced a code that has the routines I need to work together, and even with a very small code I cannot make running correctly to the end.

I would appreciate any help on the code attached below.
Thank you in advance,
Guilherme Sipahi
----------------------------------------------------------------------------------
implicit none

integer, parameter :: m = 1000, n = m, mb = 8, nb = mb
integer, parameter :: rsrc = 0, csrc = 0, dlen_ = 9

integer, parameter :: mm = 10

integer :: ictxt, rank, tproc, icroot
integer :: prow, pcol, nprow, npcol, myrow, mycol
integer :: mxlocr, mxlocc, mxlld, mxlld_root
integer :: i,j,k,info,IOFILE

integer, external :: blacs_pnum, numroc

complex, allocatable :: a(:,:), b(:,:), c(:,:),z(:,:)
real, allocatable :: w(:)
integer :: desca(dlen_), descb(dlen_), descdiago(dlen_)
integer :: IERR

complex, allocatable :: work(:)
real, allocatable :: rwork(:)
integer, allocatable :: iwork(:)
integer, allocatable :: ifail(:)
integer, allocatable :: iclustr(:)
real, allocatable :: gap(:)

integer :: lwork,lrwork,liwork
integer :: il,iu
integer :: nz
real :: vl,vu
real :: abstol
real :: orfac = -1.0
integer :: nn

real :: pslamch
external :: pslamch

include 'mpif.h'

call MPI_INIT(IERR)

call blacs_pinfo(rank, tproc)
call blacs_get(-1,0,ictxt)
nprow = int(sqrt(real(tproc+0.0001)))
npcol = tproc/nprow

call blacs_gridinit(ictxt,'r',nprow, npcol)
call blacs_gridinfo(ictxt,nprow,npcol,myrow,mycol)
print *, 'rank',rank, 'myrow', myrow, 'mycol',mycol
!print *, '(nprow,npcol)', '(',nprow,npcol,')'

mxlocr = numroc(m,mb,myrow, rsrc,nprow)
mxlocc = numroc(n,nb, mycol, csrc, npcol)
mxlld = max(1,mxlocr)

print *,'mxlocr=', mxlocr, 'mxlocc=', mxlocc, 'mxlld=',mxlld

allocate(a(mxlld,mxlocc))
allocate(z(mxlld,mxlocc))
allocate(w(m))
allocate(b(M,N), c(M,mm))
print *,'after allocate'

call descinit(desca,m,n,mb,nb,rsrc,csrc,ictxt,mxlld,info)
IF (INFO .NE. 0) print *, 'info NE zero', info

!initialize the matrix in one processor
IF (myrow .eq. 0 .and. mycol .eq. 0) then
k = 0
do j=1, N
do i=1, M
k = k+1
b(i,j) = k
enddo
enddo
!write(*,202) ((b(i,j),j=1,n),i=1,m)
do i=1,m
write(8,202) b(i,1:mm)
enddo

ENDIF
202 format(100e15.6e2)

a=0.0d0
call blacs_barrier(ictxt, 'A')

call blacs_get(-1,0,icroot)
call blacs_gridinit(icroot,'R',1,1)

IF (myrow .eq. 0 .and. mycol .eq. 0) then
mxlld_root = numroc(M,M,myrow,0,nprow)
print *, 'myrow =', myrow
print *, 'nprow =', nprow
print *, 'mxlld_root=', mxlld_root
call descinit(descb,m,n,m,n,rsrc,csrc,icroot,mxlld_root,info)
IF (INFO .NE. 0) print *, 'info NE zero', info
ELSE
descb(1:9) = 0
descb(2) = -1
ENDIF

!distribute the matrix to the different processors

call PCGEMR2D( M, N, b, 1, 1, descb, A, 1, 1, desca, desca( 2 ) )

call blacs_barrier(ictxt, 'A')

ABSTOL = PSLAMCH(ictxt, 'U')

vl = 0.0
vu = 0.0
il = 203
iu = 298

ALLOCATE(IFAIL(m),ICLUSTR(2*nprow*npcol),GAP(nprow*npcol))

LWORK = -1
LRWORK = -1
LIWORK = -1
ALLOCATE(WORK(3))
ALLOCATE(RWORK(1))
ALLOCATE(IWORK(1))

call PCHEEVX( 'V', 'I', 'U', M, a, 1, 1, desca, VL, VU, IL, IU, ABSTOL, nn, NZ, w, ORFAC, z, 1, 1, desca,&
WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO )

LWORK = INT(ABS(WORK(1)))
DEALLOCATE(WORK)

LRWORK = INT(ABS(RWORK(1)))
DEALLOCATE(RWORK)

LIWORK = INT(ABS(IWORK(1)))
DEALLOCATE(IWORK)

call blacs_barrier(ictxt, 'A')

call PCHEEVX( 'V', 'A', 'U', M, a, 1, 1, descdiago, VL, VU, IL, IU, ABSTOL, nn, NZ, w, ORFAC, z, 1, 1, descdiago,&
WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO )

call blacs_barrier(ictxt, 'A')

c=0.0d0

call PCGEMR2D( M, 10, A, 1, 1, desca, c, 1, 1, descb, desca( 2 ) )

IF (myrow .eq. 0 .and. mycol .eq. 0) then
call blacs_gridexit(icroot)
ENDIF

IOFILE = RANK+10
write(IOFILE,*)((a(i,j),j=1,mxlocc),i=1,mxlocr)

201 format(4f4.0)


IF (myrow .eq. 0 .and. mycol .eq. 0) then
do i=1,m
write(9,202) c(i,:)
enddo
ENDIF

203 format(100e15.6e2)

call blacs_gridexit(ictxt)
call blacs_exit(0)

end