LAPACK：大型矩阵的 ZHEEV 例程失败

Question

我有一个 Fortran 90 代码，它在块对角化之后找到自旋链哈密顿量的特征值。我在生成每个块时将它们对角化，这似乎工作得很好，直到块的大小变得太大，此时我得到错误。

在 spinchain.exe 中的 0x05E8A247 (mkl_avx2.dll) 抛出异常：0xC0000005：访问冲突写入位置 0x054E8000。

我在 Visual Studio 2017 年使用带有 MKL 库的英特尔 Fortran 编译器。据我所知，当块矩阵的阶数大于 60 左右时，我会收到此错误。运行在 Release 模式下（与 Debug 相对），代码将运行一直运行，没有任何抱怨，但我认为输出不正确。

我写了一个简短的代码来创建一个 N x N 厄密矩阵，并通过 zheev 例程运行它。为简单起见，我选择了全矩阵，它有 N-1 个零特征值和一个单一的特征值 N。我简单地遍历 N 的值，生成矩阵并用 zheev() 对角化它，打印 N，最大特征值，以及所有特征值的总和。我发现在 N=51 之前一切都很好，此时最大特征值 returns 52 然后当我尝试解除分配错误的数组时出现异常 Critical error detected c0000374.

命中"continue"异常变成

Eigtest.exe 中 0x77A5A879 (ntdll.dll) 的未处理异常：0xC0000374：堆已损坏（参数：0x77A95910）。

代码如下：

program Eigenvalues
    implicit none
    integer, parameter :: dp = kind(0.d0)
    integer :: N, lwork, info, i
    real(dp), allocatable :: lambda(:), work(:), rwork(:)
    complex(dp), allocatable :: H(:,:)

    do N = 1,70

        lwork = 3*N
        allocate(H(N,N))
        allocate(lambda(N))
        allocate(work(lwork))
        allocate(rwork(lwork))

        H = dcmplx(1.0_dp,0)

        call zheev('N','U',N,H,N,lambda,work,lwork,rwork,info)

        if (info == 0) then
                write(*,'(I5, F16.12, F16.12)'), N, lambda(N), sum(lambda)
        else
            print*, "diagonalization failed: info = ", info
            read*, i
            stop
        end if

        deallocate(H,lambda,work,rwork)

    end do

    read*, i

end program

我还写了另一个版本的代码，我在开头将 N 固定为参数并定义了没有 allocatable 的数组，这将运行无一例外地给出N=51 的错误特征值和 N=51 以上的大多数值。

如果我的代码不好请见谅，我是物理学家，不是程序员。意见或建议表示赞赏。

Answer 1

根据 zheev 的手册页，数组 work 应该是 complex*16，所以我们可能需要

real(dp), allocatable :: lambda(:), rwork(:)
complex(dp), allocatable :: H(:,:), work(:)

（这似乎在我的电脑上使用 gfortran）。

LAPACK：大型矩阵的 ZHEEV 例程失败

LAPACK: ZHEEV routine failing for large matrices

fortran

lapack

fortran90

intel-mkl

intel-fortran