Merge branch 'master' into linalg_solve

Author: perazz

doc/specs/stdlib_linalg.md

@@ -625,9 +625,7 @@ Expert interface:
 `x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b [, overwrite_a], err)`
 
 ### Arguments
-
-Two 
-
+ 
 `a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. 
 
 `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument.
@@ -664,6 +662,7 @@ This subroutine computes the solution to a linear matrix equation \( A \cdot x =
 
 Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. 
 An error is returned if the matrix is rank-deficient or singular to working precision. 
+If all optional arrays are provided by the user, no internal allocations take place.
 The solver is based on LAPACK's `*GESV` backends.
 
 ### Syntax
@@ -678,8 +677,6 @@ Expert (`Pure`) interface:
 
 ### Arguments
 
-Two 
-
 `a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. 
 
 `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument.
@@ -690,8 +687,6 @@ Two
 
 `overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
 
-`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
-
 ### Return value
 
 For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.
@@ -703,11 +698,133 @@ If `err` is not present, exceptions trigger an `error stop`.
 ### Example
 
 ```fortran
-{!example/linalg/example_solve1.f90!}
+{!example/linalg/example_solve3.f90!}
+```
 
-{!example/linalg/example_solve2.f90!}
+## `lstsq` - Computes the least squares solution to a linear matrix equation. 
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes the least-squares solution to a linear matrix equation \( A \cdot x = b \).
+
+Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends.
+
+### Syntax
+
+`x = ` [[stdlib_linalg(module):lstsq(interface)]] `(a, b, [, cond, overwrite_a, rank, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument.
+
+`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. 
+
+`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument.
+
+`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `A` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
+
+`rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `A`. This is an `intent(out)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
+
+### Return value
+
+Returns an array value of the same kind and rank as `b`, containing the solution(s) to the least squares system. 
+
+Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
+Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_lstsq1.f90!}
 ```
 
+## `solve_lstsq` - Compute the least squares solution to a linear matrix equation (subroutine interface). 
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the least-squares solution to a linear matrix equation \( A \cdot x = b \).
+
+Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):solve_lstsq(interface)]] `(a, b, x, [, real_storage, int_storage, [cmpl_storage, ] cond, singvals, overwrite_a, rank, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument.
+
+`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. 
+
+`x`: Shall be an array of same kind and rank as `b`, containing the solution(s) to the least squares system. It is an `intent(inout)` argument.
+
+`real_storage` (optional): Shall be a `real` rank-1 array of the same kind `a`, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument.
+
+`int_storage` (optional): Shall be an `integer` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument.
+
+`cmpl_storage` (optional): For `complex` systems, it shall be a `complex` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument.
+
+`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument.
+
+`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `minval(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` argument.
+
+`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `A` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
+
+`rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `A`. This is an `intent(out)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
+
+### Return value
+
+Returns an array value that represents the solution to the least squares system.
+
+Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
+Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_lstsq2.f90!}
+```
+
+## `lstsq_space` - Compute internal working space requirements for the least squares solver.
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the internal working space requirements for the least-squares solver, [[stdlib_linalg(module):solve_lstsq(interface)]] .
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):lstsq_space(interface)]] `(a, b, lrwork, liwork [, lcwork])`
+
+### Arguments
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the linear system coefficient matrix. It is an `intent(in)` argument.
+
+`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the system's right-hand-side vector(s). It is an `intent(in)` argument. 
+
+`lrwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `real` working storage to this system.
+
+`liwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `integer` working storage to this system.
+
+`lcwork` (`complex` `a`, `b`): For a `complex` system, shall be an `integer` scalar, that returns the minimum array size required for the `complex` working storage to this system.
+
 ## `det` - Computes the determinant of a square matrix
 
 ### Status

example/linalg/CMakeLists.txt

@@ -18,6 +18,8 @@ ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)
 ADD_EXAMPLE(blas_gemv)
 ADD_EXAMPLE(lapack_getrf)
+ADD_EXAMPLE(lstsq1)
+ADD_EXAMPLE(lstsq2)
 ADD_EXAMPLE(solve1)
 ADD_EXAMPLE(solve2)
 ADD_EXAMPLE(solve3)

example/linalg/example_lstsq1.f90

@@ -0,0 +1,26 @@
+! Least-squares solver: functional interface
+program example_lstsq1
+  use stdlib_linalg_constants, only: dp
+  use stdlib_linalg, only: lstsq
+  implicit none
+
+  integer, allocatable :: x(:),y(:)
+  real(dp), allocatable :: A(:,:),b(:),coef(:)
+
+  ! Data set
+  x = [1, 2, 2]
+  y = [5, 13, 25]
+
+  ! Fit three points using a parabola, least squares method
+  ! A = [1 x x**2]
+  A = reshape([[1,1,1],x,x**2],[3,3])
+  b = y
+
+  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
+  coef = lstsq(A,b)
+
+  print *, 'parabola: ',coef
+  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 
+
+
+end program example_lstsq1

example/linalg/example_lstsq2.f90

@@ -0,0 +1,40 @@
+! Demonstrate expert subroutine interface with pre-allocated arrays
+program example_lstsq2
+  use stdlib_linalg_constants, only: dp,ilp
+  use stdlib_linalg, only: solve_lstsq, lstsq_space, linalg_state_type
+  implicit none
+
+  integer, allocatable :: x(:),y(:)
+  real(dp), allocatable :: A(:,:),b(:),coef(:),real_space(:),singvals(:)
+  integer(ilp), allocatable :: int_space(:)
+  integer(ilp) :: lrwork,liwork,arank
+
+  ! Data set
+  x = [1, 2, 2]
+  y = [5, 13, 25]
+
+  ! Fit three points using a parabola, least squares method
+  ! A = [1 x x**2]
+  A = reshape([[1,1,1],x,x**2],[3,3])
+  b = y
+  
+  ! Get storage sizes for the arrays and pre-allocate data
+  call lstsq_space(A,b,lrwork,liwork)
+  allocate(coef(size(x)),real_space(lrwork),int_space(liwork),singvals(minval(shape(A))))
+  
+  ! Solve coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
+  ! with no internal allocations
+  call solve_lstsq(A,b,x=coef,              &
+                   real_storage=real_space, &
+                   int_storage=int_space,   &
+                   singvals=singvals,       &
+                   overwrite_a=.true.,      &
+                   rank=arank)
+  
+  print *, 'parabola: ',coef
+  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 
+  print *, 'rank: ',arank
+  ! rank: 2
+
+
+end program example_lstsq2

src/CMakeLists.txt

@@ -21,6 +21,7 @@ set(fppFiles
     stdlib_kinds.fypp
     stdlib_linalg.fypp
     stdlib_linalg_diag.fypp
+    stdlib_linalg_least_squares.fypp
     stdlib_linalg_outer_product.fypp
     stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp