Dynamic memory allocation and dgetrf/dgetri

[taylor]

Hi everyone,

Below is a subroutine I've created (the variable fit_coeffs is defined externally and is allocatable, it's not the problem). When I compile 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 are invoked from the dgetrf and dgetri routines. WHy might this be happening? I'm compiling using the command:

gfortran -Wall -cpp -std=f2003 -ffree-form -L/home/***USER***/lapack-3.4.0/ file.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

[Julien Langou]

No idea why your matmul fails.

Oh me, oh my, you are constructing a Vandermonde matrix for the interpolation.
This is what A is, right? So that's going to be seriously ill-conditionned. Do not take n too large.

Here are two better ways to do what you are doing. This will work for any A.
At some point you should probably take into account that you do polynomial interpolation / curve fitting.

Method 1 (close to yours, use the normal equations)
1) do not allocate AT, (see item #2,)
2) do not use matmul for ATA, construct ATA using the BLAS function DSYRK, this is half of the flops and is optimized BLAS routine, should be faster than matmul,
3) do not use DGETRF for factoring SPD matrix, you should use DPOTRF if DPOTRF fails, this means that your problem is ill-conditionned, this is
a bad sign so you can probably stop right there, or you can back up on DSYTRF (assume symmetric only) but I would try to have a well-conditionned problem.
4) do not use DGETRI+MATMUL (so multiply by the inverse for a solve), use DGETRS, and in this case use DPOTRS (or DSYTRS).

So I repeat three lines: 1) DSYRK, 2) DPOTRF, (check if INFO is zero), 3) DPOTRS.

Method 2:
this is much more stable, so for a Vandermonde method you should give it a thought, try to call DGELS right away. (So the code is one line.)
( You do not use the normal equations but a QR factorization of your original matrix. No matmul involved either )

Cheers,
Julien.

[taylor]

Hi, you guessed it! I as concerned about the matrix being singular, and now that I look closer at the value of my 'ok' parameter I find that it is equal to 6, indicating that U(6,6) from the LU factorization is exactly zero -> there will be a problem inverting U.
I will take your suggestions into account carefully and report back soon. Thank you!

[taylor]

Hi Julien,
I'm afraid the first part of your answer doesn't make any sense to me. I don't know anything about the functions that you mentioned, so I read about them to try to understand what they do, but I don't understand how to use them, or how to string them together to get the result that I want.

So I tried your second option, which I had actually tried earlier before resorting to the code that was in the original post. I don't remember why I couldn't get it working then, but I've reimplemented it like this:

subroutine polynomial_fit(x_array, y_array, D)
integer, intent(in) :: D
real, intent(inout), dimension(:) :: x_array, y_array
real, allocatable, dimension(:,:) :: A
real, allocatable, dimension(:) :: work
integer :: m, n, lwork, ok

allocate(work(D+1))
allocate(A(size(x_array),D+1))

do m = 1,size(x_array),1
do n = 1,D+1,1
A(m,n) = x_array(m)**(n-1)
end do
write(6,*) A(m,:)
end do

call dgels('N', size(x_array), D+1, 1, A, size(x_array), y_array, size(y_array), work, -1, ok)
lwork = D+1
call dgels('N', size(x_array), D+1, 1, A, size(x_array), y_array, size(y_array), work, lwork, ok)
end subroutine polynomial_fit

Anyways, what happens is when I do a query (the first time dgels is called), work(1) = 0, and when I pass that two the subroutine on the second call, I get the error:
** On entry to DGELS parameter number 10 had an illegal value
...which is the lwork argument. When I manually set lwork to something, such as D+1 as I have above, I get the same error. So apparently both 0 and D+1 are illegal? Am I doing something wrong?

[taylor]

