Closed loiseaujc closed 1 month ago
Seems like the issue is here: https://github.com/fortran-lang/stdlib/blob/42182b080d232fb0704ae1578a9235fe3cf57ef7/src/stdlib_linalg_least_squares.fypp#L138
x
should not be allocated using mold=b
but by extraction of size(A,dim=2), and if rank(b)>1 then also the size of b in the second dimension.
@perazz what are your thoughts?
Sorry I hadn't thought about this. *GELSD
require that x
is sized at least as b
, because they take the rhs [n]
on input and return x[m]
on output.
I totally agree that this is misleading. So, there is no way to avoid an additional allocation, in order to return an appropriately sized x
. The only thing we can do to save an allocation is maybe in the subroutine
interface, to let the user provide a larger x
, and only return the appropriate part in that case.
Awesome. I'll try it out as soon as it's merged. Closing the issue. Thx.
Hej,
I'm really happy that the
linalg
module is on its way to being production-ready. I have already played a bit with the latest additions inv0.6.0
and observed something an undesirable (in my opinion) behaviour regarding the output oflstsq
. Below is a MWE.When running this code, the vector
x_lstsq
computed bylstsq
is of sizem
rather than being of sizen
. I know that, by default,*gelsd
takesA
andb
as input and overwrites the firstn
entries ofb
with the solutionx
. I believe this is quite misleading and can potentially cause issues if the user is directly usingsolve_lstsq
. Callingsolve_lstsq
withx
being ann
-vector instead of anm
-vector actually raises an error since the following test then failsin e.g.
stdlib_linalg_s_solve_lstsq_one
. Are there any reasons why we should keep the output vector as anm
-vector rather than ann
-vector ? From a mathematical and practical point of view, I would expect the following behaviour:lstsq
takes anm
xn
matrixA
and anm
-vectorb
as entry and returns ann
-vector.solve_lstsq
takes as argumentsm
xn
matrixA
, anm
-vectorb
and ann
-vectorx
.