Strange phenomenon of ZGESVD

[chenxikane]

Hello, everyone, thank you for reading my topic.

I have encountered a strange phenomenon, when I call ZGESVD to Singular Value Decomposition, I discover it can not handle the matrix which size bigger than 132 in my program(it needs unbelievable memory >2.7G), however the similar ZGGSVD works well for the same size of matrix. So I am confused to this. Then I tried ZGGSVD and ZGESVD both, if matrix is not bigger than 132, both of ZGGSVD and ZGESVD can be functional both, however if W>=16, ZGGSVD works well, but ZGESVD crushed. I do not understand this. And I checked my program many times, I can not find what could cause this strange result. Anyway, I think even I inputted a weird matrix, the result of ZGESVD should not be crush of program.

Therefore, if anyone can help me, I will very appreciate.

My program is following, my FORTRAN environment is IVF2011+MVS+MKL(lapack). And one can input the wcs, then the later subroutine will generate a square matrix for SVD. If wcs <=15, everything is OK, if wcs>=16, ZGGSVD is OK, ZZGESVD will crush.
!=================================================
program tryZGESVD
implicit none
external zgesvd
external zggsvd

integer :: Wcs,Sz
complex(8),dimension(4) :: MtrxElmt

!input and output for ZGGSVD and ZGGSVD
integer :: lwork,info,SDK,SDL
real(8),dimension(:),allocatable :: sglV(:),rwork2(:),rwork1(:),alpha,beta
complex(8),dimension(:,:),allocatable :: SVDU,SVDVt,SVDQ,Mtrx00,Mtrx,Mtrxb
complex,allocatable,dimension(:) :: work0,work1,work2,iwork

write(*,*) "input Wcs(INTEGER), the size of matirx = (wcs*2+1)*4"
write(*,*) "if wcs<=15, everything works well, but if wcs>=16, ZGESVD directly crushed"
read(*,*) wcs

sz=(wcs*2+1)*4

!Establish Matrix, by the following subroutine, just a sz by sz matrix, which only lower-left block is valued
allocate(Mtrx00(sz,sz))
allocate(Mtrx(sz,sz))
Mtrx00=0
MtrxElmt=(/-1.0d0,-0.1d0,-0.1d0,-1.0d0/)
call MFSup1(wcs*2+1,wcs*2+1,2,2,(/Wcs,wcs+1/),(/wcs+1,Wcs/),MtrxElmt,Mtrx00(6*Wcs+4:Sz,1:wcs*2+1))


allocate(SVDU(sz,sz))
allocate(SVDvt(sz,sz))
allocate(SVDQ(sz,sz))

!zggsvd can give the SVD when wcs>16, i.e. sz>132
Mtrx=Mtrx00
allocate(alpha(sz))
allocate(beta(sz))
allocate(work1(4*sz))
allocate(Mtrxb(sz,sz))
allocate(iwork(sz))
allocate(rwork1(sz*2))
Mtrxb=0
call zggsvd( "N", "N", "Q", sz,sz,sz,SDk,SDl,Mtrx,sz,Mtrxb,sz, alpha, beta,SVDu,sz,SVDvt,sz,SVDQ,sz, work1, rwork1, iwork, info)
!write(*,*) "SVDQ",SVDQ(1:sz,1)

!zgesvd can NOT give the SVD when wcs>16, i.e. sz>132
allocate(sglV(sz))
allocate(rwork2(sz*5))
SVDU=0
SVDvt=0
Mtrx=Mtrx00
!call zgesvd( jobu,jobvt, m, n, a, lda, s, SVDU, ldu, SVDvt,ldvt, work, lwork, rwork, info)
call zgesvd( "A", "A",Sz,Sz,Mtrx, Sz,sglV, SVDU, Sz, SVDVt, Sz, rwork2, -1, rwork2, info)
lwork=rwork2(1)
allocate(work0(lwork))
write(*,*) lwork
rwork2=0
Mtrx=Mtrx00
!call zgesvd( jobu,jobvt, m, n, a, lda, s, SVDU, ldu, SVDvt,ldvt, work, lwork, rwork, info)
call zgesvd( "A", "A",Sz,Sz,Mtrx, Sz,sglV, SVDU, Sz, SVDVt, Sz, work0, lwork, rwork2, info)
deallocate(work0)
!write(*,*) "SVDvt",SVDvt(1:sz,1)
write(*,*) "sglV",sglV
end program

! ============================================================
subroutine MFSup1(SizeLTotal,SizeCTotal,BlockL,BlockC,SizeL,sizeC,BlockValue,res)
!This give all kinds of matrix block.
! ============================================================
implicit none
integer,intent(in) :: SizeLTotal,SizeCTotal,BlockL,BlockC
integer,intent(in) :: SizeL(BlockL),SizeC(BlockC)
complex(8),intent(in),dimension(BlockL*BlockC) :: BlockValue
complex(8),intent(out),dimension(SizeLTotal,SizeCTotal) :: res
integer :: StartL(BlockL),EndL(BlockL),StartC(BlockC),EndC(BlockC)
integer :: i1,i2,InitialL,FinalL,InitialC,FinalC

if (Sum(SizeL).ne.SizeLTotal.or.Sum(SizeC).ne.SizeCTotal) then
continue
stop 'MFSup1'
end if

InitialL=1
FinalL=0
StartL(1)=InitialL
EndL(BlockL)=SizeLTotal
do i1=1,BlockL-1
InitialL=InitialL+SizeL(i1)
StartL(i1+1)=InitialL
FinalL=FinalL+SizeL(i1)
EndL(i1)=FinalL
end do
InitialC=1
FinalC=0
StartC(1)=InitialC
EndC(BlockC)=SizeCTotal
do i1=1,BlockC-1
InitialC=InitialC+SizeC(i1)
StartC(i1+1)=InitialC
FinalC=FinalC+SizeC(i1)
EndC(i1)=FinalC
end do

do i1=1,BlockL
do i2=1,BlockC
call MFOri2(SizeL(i1),SizeC(i2),BlockValue((i1-1)*BlockC+i2),res(StartL(i1):EndL(i1),StartC(i2):EndC(i2)))
end do
end do

end subroutine MFSup1

! ============================================================
subroutine MFOri2(Size1,Size2,DiagE,res)
!This give the non-square matrix with only semi-diagonal elements.
! ============================================================
implicit none
integer,intent(in) :: size1,Size2
complex(8),intent(in) :: DiagE
complex(8),intent(out),dimension(size1,size2) :: res
integer :: i1,i2,temp
res=0

if (size1.gt.size2) then
temp=size1-size2
do i1=1,size2
res(i1,i1)=DiagE
res(i1+temp,i1)=DiagE
end do
else if (size1.lt.size2) then
temp=size2-size1
do i1=1,size1
res(i1,i1)=DiagE
res(i1,i1+temp)=DiagE
end do
else if (size1.eq.size2) then
do i1=1,size1
res(i1,i1)=DiagE
end do
end if

end subroutine MFOri2