PDLAHQR wrong? PZLHQR also?

[Klystron]

Hi all,

I've encountered strange results quite frequently with PDLAHQR, while PZLAHQR always seemed to do a good job. To find out whetheer there is something wrong with PDLAHQR, or I am making a mistake myself, I wrote a small test program, which computes eigenvalues of a Hessenberg matrix with DGEEV, PDLAHQR, and PZLAHQR, comparing the results. I am pasting this program below. It it meant to run on a single processor only, but in an MPI environment. So the data is cut into blocks formally, but this should have no actual effect. In the program, there is a variable "opt" which should be set to any value between 1 and 5 before compiling. This parameter introduces minor chnages in matrix and block sizes, and matrix elements. In some cases, the three results agree, in others they don't. The fact that in case opt=5, even PZLAHQR doesn't agree with DGEEV, makes me feel particularly uncomfortable, as I am relying heavily on PZLAHQR in my work, where results always seemed to be ok until now.

Could anybody, please, test this program, run it with the differen opt values, and tell me if the strange effects also occur on his/her machine, and tell me (hopefully) that the mistake is in my program and not in the ScaLAPACK routines? I cannot find anything, and also colleagues from the computing center don't see anything and do not sure that everything is ok with the ScaLAPACK routines. I've tested the program on two different platforms, IBM power5 and a Linux cluster. It has to be compiled with double precision flags.

Thanks!
Klystron



program eigen

implicit none

! BLACS stuff
integer iam,nprocs,context,pg_rows,pg_cols,myrow,mycol,info,pgrow,pgcol,opt
integer BLOCKROWS,BLOCKCOLS,SIZE

! array descriptor
integer, dimension(9) :: desc_matrix

! processor's local memory
real,allocatable :: loc_matrix(:,:),loc_matrix2(:,:),work(:)
real,allocatable :: eigre(:),eigre2(:),eigim(:),eigim2(:),loc_z(:,:)
complex,allocatable :: loc_cmatrix(:,:),loc_cz(:,:),cwork(:)
complex,allocatable :: eigval(:),eigval2(:),eigval3(:)
integer i,j,loc_rows,loc_cols,lwork
integer,allocatable :: sort(:),sort2(:),sort3(:)
real diff,maxdiff

! ScaLAPACK functions
integer iceil

! This program computes eigenvalues of an invented Hessenberg matrix by
! means of three different methods: DGEEV, PDLAHQR, and PZLAHQR. The results
! are compared. The variable opt influences the program parameters
! BLOCKROWS, BLOCKCOLS, SIZE, and the matrix elements.

opt=1

select case(opt)
case(1) ! all the three routines agree
BLOCKROWS=64
BLOCKCOLS=64
SIZE=64
case(2) ! only DGEEV und PZLAHQR agree?
BLOCKROWS=64
BLOCKCOLS=64
SIZE=65
case(3) ! all the three routines agree again
BLOCKROWS=65
BLOCKCOLS=65
SIZE=65
case(4) ! all the three routines still agree (diagonal matrix)
BLOCKROWS=64
BLOCKCOLS=64
SIZE=65
case(5) ! all routines disagree
BLOCKROWS=64
BLOCKCOLS=64
SIZE=64
case default
BLOCKROWS=64
BLOCKCOLS=64
SIZE=64
end select

! initialize process grid
call BLACS_PINFO(iam,nprocs)
if(nprocs .ne. 1) then
write(*,*) 'only one processor please...'
stop
endif
pg_rows=1 ! process grid rows
pg_cols=1 ! process grid columns
call BLACS_GET(0,0,context)
call BLACS_GRIDINIT(context,'R',pg_rows,pg_cols)
call BLACS_GRIDINFO(context,pg_rows,pg_cols,myrow,mycol)

! generate matrix descriptor
loc_cols=iceil(iceil(SIZE,BLOCKCOLS),pg_cols)*BLOCKCOLS
loc_rows=iceil(iceil(SIZE,BLOCKROWS),pg_rows)*BLOCKROWS
write(*,*) 'loc_rows =',loc_rows,'loc_cols =',loc_cols
call DESCINIT(desc_matrix,SIZE,SIZE,BLOCKROWS,BLOCKCOLS,0,0, &
context,loc_rows,info)

allocate(loc_matrix(loc_rows,loc_cols),loc_matrix2(loc_rows,loc_cols), &
loc_cmatrix(loc_rows,loc_cols),loc_z(loc_rows,loc_cols), &
loc_cz(loc_rows,loc_cols),eigval(SIZE),eigval2(SIZE),eigval3(SIZE), &
eigre(SIZE),eigre2(SIZE),eigim(SIZE),eigim2(SIZE),sort(SIZE), &
sort2(SIZE),sort3(SIZE))

! generate an upper Hessenberg matrix
loc_matrix=0.0
do j=1,SIZE
do i=1,j+1
if(opt .eq. 5) then
loc_matrix(i,j)=2.0*real(i)+real(j)
else
if(opt .ne. 4) loc_matrix(i,j)=0.6*sin(2.0*real(i)+real(j))
endif
if(i .eq. j) loc_matrix(i,i)=real(i)+50.5
enddo
enddo

