Equality-constrained least-squares solver

doc/specs/stdlib_linalg.md

@@ -860,6 +860,127 @@ 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.
+
+`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(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.
+
+`d`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the linear equality constraints.
+
+`lwork`: Shall be an `interger` 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

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)

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

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

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

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 `*GLLSE` 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 `*GLLSE` 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