A quick update: I looked carefully at the code for dgels.f and xerbla.f and determined that any value I had tried for lwork was too small, and would produce an INFO of -10 (I still don't know why doing a query with lwork = -1 does not yield a correct value for work(1) though). So, I manually changed the lwork value that I supplied dgels with so that it did not trip the INFO = -10 line in dgels.f. Then when I run the program, I get:

*** glibc detected *** ./a.out: free(): invalid next size (fast): 0x000000000122adb0 ***
======= Backtrace: =========
(stuff)........
======= Memory map: ========
(stuff)........
Aborted

This is basically the same error that I was getting earlier (when I was using dgetrf and dgetri), except that I was getting a malloc() error instead of a free() error.

[Julien Langou]

I'm afraid the first part of your answer doesn't make any sense to me. I don't know anything about the functions that you mentioned, so I read about them to try to understand what they do, but I don't understand how to use them, or how to string them together to get the result that I want.

This is fine. This is not a stable option anyhow. (This is why we do not offer in LAPACK a driver for it in the first place.)
Let us focus on the second option then.

So I am reading the code. You provide size(x_array) x values and size(x_array) y values and a degree D and you want to find the polynomial of degree D that best fits (in the least squares sense) the size(x_array) points (x,y). Fine.

In this case, here what you need to do:
1) allocate and fill the Vandermonde matrix A with size(x_array)-by-(D+1) points,
what you are doing looks OK,
2) LWORK needs to be of size: 2*size(x_array)
( I am assuming size(x_array) >= D+1 ... and I guess it better be )
3) Now you call DGELS just once with:

Code: Select all

lwork = 2*size(x_array)
allocate(work(lwork))
call dgels('N', size(x_array), D+1, 1, A, size(x_array), y_array, size(y_array), work, lwork, ok)
deallocate(lwork)


The D+1 first entries of y_array should contain the coefficients of the polynomial of degree D that best fits (in the least squares sense) your initial data points.

Some more feedback.

Anyways, what happens is when I do a query (the first time dgels is called), work(1) = 0,

mmm that's weird, if you call DGELS with LWORK=-1, and if OK=0, then WORK(0) should contain the optimal value for LWORK. You can use 2*size(x_array) for now. Please check OK as well.

and when I pass that two the subroutine on the second call, I get the error:
** On entry to DGELS parameter number 10 had an illegal value
...which is the lwork argument. When I manually set lwork to something, such as D+1 as I have above, I get the same error. So apparently both 0 and D+1 are illegal?

Yes yes both 0 and D+1 were illegal values for your settings. You need at least 2*size(x_array).

Also a deallocate of A and work at the end of your subroutine will not hurt.

( At the end, you might have a problem with assumed shape array in fortran. (The work(1)=0 is weird.) But let us try first to fix the easy stuff first. )

Julien.

[taylor]

Hi Julien, thank you for the reply, and sorry to be so late on getting back to you. I was away for the weekend.
Here is my whole subroutine, modified with your instructions:

subroutine polynomial_fit(x_array, y_array, D)
integer, intent(in) :: D
real, intent(inout), dimension(:) :: x_array, y_array
real, allocatable, dimension(:,:) :: A
real, allocatable, dimension(:) :: work
integer :: m, n, lwork, ok

allocate(fit_coeffs(D+1))
allocate(A(size(x_array),D+1))

do m = 1,size(x_array),1
do n = 1,D+1,1
A(m,n) = x_array(m)**(n-1)
end do
end do

lwork = 2*size(x_array)
allocate(work(lwork))
call dgels('N', size(x_array), D+1, 1, A, size(x_array), &
y_array, size(y_array), work, lwork, ok)

deallocate(A)
deallocate(work)
end subroutine polynomial_fit

It compiles, but when I run it I get the usual error:

