Dynamic memory allocation error in Fortran2003 using LAPACK

Question

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.

Thanks in advance!

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

What is the specific error message? Are the arrays big enough? — M. S. B.
– M. S. B., Commented Mar 9, 2012 at 17:22
Hi, thanks. The error is *** glibc detected *** ./a.out: malloc(): memory corruption: 0x00000000020dfbd0 ***, followed by a Backtrace and Memory map, which I don't understand really. I''ve checked that the arrays are large enough, according to what I understand from the LAPACK libraries. I'm fairly sure the problem lies in the external function calls, because of the error-checking statements that I inserted. — Taylor
– Taylor, Commented Mar 9, 2012 at 17:26
Have you tried using the compiler switch -fcheck=all? I would also try assigning matmul(ATA,AT) to a temporary variable first and then use that variable in the second call to matmul. — eriktous
– eriktous, Commented Mar 9, 2012 at 17:55
It may or may not help in this specific case, but -fcheck=all turns on several runtime checks. — eriktous
– eriktous, Commented Mar 9, 2012 at 20:34
I would also try valgrind. Be sure to compile for debugging. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Mar 9, 2012 at 20:40

Ian Bush · Accepted Answer · 2012-03-12 09:29:25Z

4

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

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

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

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

Does this answer your question?

answered Mar 12, 2012 at 9:29

Ian Bush

411 bronze badge

Sign up to request clarification or add additional context in comments.

4 Comments

Taylor Over a year ago

Hi, thank you for your answer. fit_coeffs is an allocatable array declared at the beginning of a module (containing public and implicit none statements), and that module is 'used' within a subroutine from which polynomial_fit is called. So an interface block interferes with the use module command. When I use dgetrf, I get a malloc(): memory corruption error, but when I use sgetrf, I simply get a segfault.

Ian Bush Over a year ago

Do you have a short, complete program that demonstrated your problem then?

Taylor Over a year ago

Good news, no further analysis is required. I think you were correct about the sgetrf vs. dgetrf. The segfault I encountered appears to be from something else, but now that it's fixed, sgetrf and sgetri work just fine, no memory allocation error! Thank you.

eriktous Over a year ago

@Ian: Well spotted! I think this is an excellent illustration of how (module assisted) interface checking was a great advance in fortran 90.

Collectives™ on Stack Overflow

Dynamic memory allocation error in Fortran2003 using LAPACK

1 Answer 1

4 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

4 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related