Ошибка динамического выделения памяти в Fortran2003 с использованием LAPACK - PullRequest
Новая
       22

Ошибка динамического выделения памяти в Fortran2003 с использованием LAPACK

1 голос
/ 09 марта 2012

Я борюсь с подпрограммами LAPACK dgetrf и dgetri.Ниже приведена подпрограмма, которую я создал (переменная fit_coeffs определяется внешне и может быть размещена, это не проблема).При запуске я получаю ошибки выделения памяти, которые появляются при назначении fit_coeffs из-за строки matmul (ATA, AT).Я знаю это, вставив кучу утверждений для печати.Кроме того, оба оператора проверки ошибок после вызовов подпрограмм 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

1 Ответ

4 голосов
/ 12 марта 2012

1) Где объявлен fit_coeffs? Я не могу понять, как вышесказанное может даже скомпилировать 1b) Неявный никто не твой друг!

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

Это работает до завершения. Если я опущу интерфейс, я получу «ЗДЕСЬ» напечатан дважды. Если я использую версии d, я получаю ошибки seg.

Это отвечает на ваш вопрос?

...