运行dgesvd时内存引用无效

Question

我正在尝试使用dgesvd计算SVD，但遇到以下错误：

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x10af44a69
#1  0x10af43e35
#2  0x7fff8ae4b529 Segmentation fault: 11

我的代码如下

program calc_pseud3
implicit none
character:: jobu*1, jobvt*1
integer:: N, M, i, j, lwork, info, LDg, LDU, LDVT
real*8:: t0,t1,t2,t3
real*8, allocatable, dimension(:) :: g_, work, S
real*8, allocatable, dimension(:,:) :: g,U,VT,S2,I_check,g_i

call cpu_time(t0)


N = 512*512
M = 2560
lwork = 280000
jobu = 'A'
jobvt = 'A'
LDg = N
LDU = N
LDVT = M

allocate(g_(N*M))
allocate(g(LDg,M))
allocate(g_i(M,N))
allocate(work(lwork))
allocate(S(M))
allocate(S2(N,M))
allocate(U(LDU,N))
allocate(VT(LDVT,M))
allocate(I_check(M,M))

open(1, file='projection.txt')

do i = 1,M
  print *, i
  do j = 1,N
    read(1,*) g_(N*(i-1)+j)
  end do
end do

call cpu_time(t1)
print *, t1-t0
g = reshape(g_, [N, M])
print *, 'Shape of matrix G is ', shape(g)

call cpu_time(t2)
print *, t2-t1

call dgesvd(jobu,jobvt,N,M,g,LDg,S,U,LDU,VT,LDVT,work,lwork,info)

print *, 'here'
S2 = 0.0*g
do j = 1,N
  do i = 1,M
    if (j == i .and. S(i) >= 10**(-4)*S(1)) then        ! If M > N, replace S(i) with S(j)
      S2(j,i) = 1/S(i)
    end if
  end do
end do

g_i = matmul( transpose(VT),matmul(transpose(S2),transpose(U)) )
I_check = matmul(g_i,g)
write (*,1) sum(sum(I_check,1))
1 format(1f10.5)
print *, 'info = ',info

open(4, file='Pseud_inverse.txt')
do i = 1,N*M
write(4,3) g_i
end do
close(4)
3 format(1f10.3)


call cpu_time(t3)
print *, t3-t2


deallocate(g_)
deallocate(g)
deallocate(g_i)
deallocate(work)
deallocate(S)
deallocate(S2)
deallocate(U)
deallocate(VT)
deallocate(I_check)


end

调用dgesvd时出现错误消息，因此我认为它一定来自此子例程。我查阅了该例程的文档，似乎没有违反任何要求。

我也用相同的输入矩阵做了与上面相同的计算，但是jobu=jobvt='N'，所以不计算U和VT，因此它们的内存分配并不重要，它起作用了。使用上面的代码仍然没有成功。

Answer 1

首先，确保您有足够的内存来分配这么大的矩阵，它大约是5 5GB，即使这不是错误的来源。无论如何，在这种情况下，我建议对矩阵执行截断奇异值分解，因为这样可以节省大量内存。

此外，文档中对lwork的定义如下：

数组乘积的维数。

LWORK (1，5* >= ( M，N))对于路径(参见代码中的注释)：- PATH 1(M远大于N，JOBU='N') - PATH 1t (N远大于M，JOBVT='N') LWORK >= MAX(1,3*MIN(M，N) + MAX(M，N),5*MIN(M，N))对于其他路径为了获得良好的性能，LWORK通常应该更大。

将N作为行传递给该子例程，将M作为列传递给子例程。请注意，在文档中，它使用N作为列，使用M作为行。也许这里有一个混淆。请仔细检查以确保。因此，对于您的情况，lwork将为269824，低于您指定的值。但是，我会尝试坚持使用文档。

最后，我前段时间写了一个Fortran模块，它有助于使用Lapack库和其他有用的计算执行SVD (截断节省内存)。你可以在here上找到它。请随时查看并使用它，如果它有帮助。注意，它适用于单精度的Lapack子例程，但对于双精度的子例程，您应该不会有任何问题。