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)
call dgels('N', size(x_array), D+1, 1, A, size(x_array), y_array, size(y_array), work, lwork, ok)
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. )