diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 89bedaa42..a081291d0 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -600,6 +600,130 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l {!example/linalg/example_is_hessenberg.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 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index f515af446..729f117d3 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -18,6 +18,7 @@ ADD_EXAMPLE(state1) ADD_EXAMPLE(state2) ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) +ADD_EXAMPLE(lstsq1) +ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) - diff --git a/example/linalg/example_lstsq1.f90 b/example/linalg/example_lstsq1.f90 new file mode 100644 index 000000000..dfbc031ea --- /dev/null +++ b/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 diff --git a/example/linalg/example_lstsq2.f90 b/example/linalg/example_lstsq2.f90 new file mode 100644 index 000000000..389b20d4c --- /dev/null +++ b/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 diff --git a/include/common.fypp b/include/common.fypp index 1683239e4..ed1cf2b4e 100644 --- a/include/common.fypp +++ b/include/common.fypp @@ -27,11 +27,20 @@ #:set REAL_KINDS = REAL_KINDS + ["qp"] #:endif +#! BLAS/LAPACK initials for each real kind +#:set REAL_INIT = ["s", "d"] +#:if WITH_XDP +#:set REAL_INIT = REAL_INIT + ["x"] +#:endif +#:if WITH_QP +#:set REAL_INIT = REAL_INIT + ["q"] +#:endif + #! Real types to be considered during templating #:set REAL_TYPES = ["real({})".format(k) for k in REAL_KINDS] #! Collected (kind, type) tuples for real types -#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES)) +#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_INIT)) #! Complex kinds to be considered during templating #:set CMPLX_KINDS = ["sp", "dp"] @@ -42,11 +51,20 @@ #:set CMPLX_KINDS = CMPLX_KINDS + ["qp"] #:endif +#! BLAS/LAPACK initials for each complex kind +#:set CMPLX_INIT = ["c", "z"] +#:if WITH_XDP +#:set CMPLX_INIT = CMPLX_INIT + ["y"] +#:endif +#:if WITH_QP +#:set CMPLX_INIT = CMPLX_INIT + ["w"] +#:endif + #! Complex types to be considered during templating #:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS] -#! Collected (kind, type) tuples for complex types -#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES)) +#! Collected (kind, type, initial) tuples for complex types +#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_INIT)) #! Integer kinds to be considered during templating #:set INT_KINDS = ["int8", "int16", "int32", "int64"] @@ -109,6 +127,17 @@ #{if rank > 0}#(${":" + ",:" * (rank - 1)}$)#{endif}# #:enddef +#! Generates an empty array rank suffix. +#! +#! Args: +#! rank (int): Rank of the variable +#! +#! Returns: +#! Empty array rank suffix string (e.g. (0,0) if rank = 2) +#! +#:def emptyranksuffix(rank) +#{if rank > 0}#(${"0" + ",0" * (rank - 1)}$)#{endif}# +#:enddef #! Joins stripped lines with given character string #! diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 38d0ab569..b58ba0b9c 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 2f837f15a..460aa0593 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1,11 +1,15 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES #:set RCI_KINDS_TYPES = RC_KINDS_TYPES + INT_KINDS_TYPES +#:set RHS_SUFFIX = ["one","many"] +#:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]] +#:set RHS_EMPTY = [emptyranksuffix(r) for r in [1,2]] +#:set ALL_RHS = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY)) module stdlib_linalg !!Provides a support for various linear algebra procedures !! ([Specification](../page/specs/stdlib_linalg.html)) - use stdlib_kinds, only: sp, dp, xdp, qp, lk, & - int8, int16, int32, int64 + use stdlib_kinds, only: xdp, int8, int16, int32, int64 + use stdlib_linalg_constants, only: sp, dp, qp, lk, ilp use stdlib_error, only: error_stop use stdlib_optval, only: optval use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling @@ -16,6 +20,9 @@ module stdlib_linalg public :: operator(.det.) public :: diag public :: eye + public :: lstsq + public :: lstsq_space + public :: solve_lstsq public :: trace public :: outer_product public :: kronecker_product @@ -221,6 +228,131 @@ module stdlib_linalg #:endfor end interface is_hessenberg + ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. + interface lstsq + !! version: experimental + !! + !! Computes the squares solution to system \( A \cdot x = b \). + !! ([Specification](../page/specs/stdlib_linalg.html#lstsq-computes-the-least-squares-solution-to-a-linear-matrix-equation)) + !! + !!### Summary + !! Interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution. + !! + !!### Description + !! + !! This interface provides methods for computing the least squares of a linear matrix system. + !! Supported data types include `real` and `complex`. + !! + !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) + !> Input matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, allocatable, target :: x${nd}$ + end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface lstsq + + ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. + interface solve_lstsq + !! version: experimental + !! + !! Computes the squares solution to system \( A \cdot x = b \). + !! ([Specification](../page/specs/stdlib_linalg.html#solve-lstsq-compute-the-least-squares-solution-to-a-linear-matrix-equation-subroutine-interface)) + !! + !!### Summary + !! Subroutine interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution. + !! + !!### Description + !! + !! This interface provides methods for computing the least squares of a linear matrix system using + !! a subroutine. Supported data types include `real` and `complex`. If pre-allocated work spaces + !! are provided, no internal memory allocations take place when using this interface. + !! + !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage,& + #{if rt.startswith('c')}#cmpl_storage,#{endif}#cond,singvals,overwrite_a,rank,err) + !> Input matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] real working storage space + real(${rk}$), optional, intent(inout), target :: real_storage(:) + !> [optional] integer working storage space + integer(ilp), optional, intent(inout), target :: int_storage(:) + #:if rt.startswith('complex') + !> [optional] complex working storage space + ${rt}$, optional, intent(inout), target :: cmpl_storage(:) + #:endif + !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] list of singular values [min(m,n)], in descending magnitude order, returned by the SVD + real(${rk}$), optional, intent(out), target :: singvals(:) + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface solve_lstsq + + ! Return the working array space required by the least squares solver + interface lstsq_space + !! version: experimental + !! + !! Computes the integer, real [, complex] working space required by the least-squares solver + !! ([Specification](../page/specs/stdlib_linalg.html#lstsq-space-compute-internal-working-space-requirements-for-the-least-squares-solver)) + !! + !!### Description + !! + !! This interface provides sizes of integer, real [, complex] working spaces required by the + !! least-squares solver. These sizes can be used to pre-allocated working arrays in case several + !! repeated least-squares solutions to a same system are sought. If pre-allocated working arrays + !! are provided, no internal allocations will take place. + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Size of the working space arrays + integer(ilp), intent(out) :: lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}# + end subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface lstsq_space interface det !! version: experimental diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp new file mode 100644 index 000000000..0b7c2e139 --- /dev/null +++ b/src/stdlib_linalg_least_squares.fypp @@ -0,0 +1,353 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set RHS_SUFFIX = ["one","many"] +#:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]] +#:set RHS_EMPTY = [emptyranksuffix(r) for r in [1,2]] +#:set ALL_RHS = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY)) +submodule (stdlib_linalg) stdlib_linalg_least_squares +!! Least-squares solution to Ax=b + use stdlib_linalg_constants + use stdlib_linalg_lapack, only: gelsd, stdlib_ilaenv + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + implicit none(type,external) + + character(*), parameter :: this = 'lstsq' + + contains + + elemental subroutine handle_gelsd_info(info,lda,n,ldb,nrhs,err) + integer(ilp), intent(in) :: info,lda,n,ldb,nrhs + type(linalg_state_type), intent(out) :: err + + ! Process output + select case (info) + case (0) + ! Success + case (:-1) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=',[lda,n], & + ', b=',[ldb,nrhs]) + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select + + end subroutine handle_gelsd_info + + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + ! Workspace needed by gelsd + elemental subroutine ${ri}$gelsd_space(m,n,nrhs,lrwork,liwork,lcwork) + integer(ilp), intent(in) :: m,n,nrhs + integer(ilp), intent(out) :: lrwork,liwork,lcwork + + integer(ilp) :: smlsiz,mnmin,nlvl + + mnmin = min(m,n) + + ! Maximum size of the subproblems at the bottom of the computation (~25) + smlsiz = stdlib_ilaenv(9,'${ri}$gelsd',' ',0,0,0,0) + + ! The exact minimum amount of workspace needed depends on M, N and NRHS. + nlvl = max(0, ilog2(mnmin/(smlsiz+1))+1) + + ! Real space + #:if rt.startswith('complex') + lrwork = 10*mnmin+2*mnmin*smlsiz+8*mnmin*nlvl+3*smlsiz*nrhs+max((smlsiz+1)**2,n*(1+nrhs)+2*nrhs) + #:else + lrwork = 12*mnmin+2*mnmin*smlsiz+8*mnmin*nlvl+mnmin*nrhs+(smlsiz+1)**2 + #:endif + lrwork = max(1,lrwork) + + ! Complex space + lcwork = 2*mnmin + nrhs*mnmin + + ! Integer space + liwork = max(1, 3*mnmin*nlvl+11*mnmin) + + ! For good performance, the workspace should generally be larger. + ! Allocate 25% more space than strictly needed. + lrwork = ceiling(1.25*lrwork,kind=ilp) + lcwork = ceiling(1.25*lcwork,kind=ilp) + liwork = ceiling(1.25*liwork,kind=ilp) + + end subroutine ${ri}$gelsd_space + + #:endif + #:endfor + + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + + ! Compute the integer, real, [complex] working space requested byu the least squares procedure + pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Size of the working space arrays + integer(ilp), intent(out) :: lrwork,liwork + integer(ilp) #{if rt.startswith('c')}#, intent(out)#{endif}# :: lcwork + + integer(ilp) :: m,n,nrhs + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + nrhs = size(b,kind=ilp)/size(b,1,kind=ilp) + + call ${ri}$gelsd_space(m,n,nrhs,lrwork,liwork,lcwork) + + end subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$ + + ! Compute the least-squares solution to a real system of linear equations Ax = b + module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) + !!### Summary + !! Compute least-squares solution to a real system of linear equations \( Ax = b \) + !! + !!### Description + !! + !! This function computes the least-squares solution of a linear matrix problem. + !! + !! param: a Input matrix of size [m,n]. + !! param: b Right-hand-side vector of size [n] or matrix of size [n,nrhs]. + !! param: cond [optional] Real input threshold indicating that singular values `s_i <= cond*maxval(s)` + !! do not contribute to the matrix rank. + !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten. + !! param: rank [optional] integer flag returning matrix rank. + !! param: err [optional] State return flag. + !! return: x Solution vector of size [n] or solution matrix of size [n,nrhs]. + !! + !> Input matrix a[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, allocatable, target :: x${nd}$ + + ! Initialize solution with the shape of the rhs + allocate(x,mold=b) + + call stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,& + cond=cond,overwrite_a=overwrite_a,rank=rank,err=err) + + end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + + ! Compute the least-squares solution to a real system of linear equations Ax = b + module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage, & + #{if rt.startswith('c')}#cmpl_storage,#{endif}# cond,singvals,overwrite_a,rank,err) + + !!### Summary + !! Compute least-squares solution to a real system of linear equations \( Ax = b \) + !! + !!### Description + !! + !! This function computes the least-squares solution of a linear matrix problem. + !! + !! param: a Input matrix of size [m,n]. + !! param: b Right-hand-side vector of size [n] or matrix of size [n,nrhs]. + !! param: x Solution vector of size [n] or solution matrix of size [n,nrhs]. + !! param: real_storage [optional] Real working space + !! param: int_storage [optional] Integer working space + #:if rt.startswith('c') + !! param: cmpl_storage [optional] Complex working space + #:endif + !! param: cond [optional] Real input threshold indicating that singular values `s_i <= cond*maxval(s)` + !! do not contribute to the matrix rank. + !! param: singvals [optional] Real array of size [min(m,n)] returning a list of singular values. + !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten. + !! param: rank [optional] integer flag returning matrix rank. + !! param: err [optional] State return flag. + !! + !> Input matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] real working storage space + real(${rk}$), optional, intent(inout), target :: real_storage(:) + !> [optional] integer working storage space + integer(ilp), optional, intent(inout), target :: int_storage(:) + #:if rt.startswith('c') + !> [optional] complex working storage space + ${rt}$, optional, intent(inout), target :: cmpl_storage(:) + #:endif + !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] list of singular values [min(m,n)], in descending magnitude order, returned by the SVD + real(${rk}$), optional, intent(out), target :: singvals(:) + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + !! Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: m,n,lda,ldb,nrhs,ldx,nrhsx,info,mnmin,mnmax,arank,lrwork,liwork,lcwork + integer(ilp) :: nrs,nis,ncs,nsvd + integer(ilp), pointer :: iwork(:) + logical(lk) :: copy_a + real(${rk}$) :: acond,rcond + real(${rk}$), pointer :: rwork(:),singular(:) + ${rt}$, pointer :: xmat(:,:),amat(:,:),cwork(:) + + ! Problem sizes + m = size(a,1,kind=ilp) + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + ldb = size(b,1,kind=ilp) + nrhs = size(b ,kind=ilp)/ldb + ldx = size(x,1,kind=ilp) + nrhsx = size(x ,kind=ilp)/ldx + mnmin = min(m,n) + mnmax = max(m,n) + + 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 + + ! Can A be overwritten? By default, do not overwrite + if (present(overwrite_a)) then + copy_a = .not.overwrite_a + else + copy_a = .true._lk + endif + + ! Initialize a matrix temporary + if (copy_a) then + allocate(amat(lda,n),source=a) + else + amat => a + endif + + ! Initialize solution with the rhs + x = b + xmat(1:n,1:nrhs) => x + + ! Singular values array (in decreasing order) + if (present(singvals)) then + singular => singvals + nsvd = size(singular,kind=ilp) + else + allocate(singular(mnmin)) + nsvd = mnmin + endif + + ! rcond is used to determine the effective rank of A. + ! Singular values S(i) <= RCOND*maxval(S) are treated as zero. + ! Use same default value as NumPy + if (present(cond)) then + rcond = cond + else + rcond = epsilon(0.0_${rk}$)*mnmax + endif + if (rcond<0) rcond = epsilon(0.0_${rk}$)*mnmax + + ! Get working space size + call ${ri}$gelsd_space(m,n,nrhs,lrwork,liwork,lcwork) + + ! Real working space + if (present(real_storage)) then + rwork => real_storage + else + allocate(rwork(lrwork)) + endif + nrs = size(rwork,kind=ilp) + + ! Integer working space + if (present(int_storage)) then + iwork => int_storage + else + allocate(iwork(liwork)) + endif + nis = size(iwork,kind=ilp) + + #:if rt.startswith('complex') + ! Complex working space + if (present(cmpl_storage)) then + cwork => cmpl_storage + else + allocate(cwork(lcwork)) + endif + ncs = size(iwork,kind=ilp) + #:endif + + if (nrs=',lrwork, & + ', int=',nis,' should be >=',liwork, & + #{if rt.startswith('complex')}#', cmplx=',ncs,' should be >=',lcwork, &#{endif}# + ', singv=',nsvd,' should be >=',mnmin) + + else + + ! Solve system using singular value decomposition + call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank, & + #:if rt.startswith('complex') + cwork,nrs,rwork,iwork,info) + #:else + rwork,nrs,iwork,info) + #:endif + + ! The condition number of A in the 2-norm = S(1)/S(min(m,n)). + acond = singular(1)/singular(mnmin) + + ! Process output + call handle_gelsd_info(info,lda,n,ldb,nrhs,err0) + + endif + + ! Process output and return + if (copy_a) deallocate(amat) + if (present(rank)) rank = arank + if (.not.present(real_storage)) deallocate(rwork) + if (.not.present(int_storage)) deallocate(iwork) + #:if rt.startswith('complex') + if (.not.present(cmpl_storage)) deallocate(cwork) + #:endif + if (.not.present(singvals)) deallocate(singular) + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$ + + #:endif + #:endfor + #:endfor + + ! Simple integer log2 implementation + elemental integer(ilp) function ilog2(x) + integer(ilp), intent(in) :: x + + integer(ilp) :: remndr + + if (x>0) then + remndr = x + ilog2 = -1_ilp + do while (remndr>0) + ilog2 = ilog2 + 1_ilp + remndr = shiftr(remndr,1) + end do + else + ilog2 = -huge(0_ilp) + endif + end function ilog2 + +end submodule stdlib_linalg_least_squares diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 92006323e..5a07b5bc4 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -2,6 +2,7 @@ set( fppFiles "test_linalg.fypp" "test_blas_lapack.fypp" + "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" "test_linalg_matrix_property_checks.fypp" ) @@ -10,4 +11,5 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_determinant) ADDTEST(linalg_matrix_property_checks) +ADDTEST(linalg_lstsq) ADDTEST(blas_lapack) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp new file mode 100644 index 000000000..29cd07519 --- /dev/null +++ b/test/linalg/test_linalg_lstsq.fypp @@ -0,0 +1,122 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test least squares solver +module test_linalg_least_squares + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: lstsq + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + private + + public :: test_least_squares + + contains + + !> Solve sample least squares problems + subroutine test_least_squares(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + tests = [tests,new_unittest("least_squares_${ri}$",test_lstsq_one_${ri}$), & + new_unittest("least_squares_randm_${ri}$",test_lstsq_random_${ri}$)] + #:endif + #:endfor + + end subroutine test_least_squares + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + !> Simple polynomial fit + subroutine test_lstsq_one_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp) :: rank + + !> Example scattered data + ${rt}$, parameter :: x(*) = real([1.0, 2.5, 3.5, 4.0, 5.0, 7.0, 8.5], ${rk}$) + ${rt}$, parameter :: y(*) = real([0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6], ${rk}$) + ${rt}$, parameter :: ab(*) = real([0.20925829, 0.12013861], ${rk}$) + + ${rt}$ :: M(size(x),2),p(2) + + ! Coefficient matrix for polynomial y = a + b*x**2 + M(:,1) = x**0 + M(:,2) = x**2 + + ! Find polynomial + p = lstsq(M,y,rank=rank,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(p-ab)<1.0e-6_${rk}$), 'data converged') + if (allocated(error)) return + + call check(error, rank==2, 'matrix rank == 2') + if (allocated(error)) return + + end subroutine test_lstsq_one_${ri}$ + + !> Fit from random array + subroutine test_lstsq_random_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp), parameter :: n = 12, m = 3 + ${rt}$ :: xsol(m),x(m),y(n),A(n,m) + + ! Random coefficient matrix and solution + call random_number(A) + call random_number(xsol) + + ! Compute rhs + y = matmul(A,xsol) + + ! Find polynomial + x = lstsq(A,y,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-xsol)<1.0e-6_${rk}$), 'data converged') + if (allocated(error)) return + + end subroutine test_lstsq_random_${ri}$ + + #:endif + #:endfor + +end module test_linalg_least_squares + +program test_lstsq + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_least_squares, only : test_least_squares + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_least_squares", test_least_squares) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_lstsq