! copy the matrix
do j=1,SIZE
do i=1,SIZE
loc_matrix2(i,j)=loc_matrix(i,j)
loc_cmatrix(i,j)=cmplx(loc_matrix(i,j),0.0)
enddo
enddo

!
! compute eigenvalues (dgeev), result in eigval
!
lwork=10*SIZE
allocate(work(lwork))
call DGEEV('N','N',SIZE,loc_matrix,loc_rows,eigre,eigim,0,1,0,1,work,lwork, &
info)
do i=1,SIZE
eigval(i)=cmplx(eigre(i),eigim(i))
enddo

! sort eigval
call evsort(SIZE,eigval,sort)

! unit matrix in loc_z, and loc_cz
loc_z=0.0
loc_cz=(0.0,0.0)
do i=1,SIZE
loc_z(i,i)=1.0
loc_cz(i,i)=(1.0,0.0)
enddo

!
! compute eigenvalues (pdlahqr), result in eigval2
!
!!$ ! workspace query
!!$ lwork=-1
!!$ call PDLAHQR(.true.,.true.,SIZE,1,SIZE,loc_matrix2,desc_matrix,eigre2, &
!!$ eigim2,1,size,loc_z,desc_matrix,work,lwork,0,0,info)
!!$ lwork=int(work(1))
lwork=1000000
write(*,*) 'PDLAHQR: lwork =',lwork
deallocate(work)
allocate(work(lwork))
! the actual PDLAHQR computation
call PDLAHQR(.true.,.true.,SIZE,1,SIZE,loc_matrix2,desc_matrix,eigre2, &
eigim2,1,size,loc_z,desc_matrix,work,lwork,0,0,info)
do i=1,SIZE
eigval2(i)=cmplx(eigre2(i),eigim2(i))
enddo

! sort eigval2
call evsort(SIZE,eigval2,sort2)

!
! compute eigenvalues (pzlahqr), result in eigval3
!

! workspace query
lwork=-1
allocate(cwork(1))
call PZLAHQR('true','true',SIZE,1,SIZE,loc_cmatrix,desc_matrix,eigval3,1, &
size,loc_cz,desc_matrix,cwork,lwork,0,0,info)
lwork=int(cwork(1))
deallocate(cwork)
allocate(cwork(lwork))
! The actual PZLAHQR computation. Attention: the first two parameters
! should be always set to true, otherwise the results cannot be trusted!!!
call PZLAHQR('true','true',SIZE,1,SIZE,loc_cmatrix,desc_matrix,eigval3,1, &
size,loc_cz,desc_matrix,cwork,lwork,0,0,info)

! sort eigval3
call evsort(SIZE,eigval3,sort3)

! check differences between eigval and eigval2, print eigval and eigval2
maxdiff=0.0
do i=1,size
diff=abs(abs(real(eigval(sort(i)))) &
-abs(real(eigval2(sort2(i))))) &
+abs(abs(aimag(eigval(sort(i)))) &
-abs(aimag(eigval2(sort2(i)))))
if(diff .gt. maxdiff) maxdiff=diff
write(*,*) i,eigval(sort(i)),eigval2(sort2(i)),diff
enddo
write(*,*) 'DGEEV and PDLAHQR: maxdiff =',maxdiff

! check differences between eigval and eigval3
maxdiff=0.0
do i=1,size
diff=abs(abs(real(eigval(sort(i)))) &
-abs(real(eigval3(sort3(i))))) &
+abs(abs(aimag(eigval(sort(i)))) &
-abs(aimag(eigval3(sort3(i)))))
if(diff .gt. maxdiff) maxdiff=diff
enddo
write(*,*) 'DGEEV and PZLAHQR: maxdiff =',maxdiff

! check differences between eigval2 and eigval3
maxdiff=0.0
do i=1,size
diff=abs(abs(real(eigval2(sort2(i)))) &
-abs(real(eigval3(sort3(i))))) &
+abs(abs(aimag(eigval2(sort2(i)))) &
-abs(aimag(eigval3(sort3(i)))))
if(diff .gt. maxdiff) maxdiff=diff
enddo
write(*,*) 'PDLAHQR and PZLAHQR: maxdiff =',maxdiff

call BLACS_GRIDEXIT(context)
call BLACS_EXIT(0)

end program eigen



subroutine evsort(num,val,sort)

! sort values (rather primitive)
! on return, array sort contains the indices for the entries of array val
! such that they are sorted in descending order of the real parts

implicit none

integer num
complex val(num)
integer sort(num)

integer sortlen,i,j,k

sortlen=1
sort(1)=1
do i=2,num
j=sortlen
do while(real(val(i)) .gt. real(val(sort(j))) .and. (j .gt. 0))
j=j-1
if(j .eq. 0) exit
enddo
do k=sortlen,j+1,-1
sort(k+1)=sort(k)
enddo
sort(j+1)=i
sortlen=sortlen+1
enddo
return

end subroutine evsort

[Klystron]

Please, forget case "opt=5". I've verified that, in that case, the eigenvalues are just extremely ill-conditioned. So I think PZLAHQR is ok. However, I'm still concerned with PDLAHQR.