LAPACKを使用したFortran2003での動的メモリ割り当てエラー

Question

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 

Answer 1

1）fit_coeffsはどこで宣言されていますか？私は上記のコンパイル方法を見ることができません 1b）暗黙のうちにあなたの友人はいません！

2）コールポイントのスコープ内にインターフェイスがありますか？

3）dgertfとdgetriは、単精度の間に倍精度を求めます。だから、私は

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 

を取得sgetrfとsgetri

これらすべてを「修正」して、プログラムをcompleteingを必要とするこれが完了するまで実行されます。インターフェイスを省略すると、「ここに」が2回印刷されます。私がdバージョンを使用すると、segフォルトが発生します。

これはあなたの質問にお答えしますか？