*** glibc detected *** ./a.out: free(): invalid next size (normal): 0x0000000002501950 ***
======= Backtrace: =========
/lib/x86_64-linux-gnu/libc.so.6(+0x7a6e6)[0x7f841ed636e6]
/lib/x86_64-linux-gnu/libc.so.6(cfree+0x6c)[0x7f841ed679cc]
./a.out[0x40311e]
./a.out[0x4062c9]
./a.out[0x40ac70]
./a.out[0x40adff]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xed)[0x7f841ed0a30d]
./a.out[0x4011e9]
======= Memory map: ========
00400000-00420000 r-xp 00000000 08:05 8129094 /home/***USER***/code/a.out
00620000-00621000 r--p 00020000 08:05 8129094 /home/binningtont/code/a.out
00621000-00622000 rw-p 00021000 08:05 8129094 /home/binningtont/code/a.out
00622000-00f0c000 rw-p 00000000 00:00 0
024fa000-0253d000 rw-p 00000000 00:00 0 [heap]
7f8418000000-7f8418021000 rw-p 00000000 00:00 0
7f8418021000-7f841c000000 ---p 00000000 00:00 0
7f841eab3000-7f841eae8000 r-xp 00000000 08:01 6819326 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f841eae8000-7f841ece7000 ---p 00035000 08:01 6819326 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f841ece7000-7f841ece8000 r--p 00034000 08:01 6819326 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f841ece8000-7f841ece9000 rw-p 00035000 08:01 6819326 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f841ece9000-7f841ee80000 r-xp 00000000 08:01 7344398 /lib/x86_64-linux-gnu/libc-2.13.so
7f841ee80000-7f841f07f000 ---p 00197000 08:01 7344398 /lib/x86_64-linux-gnu/libc-2.13.so
7f841f07f000-7f841f083000 r--p 00196000 08:01 7344398 /lib/x86_64-linux-gnu/libc-2.13.so
7f841f083000-7f841f084000 rw-p 0019a000 08:01 7344398 /lib/x86_64-linux-gnu/libc-2.13.so
7f841f084000-7f841f08a000 rw-p 00000000 00:00 0
7f841f08a000-7f841f09f000 r-xp 00000000 08:01 7344514 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f841f09f000-7f841f29e000 ---p 00015000 08:01 7344514 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f841f29e000-7f841f29f000 r--p 00014000 08:01 7344514 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f841f29f000-7f841f2a0000 rw-p 00015000 08:01 7344514 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f841f2a0000-7f841f323000 r-xp 00000000 08:01 7344407 /lib/x86_64-linux-gnu/libm-2.13.so
7f841f323000-7f841f522000 ---p 00083000 08:01 7344407 /lib/x86_64-linux-gnu/libm-2.13.so
7f841f522000-7f841f523000 r--p 00082000 08:01 7344407 /lib/x86_64-linux-gnu/libm-2.13.so
7f841f523000-7f841f524000 rw-p 00083000 08:01 7344407 /lib/x86_64-linux-gnu/libm-2.13.so
7f841f524000-7f841f638000 r-xp 00000000 08:01 6819331 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f841f638000-7f841f837000 ---p 00114000 08:01 6819331 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f841f837000-7f841f838000 r--p 00113000 08:01 6819331 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f841f838000-7f841f83a000 rw-p 00114000 08:01 6819331 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f841f83a000-7f841f85b000 r-xp 00000000 08:01 7340055 /lib/x86_64-linux-gnu/ld-2.13.so
7f841fa49000-7f841fa4d000 rw-p 00000000 00:00 0
7f841fa58000-7f841fa5a000 rw-p 00000000 00:00 0
7f841fa5a000-7f841fa5b000 r--p 00020000 08:01 7340055 /lib/x86_64-linux-gnu/ld-2.13.so
7f841fa5b000-7f841fa5d000 rw-p 00021000 08:01 7340055 /lib/x86_64-linux-gnu/ld-2.13.so
7fffbe57e000-7fffbe59f000 rw-p 00000000 00:00 0 [stack]
7fffbe5de000-7fffbe5df000 r-xp 00000000 00:00 0 [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0 [vsyscall]
Aborted

I think the work(1) = 0 is very strange as well. On Friday I spent some time tracing through dgels.f to see where the control went when I passed lwork = -1 to it, and thought I understaood things well, but still don't see how I could get a work(1) = 0 value.

[taylor]

Good news! There were actually three bugs together, which meant that even when I'd fixed the actual memory problem, one of the others persisted so that I didn't think I'd fixed it. The problem is fixed by substituting sgetrf and sgetri for the routines dgetrf and dgetri. Very difficult to trace, I expected a more transparent error, something like 'expected real(8) but got real(4)' maybe.
Come to think of it, this may have been an issue for calling dgels as well. I'll investigate and report back. Anyways thank you for all your help!

[Julien Langou]

Yet, I really recommend you to give a shot to sgles as opposed to sgetrf and sgetri. Cheers, Julien.