使用LAPACK的Fortran2003中的动态内存分配错误 [英] Dynamic memory allocation error in Fortran2003 using LAPACK

问题描述

我在LAPACK的 dgetrf 和 dgetri 例程中苦苦挣扎。下面是我创建的子例程（变量fit_coeffs是在外部定义的并且可以分配，这不是问题）。运行时，由于matmul（ATA，AT）行，我得到内存分配错误，这些错误在我分配fit_coeffs时出现。我从插入一堆打印语句知道这一点。另外，在对LAPACK子例程的调用之后，两个错误检查语句都被打印出来，提示错误。
有人知道这是哪里来的吗？我正在使用以下命令进行编译：

I'm struggling with LAPACK's dgetrf and dgetri routines. Below is a subroutine I've created (the variable fit_coeffs is defined externally and is allocatable, it's not the problem). When I run I get memory allocation errors, that appear when I assign fit_coeffs, due to the matmul(ATA,AT) line. I know this from inserting a bunch of print statements. Also, both error checking statements after calls to LAPACK subroutines are printed, suggesting an error. Does anyone understand where this comes from? I'm compiling using the command:

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

推荐答案

1）在哪里声明fit_coeffs？我看不出上面的代码甚至不能编译
1b）隐式没有一个是您的朋友！

1) Where is fit_coeffs declared? I can't see how the above can even compile 1b) Implicit None is your friend!

2）您在调用时确实在作用域内有一个接口点，不是吗？

2) You do have an interface in scope at the calling point, don't you?

3）dgertf和dgetri在单身时希望双精度。因此，您需要sgetrf和sgetri

3) dgertf and dgetri want "double precision" while you have single. So you need sgetrf and sgetri

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

"Fixing" all these and completeing the program I get

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版本，则会出现段错误。

This runs to completion. If I omit the interface I get "HERE" printed twice. If I use the d versions I get seg faults.

这能回答您的问题吗？

这篇关于使用LAPACK的Fortran2003中的动态内存分配错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！