diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 7ed224e21..1a9b852a9 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -860,6 +860,122 @@ This subroutine computes the internal working space requirements for the least-s `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. +## `constrained_lstsq` - Compute the solution of the equality-constrained least-squares problem {#constrained-lstsq} + +### Status + +Experimental + +### Description + +This function computes the solution \(x\) of the equality-constrained linear least-squares problem +$$ +\begin{aligned} + \mathrm{minimize} & \quad \| Ax - b \|^2 \\ + \mathrm{subject~to} & \quad Cx = d, +\end{aligned} +$$ +where \(A\) is an \( m \times n \) matrix (with \(m \geq n\)) and \(C\) a \( p \times n\) matrix (with \(p \leq n\)). The solver is based on LAPACK's `*GLSE` backends. + +### Syntax + +`x = ` [[stdlib_linalg(module):constrained_lstsq(interface)]] `(A, b, C, d[, overwrite_matrices, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(inout)` argument. + +`b`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the least-squares cost. It is an `intent(inout)` argument. + +`c`: Shall be a rank-2 `real` or `complex` array of the same kind as `a` defining the linear equality constraints. It is an `intent(inout)` argument. + +`d`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the linear equality constraints. + +`overwrite_matrices` (optional): Shall be an input `logical` flag. If `.true.`, the input matrices and vectors will be overwritten during the computation of the solution. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns an array of the same kind as `a` containing the solution of the equality constrained least-squares problem. + +Raises `LINALG_ERROR` if the underlying constrained least-squares solver did not converge. +Raises `LINALG_VALUE_ERROR` if the matrices and vectors have invalid/incompatible dimensions. +Exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_constrained_lstsq1.f90!} +``` + +## `solve_constrained_lstsq` - Compute the solution of the equality-constrained least squares problem (subroutine interface) {#solve-constrained-lstsq} + +### Status + +Experimental + +### Description + +This subroutine computes the solution \(x\) of the equality-constrained linear least-squares problem +$$ +\begin{aligned} + \mathrm{minimize} & \quad \| Ax - b \|^2 \\ + \mathrm{subject~to} & \quad Cx = d, +\end{aligned} +$$ +where \(A\) is an \( m \times n \) matrix (with \(m \geq n\)) and \(C\) a \( p \times n\) matrix (with \(p \leq n\)). The solver is based on LAPACK's `*GLSE` backends. + +### Syntax + + +`call ` [[stdlib_linalg(module):solve_constrained_lstsq(interface)]] `(a, b, c, d, x [, storage, overwrite_matrices, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(inout)` argument. + +`b`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the least-squares cost. It is an `intent(inout)` argument. + +`c`: Shall be a rank-2 `real` or `complex` array of the same kind as `a` defining the linear equality constraints. It is an `intent(inout)` argument. + +`d`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the linear equality constraints. + +`x`: Shall be a rank-1 array of the same kind as `a`. On exit, it contains the solution of the constrained least-squares problem. It is an `intent(out)` argument. + +`storage` (optional): Shall a rank-1 array of the same kind as `a` providing working storage for the solver. Its minimum size can be determined with a call to [stdlib_linalg(module):constrained_lstsq_space(interface)]. It is an `intent(out)` argument. + +`overwrite_matrices` (optional): Shall be an input `logical` flag. If `.true.`, the input matrices and vectors will be overwritten during the computation of the solution. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Example + +```fortran +{!example/linalg/example_constrained_lstsq2.f90!} +``` + +## `constrained_lstsq_space` - Compute internal workspace requirements for the constrained least-square solver {#constrained-lstsq-space} + +### Status + +Experimental + +### Description + +This subroutine computes the internal workspace requirements for the constrained least-squares solver, [stdlib_linalg(module):solve_constrained_lstsq(interface)]. + +### Syntax + +call [stdlib_linalg(module):constrained_lstsq_space(interface)]`(a, c, lwork [, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(in)` argument. + +`c`: Shall be a rank-2 `real` or `complex` array of the same kind as `a` defining the linear equality constraints. It is an `intent(in)` argument. +`lwork`: Shall be an `integer` scalar returning the optimal size required for the workspace array to solve the constrained least-squares problem. + ## `det` - Computes the determinant of a square matrix ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index b00cfe1d4..d99d91f6f 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -35,6 +35,8 @@ ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) +ADD_EXAMPLE(constrained_lstsq1) +ADD_EXAMPLE(constrained_lstsq2) ADD_EXAMPLE(norm) ADD_EXAMPLE(mnorm) ADD_EXAMPLE(get_norm) diff --git a/example/linalg/example_constrained_lstsq1.f90 b/example/linalg/example_constrained_lstsq1.f90 new file mode 100644 index 000000000..4884b2790 --- /dev/null +++ b/example/linalg/example_constrained_lstsq1.f90 @@ -0,0 +1,40 @@ +! Constrained least-squares solver: functional interface +program example_constrained_lstsq1 + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: constrained_lstsq + implicit none + integer, parameter :: m = 5, n = 4, p = 3 + !> Least-squares cost. + real(dp) :: A(m, n), b(m) + !> Equality constraints. + real(dp) :: C(p, n), d(p) + !> Solution. + real(dp) :: x(n), x_true(n) + + !> Least-squares cost. + A(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp] + A(2, :) = [1.0_dp, 3.0_dp, 1.0_dp, 1.0_dp] + A(3, :) = [1.0_dp, -1.0_dp, 3.0_dp, 1.0_dp] + A(4, :) = [1.0_dp, 1.0_dp, 1.0_dp, 3.0_dp] + A(5, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp] + + b = [2.0_dp, 1.0_dp, 6.0_dp, 3.0_dp, 1.0_dp] + + !> Equality constraints. + C(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp] + C(2, :) = [1.0_dp, -1.0_dp, 1.0_dp, 1.0_dp] + C(3, :) = [1.0_dp, 1.0_dp, -1.0_dp, 1.0_dp] + + d = [1.0_dp, 3.0_dp, -1.0_dp] + + ! Compute the solution. + x = constrained_lstsq(A, b, C, d) + x_true = [0.5_dp, -0.5_dp, 1.5_dp, 0.5_dp] + + print *, "Exact solution :" + print *, x_true + print * + print *, "Computed solution :" + print *, x + +end program example_constrained_lstsq1 diff --git a/example/linalg/example_constrained_lstsq2.f90 b/example/linalg/example_constrained_lstsq2.f90 new file mode 100644 index 000000000..997181f72 --- /dev/null +++ b/example/linalg/example_constrained_lstsq2.f90 @@ -0,0 +1,48 @@ +! Demonstrate expert subroutine interface with pre-allocated arrays +program example_constrained_lstsq2 + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: solve_constrained_lstsq, constrained_lstsq_space + implicit none + integer, parameter :: m = 5, n = 4, p = 3 + !> Least-squares cost. + real(dp) :: A(m, n), b(m) + !> Equality constraints. + real(dp) :: C(p, n), d(p) + !> Solution. + real(dp) :: x(n), x_true(n) + !> Workspace array. + integer :: lwork + real(dp), allocatable :: work(:) + + !> Least-squares cost. + A(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp] + A(2, :) = [1.0_dp, 3.0_dp, 1.0_dp, 1.0_dp] + A(3, :) = [1.0_dp, -1.0_dp, 3.0_dp, 1.0_dp] + A(4, :) = [1.0_dp, 1.0_dp, 1.0_dp, 3.0_dp] + A(5, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp] + + b = [2.0_dp, 1.0_dp, 6.0_dp, 3.0_dp, 1.0_dp] + + !> Equality constraints. + C(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp] + C(2, :) = [1.0_dp, -1.0_dp, 1.0_dp, 1.0_dp] + C(3, :) = [1.0_dp, 1.0_dp, -1.0_dp, 1.0_dp] + + d = [1.0_dp, 3.0_dp, -1.0_dp] + + !> Optimal workspace size. + call constrained_lstsq_space(A, C, lwork) + allocate (work(lwork)) + + ! Compute the solution. + call solve_constrained_lstsq(A, b, C, d, x, & + storage=work, & + overwrite_matrices=.true.) + x_true = [0.5_dp, -0.5_dp, 1.5_dp, 0.5_dp] + + print *, "Exact solution :" + print *, x_true + print * + print *, "Computed solution :" + print *, x +end program example_constrained_lstsq2 diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 442999ca3..7142302db 100644 --- a/src/lapack/stdlib_linalg_lapack_aux.fypp +++ b/src/lapack/stdlib_linalg_lapack_aux.fypp @@ -51,6 +51,7 @@ module stdlib_linalg_lapack_aux public :: handle_geev_info public :: handle_ggev_info public :: handle_heev_info + public :: handle_gglse_info ! SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments ! used to select eigenvalues to sort to the top left of the Schur form. @@ -1607,4 +1608,28 @@ module stdlib_linalg_lapack_aux end subroutine handle_heev_info + elemental subroutine handle_gglse_info(this, info, m, n, p, err) + character(len=*), intent(in) :: this + integer(ilp), intent(in) :: info, m, n, p + type(linalg_state_type), intent(out) :: err + ! Process output. + select case (info) + case(2) + err = linalg_state_type(this, LINALG_ERROR, "rank([A, B]) < n, the least-squares solution cannot be computed.") + case(1) + err = linalg_state_type(this, LINALG_ERROR, "rank(C) < p, the least-squares solution cannot be computed.") + case(0) + ! Success. + err%state = LINALG_SUCCESS + case(-1) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of rows for A, m=', m) + case(-2) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of columns for A and C, n=', n) + case(-3) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of rows for C, p=', p) + case default + err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'catastrophic error.') + end select + end subroutine handle_gglse_info + end module stdlib_linalg_lapack_aux diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 98fdfa172..fe0fa4dca 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -38,12 +38,15 @@ module stdlib_linalg public :: operator(.pinv.) public :: lstsq public :: lstsq_space + public :: constrained_lstsq + public :: constrained_lstsq_space public :: norm public :: mnorm public :: get_norm public :: solve public :: solve_lu public :: solve_lstsq + public :: solve_constrained_lstsq public :: trace public :: svd public :: svdvals @@ -576,6 +579,106 @@ module stdlib_linalg #:endfor end interface lstsq_space + ! Equality-constrained least-squares, i.e. minimize the sum of squares + ! cost || Ax - b ||^2 subject to the equality constraint Cx = d. + interface constrained_lstsq + !! version: experimental + !! + !! Computes the solution of the equality constrained least-squares problem + !! + !! minimize || Ax - b ||² + !! subject to Cx = d + !! + !! where A is of size `m x n` and C of size `p x n`. + !! ([Specification](../page/specs/stdlib_linalg.html#constrained-lstsq)) + !! + !! ### Description + !! + !! This interface provides methods for computing the solution of an equality-constrained + !! least-squares problem using a function. Supported data types include `real` and + !! `complex`. + !! + !! @note The solution is based on LAPACK's `*GGLSE` methods. + #:for rk, rt, ri in RC_KINDS_TYPES + module function stdlib_linalg_${ri}$_constrained_lstsq(A, b, C, d, overwrite_matrices, err) result(x) + !> Least-squares cost + ${rt}$, intent(inout), target :: A(:, :), b(:) + !> Equality constraints. + ${rt}$, intent(inout), target :: C(:, :), d(:) + !> [optional] Can A and C be overwritten? + logical(lk), optional, intent(in) :: overwrite_matrices + !> [optional] State return flag. + type(linalg_state_type), optional, intent(out) :: err + !> Solution of the constrained least-squares problem. + ${rt}$, allocatable, target :: x(:) + end function stdlib_linalg_${ri}$_constrained_lstsq + #:endfor + end interface + + ! Equality-constrained least-squares, i.e. minimize the sum of squares + ! cost || Ax - b ||^2 subject to the equality constraint Cx = d. + interface solve_constrained_lstsq + !! version: experimental + !! + !! Computes the solution of the equality constrained least-squares problem + !! + !! minimize || Ax - b ||² + !! subject to Cx = d + !! + !! where A is of size `m x n` and C of size `p x n`. + !! ([Specification](../page/specs/stdlib_linalg.html#solve-constrained-lstsq)) + !! + !! ### Description + !! + !! This interface provides methods for computing the solution of an equality-constrained + !! least-squares problem using a subroutine. Supported data types include `real` and + !! `complex`. If a pre-allocated workspace is provided, no internal memory allocation takes + !! place. + !! + !! @note The solution is based on LAPACK's `*GGLSE` methods. + #:for rk, rt, ri in RC_KINDS_TYPES + module subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq(A, b, C, d, x, storage, overwrite_matrices, err) + !> Least-squares cost. + ${rt}$, intent(inout), target :: A(:, :), b(:) + !> Equality constraints. + ${rt}$, intent(inout), target :: C(:, :), d(:) + !> Solution vector. + ${rt}$, intent(out) :: x(:) + !> [optional] Storage. + ${rt}$, optional, intent(out), target :: storage(:) + !> [optional] Can A and C be overwritten? + logical(lk), optional, intent(in) :: overwrite_matrices + !> [optional] State return flag. On error if not requested, the code stops. + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq + #:endfor + end interface + + interface constrained_lstsq_space + !! version: experimental + !! + !! Computes the size of the workspace required by the constrained least-squares solver. + !! ([Specification](../page/specs/stdlib_linalg.html#constrained-lstsq-space)) + !! + !! ### Description + !! + !! This interface provides the size of the workspace array required by the constrained + !! least-squares solver. It can be used to pre-allocate the working array in + !! case several repeated solutions to a same system are sought. If pre-allocated, + !! working arrays are provided, no internal allocation will take place. + !! + #:for rk, rt, ri in RC_KINDS_TYPES + module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, C, lwork, err) + !> Least-squares cost. + ${rt}$, intent(in) :: A(:, :) + !> Equality constraints. + ${rt}$, intent(in) :: C(:, :) + integer(ilp), intent(out) :: lwork + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space + #:endfor + end interface + ! QR factorization of rank-2 array A interface qr !! version: experimental diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index 8cb374689..bee76e40f 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -7,8 +7,8 @@ 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_lapack_aux, only: handle_gelsd_info + use stdlib_linalg_lapack, only: gelsd, gglse, stdlib_ilaenv + use stdlib_linalg_lapack_aux, only: handle_gelsd_info, handle_gglse_info use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR implicit none @@ -60,7 +60,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - ! Compute the integer, real, [complex] working space requested byu the least squares procedure + ! Compute the integer, real, [complex] working space requested by 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(:,:) @@ -170,7 +170,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:if rt.startswith('c') !> [optional] complex working storage space ${rt}$, optional, intent(inout), target :: cmpl_storage(:) - #:endif +#: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 @@ -363,4 +363,205 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares endif end function ilog2 + !------------------------------------------------------------- + !----- Equality-constrained Least-Squares solver ----- + !------------------------------------------------------------- + + pure subroutine check_problem_size(ma, na, mb, mc, nc, md, mx, err) + integer(ilp), intent(in) :: ma, na, mb, mc, nc, md, mx + type(linalg_state_type), intent(out) :: err + + ! Check sizes. + if (ma < 1 .or. na < 1) then + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid matrix size a(m, n) =', [ma, na]) + return + else if (mc < 1 .or. nc < 1) then + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid matrix size c(p, n) =', [mc, nc]) + else if (na /= nc) then + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix A and matrix C have inconsistent number of columns.') + else if (mb /= ma) then + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Size(b) inconsistent with number of rows in a, size(b) =', mb) + else if (md /= mc) then + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Size(d) inconsistent with number of rows in c, size(d) =', md) + else if (na /= mx) then + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Size(x) inconsistent with number of columns of a, size(x) =', mx) + endif + end subroutine check_problem_size + + #:for rk, rt, ri in RC_KINDS_TYPES + ! Compute the size of the workspace requested by the constrained least-squares procedure. + module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, C, lwork, err) + !> Least-squares cost. + ${rt}$, intent(in) :: A(:, :) + !> Equality constraints. + ${rt}$, intent(in) :: C(:, :) + !> Size of the workspace array. + integer(ilp), intent(out) :: lwork + !> [optional] State return flag. + type(linalg_state_type), optional, intent(out) :: err + !> Local variables. + integer(ilp) :: m, n, p, info + ${rt}$ :: a_dummy(1, 1), b_dummy(1) + ${rt}$ :: c_dummy(1, 1), d_dummy(1) + ${rt}$ :: work(1), x(1) + type(linalg_state_type) :: err0 + !> Problem dimensions. + m = size(A, 1) ; n = size(A, 2) ; p = size(C, 1) + lwork = -1_ilp + !> Workspace query. + call gglse(m, n, p, a_dummy, m, c_dummy, p, b_dummy, d_dummy, x, work, lwork, info) + call handle_gglse_info(this, info, m, n, p, err0) + !> Optimal workspace size. + lwork = ceiling(real(work(1), kind=${rk}$), kind=ilp) + + call linalg_error_handling(err0, err) + end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space + + ! Constrained least-squares solver. + module subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq(A, b, C, d, x, storage, overwrite_matrices, err) + !! ### Summary + !! Compute the solution of the equality constrained least-squares problem + !! + !! minimize || Ax - b ||² + !! subject to Cx = d + !! + !! ### Description + !! + !! This function computes the solution of an equality constrained linear least-squares + !! problem. + !! + !! param: a Input matrix of size [m, n] (with m > n). + !! param: b Right-hand side vector of size [m] in the least-squares cost. + !! param: c Input matrix of size [p, n] (with p < n) defining the equality constraints. + !! param: d Right-hand side vector of size [p] in the equality constraints. + !! param: x Vector of size [n] solution to the problem. + !! param: storage [optional] Working array. + !! param: overwrite_matrices [optional] Boolean flag indicating whether the matrices + !! and vectors can be overwritten. + !! param: err [optional] State return flag. + !! + !> Least-squares cost. + ${rt}$, intent(inout), target :: A(:, :), b(:) + !> Equality constraints. + ${rt}$, intent(inout), target :: C(:, :), d(:) + !> Solution vector. + ${rt}$, intent(out) :: x(:) + !> [optional] Storage. + ${rt}$, optional, intent(out), target :: storage(:) + !> [optional] Can A, b, C, and d be overwritten? + logical(lk), optional, intent(in) :: overwrite_matrices + !> [optional] State return flag. + type(linalg_state_type), optional, intent(out) :: err + + ! Local variables. + type(linalg_state_type) :: err0 + integer(ilp) :: ma, na, mb + integer(ilp) :: mc, nc, md + integer(ilp) :: mx + logical(lk) :: overwrite_matrices_ + ${rt}$, pointer :: amat(:, :), bvec(:) + ${rt}$, pointer :: cmat(:, :), dvec(:) + ! LAPACK related. + integer(ilp) :: lwork, info + ${rt}$, pointer :: work(:) + + !> Check dimensions. + ma = size(A, 1, kind=ilp) ; na = size(A, 2, kind=ilp) + mc = size(C, 1, kind=ilp) ; nc = size(C, 2, kind=ilp) + mb = size(b, kind=ilp) ; md = size(d, kind=ilp) ; mx = size(x, kind=ilp) + call check_problem_size(ma, na, mb, mc, nc, md, mx, err0) + if (err0%error()) then + call linalg_error_handling(err0, err) + return + endif + + !> Check if matrices can be overwritten. + overwrite_matrices_ = optval(overwrite_matrices, .false._lk) + + !> Allocate matrices. + if (overwrite_matrices_) then + amat => a + bvec => b + cmat => c + dvec => d + else + allocate(amat(ma, na), source=a) + allocate(bvec(mb), source=b) + allocate(cmat(mc, nc), source=c) + allocate(dvec(md), source=d) + endif + + !> Retrieve workspace size. + call stdlib_linalg_${ri}$_constrained_lstsq_space(A, C, lwork, err0) + + if (err0%ok()) then + !> Workspace. + if (present(storage)) then + work => storage + else + allocate(work(lwork)) + endif + if (size(work, kind=ilp) < lwork) then + err0 = linalg_state_type(this, LINALG_ERROR, 'Insufficient workspace. Should be at least ', lwork) + call linalg_error_handling(err0, err) + return + endif + + !> Compute constrained lstsq solution. + call gglse(ma, na, mc, amat, ma, cmat, mc, bvec, dvec, x, work, lwork, info) + call handle_gglse_info(this, info, ma, na, mc, err0) + + !> Deallocate. + if (.not. present(storage)) deallocate(work) + endif + + if (.not. overwrite_matrices_) then + deallocate(amat, bvec, cmat, dvec) + endif + + call linalg_error_handling(err0, err) + + end subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq + + module function stdlib_linalg_${ri}$_constrained_lstsq(A, b, C, d, overwrite_matrices, err) result(x) + !! ### Summary + !! Compute the solution of the equality constrained least-squares problem + !! + !! minimize || Ax - b ||² + !! subject to Cx = d + !! + !! ### Description + !! + !! This function computes the solution of an equality constrained linear least-squares + !! problem. + !! + !! param: a Input matrix of size [m, n] (with m > n). + !! param: b Right-hand side vector of size [m] in the least-squares cost. + !! param: c Input matrix of size [p, n] (with p < n) defining the equality constraints. + !! param: d Right-hand side vector of size [p] in the equality constraints. + !! param: x Vector of size [n] solution to the problem. + !! param: overwrite_matrices [optional] Boolean flag indicating whether the matrices + !! and vectors can be overwritten. + !! param: err [optional] State return flag. + !! + !> Least-squares cost. + ${rt}$, intent(inout), target :: A(:, :), b(:) + !> Equality constraints. + ${rt}$, intent(inout), target :: C(:, :), d(:) + !> [optional] Can A, b, C, d be overwritten? + logical(lk), optional, intent(in) :: overwrite_matrices + !> [optional] State return flag. + type(linalg_state_type), optional, intent(out) :: err + !> Solution of the constrained least-squares problem. + ${rt}$, allocatable, target :: x(:) + + ! Local variables. + integer(ilp) :: n + + n = size(A, 2, kind=ilp) + allocate(x(n)) + call stdlib_linalg_${ri}$_solve_constrained_lstsq(A, b, C, d, x, overwrite_matrices=overwrite_matrices, err=err) + end function stdlib_linalg_${ri}$_constrained_lstsq + #:endfor + end submodule stdlib_linalg_least_squares diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index a03c031ec..5f502dd1a 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -6,6 +6,7 @@ set( "test_linalg_inverse.fypp" "test_linalg_pseudoinverse.fypp" "test_linalg_lstsq.fypp" + "test_linalg_constrained_lstsq.fypp" "test_linalg_norm.fypp" "test_linalg_mnorm.fypp" "test_linalg_determinant.fypp" @@ -41,6 +42,7 @@ ADDTEST(linalg_norm) ADDTEST(linalg_mnorm) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) +ADDTEST(linalg_constrained_lstsq) ADDTEST(linalg_qr) ADDTEST(linalg_schur) ADDTEST(linalg_solve_iterative) diff --git a/test/linalg/test_linalg_constrained_lstsq.fypp b/test/linalg/test_linalg_constrained_lstsq.fypp new file mode 100644 index 000000000..215663362 --- /dev/null +++ b/test/linalg/test_linalg_constrained_lstsq.fypp @@ -0,0 +1,154 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test least squares solver +module test_linalg_constrained_least_squares + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: constrained_lstsq, solve_constrained_lstsq, constrained_lstsq_space + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + private + + public :: test_constrained_least_squares + + contains + + !> Solve sample least squares problems + subroutine test_constrained_least_squares(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in REAL_KINDS_TYPES + call add_test(tests,new_unittest("constrained_least_squares_randm_${ri}$",test_constrained_lstsq_random_${ri}$)) + #:endfor + + end subroutine test_constrained_least_squares + + #:for rk,rt,ri in REAL_KINDS_TYPES + !> Fit from random array + subroutine test_constrained_lstsq_random_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + type(linalg_state_type) :: state + integer(ilp), parameter :: m=5, n=4, p=3 + !> Least-squares cost. + ${rt}$ :: A(m, n), b(m) + !> Equality constraints. + ${rt}$ :: C(p, n), d(p) + !> Solution. + ${rt}$ :: x(n), x_true(n) + !> Workspace. + integer(ilp) :: lwork + ${rt}$, allocatable :: work(:) + + !> Least-squares cost. + A(1, :) = [1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$] + A(2, :) = [1.0_${rk}$, 3.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$] + A(3, :) = [1.0_${rk}$, -1.0_${rk}$, 3.0_${rk}$, 1.0_${rk}$] + A(4, :) = [1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$, 3.0_${rk}$] + A(5, :) = [1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$, -1.0_${rk}$] + + b = [2.0_${rk}$, 1.0_${rk}$, 6.0_${rk}$, 3.0_${rk}$, 1.0_${rk}$] + + !> Equality constraints. + C(1, :) = [1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$, -1.0_${rk}$] + C(2, :) = [1.0_${rk}$, -1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$] + C(3, :) = [1.0_${rk}$, 1.0_${rk}$, -1.0_${rk}$, 1.0_${rk}$] + + d = [1.0_${rk}$, 3.0_${rk}$, -1.0_${rk}$] + + !----- Function interface ----- + x = constrained_lstsq(A, b, C, d, err=state) + x_true = [0.5_${rk}$, -0.5_${rk}$, 1.5_${rk}$, 0.5_${rk}$] + + call check(error, state%ok(), state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-x_true) < 1.0e-4_${rk}$), 'Solver converged') + if (allocated(error)) return + + !----- Subroutine interface ----- + call solve_constrained_lstsq(A, b, C, d, x, err=state) + + call check(error, state%ok(), state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-x_true) < 1.0e-4_${rk}$), 'Solver converged') + if (allocated(error)) return + + !----- Pre-allocated storage ----- + call constrained_lstsq_space(A, C, lwork, err=state) + call check(error, state%ok(), state%print()) + if (allocated(error)) return + allocate(work(lwork)) + call solve_constrained_lstsq(A, b, C, d, x, storage=work, err=state) + + call check(error, state%ok(), state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-x_true) < 1.0e-4_${rk}$), 'Solver converged') + if (allocated(error)) return + + !----- Overwrite matrices (performance) ----- + call solve_constrained_lstsq(A, b, C, d, x, storage=work, overwrite_matrices=.true., err=state) + + call check(error, state%ok(), state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-x_true) < 1.0e-4_${rk}$), 'Solver converged') + if (allocated(error)) return + + end subroutine test_constrained_lstsq_random_${ri}$ + + #:endfor + + ! gcc-15 bugfix utility + subroutine add_test(tests,new_test) + type(unittest_type), allocatable, intent(inout) :: tests(:) + type(unittest_type), intent(in) :: new_test + + integer :: n + type(unittest_type), allocatable :: new_tests(:) + + if (allocated(tests)) then + n = size(tests) + else + n = 0 + end if + + allocate(new_tests(n+1)) + if (n>0) new_tests(1:n) = tests(1:n) + new_tests(1+n) = new_test + call move_alloc(from=new_tests,to=tests) + + end subroutine add_test + +end module test_linalg_constrained_least_squares + +program test_constrained_lstsq + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_constrained_least_squares, only : test_constrained_least_squares + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_constrained_least_squares", test_constrained_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_constrained_lstsq