您所在位置：主题：请教一个关于LAPACK的问题（环境IVF+MVS+MKL）
kabalaka
[专家分：0] 发布于 2012-02-17 02:33:00
大家好，请教一个问题。
我想用LAPACK里面的矩阵奇点分解子程序ZGESVD分解一个矩阵（132*132），但是直接崩溃（不是计算不出结果，而是直接退出程序"segmentation fault"），而如果分解稍小一点的124*124的，那么还可以进行，但是也比较慢（内存占用也非常大）。我认为132*132的矩阵并不算大，真正要分解的矩阵应该是上千的大小。但是132*132的ZGESVD矩阵分解却不是慢而是程序崩溃却很奇怪。
我又实验了类似的子程序ZGGSVD，分解同样的矩阵，却工作正常。甚至更大的也计算正常。让我很不理解，因此上来请教各位，哪里出了问题，我的程序附在后面了。
请各位不吝赐教，先谢谢了！
程序如下运行开始的时候会需要输入wcs的值，矩阵的size=(wcs*2+1)*4，程序会自动生成矩阵并做奇点分解，我调用了ZGGSVD和ZGESVD，wcs<15时一切正常，wcs>=16时ZGGSVD正常，ZGESVD就崩溃了。
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
我想用LAPACK里面的矩阵奇点分解子程序ZGESVD分解一个矩阵（132*132），但是直接崩溃（不是计算不出结果，而是直接退出程序"segmentation fault"），而如果分解稍小一点的124*124的，那么还可以进行，但是也比较慢（内存占用也非常大）。我认为132*132的矩阵并不算大，真正要分解的矩阵应该是上千的大小。但是132*132的ZGESVD矩阵分解却不是慢而是程序崩溃却很奇怪。
我又实验了类似的子程序ZGGSVD，分解同样的矩阵，却工作正常。甚至更大的也计算正常。让我很不理解，因此上来请教各位，哪里出了问题，我的程序附在后面了。
请各位不吝赐教，先谢谢了！
程序如下运行开始的时候会需要输入wcs的值，矩阵的size=(wcs*2+1)*4，程序会自动生成矩阵并做奇点分解，我调用了ZGGSVD和ZGESVD，wcs<15时一切正常，wcs>=16时ZGGSVD正常，ZGESVD就崩溃了。
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