使用 LAPACK 的 Fortran2003 中的动态内存分配错误

Posted 2023-02-21

【中文标题】使用 LAPACK 的 Fortran2003 中的动态内存分配错误

【英文标题】：Dynamic memory allocation error in Fortran2003 using LAPACK

【发布时间】：2012-03-27 03:24:32

【问题描述】：

我正在为 LAPACK 的 dgetrf 和 dgetri 例程而苦苦挣扎。下面是我创建的一个子例程（变量 fit_coeffs 是在外部定义的并且是可分配的，这不是问题）。当我运行时，由于 matmul(ATA,AT) 行，当我分配 fit_coeffs 时出现内存分配错误。我通过插入一堆打印语句知道这一点。此外，调用 LAPACK 子例程后的两个错误检查语句都会打印出来，提示错误。 有谁明白这是从哪里来的？我正在使用命令编译：

gfortran -Wall -cpp -std=f2003 -ffree-form -L/home/binningtont/lapack-3.4.0/ read_grib.f -llapack -lrefblas。

提前致谢！

subroutine polynomial_fit(x_array, y_array, D)
    integer, intent(in) :: D
    real, intent(in), dimension(:) :: x_array, y_array
    real, allocatable, dimension(:,:) :: A, AT, ATA
    real, allocatable, dimension(:) :: work
    integer, dimension(:), allocatable :: pivot
    integer :: l, m, n, lda, lwork, ok

    l = D + 1
    lda = l
    lwork = l

    allocate(fit_coeffs(l))
    allocate(pivot(l))
    allocate(work(l))
    allocate(A(size(x_array),l))
    allocate(AT(l,size(x_array)))
    allocate(ATA(l,l))

    do m = 1,size(x_array),1
      do n = 1,l,1
        A(m,n) = x_array(m)**(n-1)
      end do
    end do

    AT = transpose(A)
    ATA = matmul(AT,A)

    call dgetrf(l, l, ATA, lda, pivot, ok)
    ! ATA is now represented as PLU (permutation, lower, upper)
    if (ok /= 0) then
      write(6,*) "HERE"
    end if
    call dgetri(l, ATA, lda, pivot, work, lwork, ok)
    ! ATA now contains the inverse of the matrix ATA
    if (ok /= 0) then
      write(6,*) "HERE"
    end if

    fit_coeffs = matmul(matmul(ATA,AT),y_array)

    deallocate(pivot)
    deallocate(fit_coeffs)
    deallocate(work)
    deallocate(A)
    deallocate(AT)
    deallocate(ATA)
  end subroutine polynomial_fit

【问题讨论】：

具体报错信息是什么？数组够大吗？

嗨，谢谢。错误是 *** glibc detected *** ./a.out: malloc(): memory corruption: 0x00000000020dfbd0 ***，后面是 Backtrace 和 Memory map，我真的不明白。根据我从 LAPACK 库中了解的内容，我已经检查了数组是否足够大。我很确定问题出在外部函数调用上，因为我插入了错误检查语句。

您是否尝试过使用编译器开关-fcheck=all？我还会尝试先将matmul(ATA,AT) 分配给一个临时变量，然后在第二次调用matmul 时使用该变量。

在这种特定情况下它可能会或可能不会有帮助，但-fcheck=all 会打开几个 runtime 检查。

我也会尝试valgrind。一定要编译调试。

【参考方案1】：

1) fit_coeffs 在哪里声明？我看不出上面的内容是如何编译的 1b) Implicit None 是你的朋友！

2) 你在调用点的范围内确实有一个接口，不是吗？

3) dgertf 和 dgetri 想要“双精度”，而你只有一个。所以你需要 sgetrf 和 sgetri

“修复”所有这些并完成我得到的程序

Program testit

  Implicit None

  Real, Dimension( 1:100 ) :: x, y

  Integer :: D

  Interface 
     subroutine polynomial_fit(x_array, y_array, D)
       Implicit None ! Always use this!!
       integer, intent(in) :: D
       real, intent(in), dimension(:) :: x_array, y_array
     End subroutine polynomial_fit
  End Interface

  Call Random_number( x )
  Call Random_number( y )

  D = 6

  Call polynomial_fit( x, y, D )

End Program testit

subroutine polynomial_fit(x_array, y_array, D)

  Implicit None ! Always use this!!

    integer, intent(in) :: D
    real, intent(in), dimension(:) :: x_array, y_array
    real, allocatable, dimension(:,:) :: A, AT, ATA
    real, allocatable, dimension(:) :: work, fit_coeffs
    integer, dimension(:), allocatable :: pivot
    integer :: l, m, n, lda, lwork, ok

    l = D + 1
    lda = l
    lwork = l

    allocate(fit_coeffs(l))
    allocate(pivot(l))
    allocate(work(l))
    allocate(A(size(x_array),l))
    allocate(AT(l,size(x_array)))
    allocate(ATA(l,l))

    do m = 1,size(x_array),1
      do n = 1,l,1
        A(m,n) = x_array(m)**(n-1)
      end do
    end do

    AT = transpose(A)
    ATA = matmul(AT,A)

    call sgetrf(l, l, ATA, lda, pivot, ok)
    ! ATA is now represented as PLU (permutation, lower, upper)
    if (ok /= 0) then
      write(6,*) "HERE"
    end if
    call sgetri(l, ATA, lda, pivot, work, lwork, ok)
    ! ATA now contains the inverse of the matrix ATA
    if (ok /= 0) then
      write(6,*) "HERE"
    end if

    fit_coeffs = matmul(matmul(ATA,AT),y_array)

    deallocate(pivot)
    deallocate(fit_coeffs)
    deallocate(work)
    deallocate(A)
    deallocate(AT)
    deallocate(ATA)
  end subroutine polynomial_fit

这将运行到完成。如果我省略界面，我会打印两次“HERE”。如果我使用 d 版本，我会遇到段错误。

这能回答你的问题吗？

【讨论】：

您好，谢谢您的回答。 fit_coeffs 是在模块开头声明的可分配数组（包含公共和隐式 none 语句），并且该模块在调用 polynomial_fit 的子例程中“使用”。所以接口块会干扰使用模块命令。当我使用 dgetrf 时，我得到一个 malloc(): 内存损坏错误，但是当我使用 sgetrf 时，我只是得到一个段错误。

您是否有一个简短而完整的程序来证明您的问题？

好消息，无需进一步分析。我认为您对 sgetrf 与 dgetrf 的看法是正确的。我遇到的段错误似乎来自其他原因，但现在它已修复，sgetrf 和 sgetri 工作正常，没有内存分配错误！谢谢。

@Ian：很好看！我认为这很好地说明了（模块辅助）接口检查是如何在 fortran 90 中取得巨大进步的。