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

Question

我正在努力处理LAPACK的dgetrf和dgetri例程。下面是我创建的一个子例程(变量fit_coeffs是在外部定义的，并且是可分配的，这不是问题)。运行时，由于matmul(ATA，AT)行的原因，在分配fit_coeffs时会出现内存分配错误。我从插入一堆print语句中知道了这一点。此外，还会打印调用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

Answer 1

1) fit_coeffs在哪里声明？我不明白上面的代码是怎么编译的:1)隐含的无人是你的朋友！

2)您在调用点的作用域中确实有一个接口，不是吗？

3) dgertf和dgetri需要“双精度”，而您使用的是single。所以你需要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个版本，我会得到seg错误。

这回答了你的问题吗？