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 in v0.6.0 and observed something an undesirable (in my opinion) behaviour regarding the output of lstsq. Below is a MWE.

program main use iso_fortran_env, only: output_unit use stdlib_kinds, only: dp use stdlib_linalg, only: lstsq, solve_lstsq implicit none integer :: m, n real(dp), allocatable :: A(:, :), b(:), x(:), x_lstsq(:) ! Create the problem. m = 10 ; n = 5 allocate(A(m, n)) ; call random_number(A) allocate(x(n)) ; call random_number(x) b = matmul(A, x) ! Solve the least-squares problem. x_lstsq = lstsq(A, b) write(output_unit, *) shape(A), shape(b), shape(x), shape(x_lstsq) write(output_unit, *) "Norm of the error :", norm2(x_lstsq(1:n)-x) end program main 

When running this code, the vector x_lstsq computed by lstsq is of size m rather than being of size n. I know that, by default, *gelsd takes A and b as input and overwrites the first n entries of b with the solution x. I believe this is quite misleading and can potentially cause issues if the user is directly using solve_lstsq. Calling solve_lstsq with x being an n-vector instead of an m-vector actually raises an error since the following test then fails

 if (lda<1 .or. n<1 .or. ldb<1 .or. ldb/=m .or. ldx/=m) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & 'b=',[ldb,nrhs],' x=',[ldx,nrhsx]) call linalg_error_handling(err0,err) if (present(rank)) rank = 0 return end if 

in e.g. stdlib_linalg_s_solve_lstsq_one. Are there any reasons why we should keep the output vector as an m-vector rather than an n-vector ? From a mathematical and practical point of view, I would expect the following behaviour:

• lstsq takes an m x n matrix A and an m-vector b as entry and returns an n-vector.
• solve_lstsq takes as arguments m x n matrix A, an m-vector b and an n-vector x.