From 1e02559c39a637ced5a897aa8cafcd7f324c24af Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Wed, 22 Oct 2025 14:31:01 +0200
Subject: [PATCH 01/17] Interfaces for the constrained lstsq solver.

---
 src/stdlib_linalg.fypp | 50 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 98fdfa172..5597eb239 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -576,6 +576,56 @@ 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
+    #:for rk, rt, ri in RC_KINDS_TYPES
+    module function stdlib_linalg_${ri}$_constrained_lstsq(A, b, C, d, overwrite_matrices, err) result(x)
+        !> Input matrices.
+        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
+        !> Right-hand side vectors.
+        ${rt}$, intent(inout), target :: b(:), 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
+    #: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)
+        !> Input matrices.
+        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
+        !> Right-hand side vectors.
+        ${rt}$, intent(inout), target :: b(:), d(:)
+        !> Solution vector.
+        ${rt}$, intent(out) :: x(:)
+        !> [optional] Storage.
+        ${rt}$, optional, intent(out) :: 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
+    #:for rk, rt, ri in RC_KINDS_TYPES
+    pure module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, b, C, d, lwork)
+        ${rt}$, intent(in), target :: A(:, :), C(:, :)
+        ${rt}$, intent(in), target :: b(:), d(:)
+        integer(ilp), intent(out) :: lwork
+    end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space
+    #:endfor
+  end interface
+
   ! QR factorization of rank-2 array A
   interface qr
     !! version: experimental 

From a8873d9bf548e5cb3876e85ff0a8cabe04ce36b9 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Wed, 22 Oct 2025 15:21:12 +0200
Subject: [PATCH 02/17] Info handling function for gglse.

---
 src/lapack/stdlib_linalg_lapack_aux.fypp | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp
index 442999ca3..5670bfc95 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,27 @@ 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.
+            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

From 728d221923654c22c38ea04195801e92d6a8774e Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Wed, 22 Oct 2025 17:20:25 +0200
Subject: [PATCH 03/17] Compilable implementation.

---
 src/stdlib_linalg.fypp               |   8 +-
 src/stdlib_linalg_least_squares.fypp | 144 ++++++++++++++++++++++++++-
 2 files changed, 147 insertions(+), 5 deletions(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 5597eb239..11bfb752d 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
@@ -607,7 +610,7 @@ module stdlib_linalg
         !> Solution vector.
         ${rt}$, intent(out) :: x(:)
         !> [optional] Storage.
-        ${rt}$, optional, intent(out) :: 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.
@@ -618,10 +621,11 @@ module stdlib_linalg
 
   interface constrained_lstsq_space
     #:for rk, rt, ri in RC_KINDS_TYPES
-    pure module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, b, C, d, lwork)
+    module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, b, C, d, lwork, err)
         ${rt}$, intent(in), target :: A(:, :), C(:, :)
         ${rt}$, intent(in), target :: b(:), d(:)
         integer(ilp), intent(out) :: lwork
+        type(linalg_state_type), optional, intent(out) :: err
     end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space
     #:endfor
   end interface
diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp
index 8cb374689..6a51b1cce 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
@@ -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,142 @@ 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(m, 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
+    module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, b, C, d, lwork, err)
+        ${rt}$, intent(in), target :: A(:, :), C(:, :)
+        ${rt}$, intent(in), target :: b(:), d(:)
+        integer(ilp), intent(out)  :: lwork
+        type(linalg_state_type), optional, intent(out) :: err
+    end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space
+
+    module subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq(A, b, C, d, x, storage, overwrite_matrices, err)
+        !> Input matrices.
+        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
+        !> Right-hand side vectors.
+        ${rt}$, intent(inout), target :: b(:), 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, b, C, d, 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.
+            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)
+        !> Input matrices.
+        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
+        !> Right-hand side vectors.
+        ${rt}$, intent(inout), target :: b(:), 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

From 36c17fb182db92fb2209f56e38e38d736251e8c5 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 06:58:53 +0200
Subject: [PATCH 04/17] Working implementation of constrained lstsq.

---
 src/stdlib_linalg_least_squares.fypp | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp
index 6a51b1cce..ce7b6c697 100644
--- a/src/stdlib_linalg_least_squares.fypp
+++ b/src/stdlib_linalg_least_squares.fypp
@@ -394,6 +394,19 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
         ${rt}$, intent(in), target :: b(:), d(:)
         integer(ilp), intent(out)  :: lwork
         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)
+        !> Problem dimensions.
+        m = size(A, 1) ; n = size(A, 2) ; p = size(C, 1)
+        lwork = -1_ilp
+        !> Compute constrained lstsq solution.
+        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, err)
+        !> Optimal workspace size.
+        lwork = ceiling(real(work(1), kind=${rk}$), kind=ilp)
     end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space
 
     module subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq(A, b, C, d, x, storage, overwrite_matrices, err)
@@ -469,7 +482,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
             call handle_gglse_info(this, info, ma, na, mc, err0)
 
             !> Deallocate.
-            deallocate(work)
+            if (.not. present(storage)) deallocate(work)
         endif
 
         if (.not. overwrite_matrices_) then

From fc546449a515b23ba70bb953021bf01459eea09a Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 06:59:00 +0200
Subject: [PATCH 05/17] Added test for constrained lstsq.

---
 test/linalg/CMakeLists.txt                    |   2 +
 .../linalg/test_linalg_constrained_lstsq.fypp | 154 ++++++++++++++++++
 2 files changed, 156 insertions(+)
 create mode 100644 test/linalg/test_linalg_constrained_lstsq.fypp

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..f20ed6aac
--- /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, b, C, d, 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 (performances) -----
+        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

From 9582df2073c128298a9cf03d49878a2980a4073d Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 11:20:40 +0200
Subject: [PATCH 06/17] In-code documentation.

---
 src/stdlib_linalg.fypp               | 69 ++++++++++++++++++++++----
 src/stdlib_linalg_least_squares.fypp | 73 +++++++++++++++++++++++-----
 2 files changed, 119 insertions(+), 23 deletions(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 11bfb752d..7e9d3a5ea 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -582,12 +582,29 @@ module stdlib_linalg
   ! 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]())
+    !!
+    !!  ### 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)
-        !> Input matrices.
-        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
-        !> Right-hand side vectors.
-        ${rt}$, intent(inout), target :: b(:), d(:)
+        !> 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.
@@ -601,12 +618,30 @@ module stdlib_linalg
   ! 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]())
+    !!
+    !!  ### 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)
-        !> Input matrices.
-        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
-        !> Right-hand side vectors.
-        ${rt}$, intent(inout), target :: b(:), d(:)
+        !> 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.
@@ -620,10 +655,24 @@ module stdlib_linalg
   end interface
 
   interface constrained_lstsq_space
+    !!  version: experimental
+    !!
+    !!  Computes the size of the workspace required by the constrained least-squares solver.
+    !!  ([Specification]())
+    !!
+    !!  ### 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, b, C, d, lwork, err)
-        ${rt}$, intent(in), target :: A(:, :), C(:, :)
-        ${rt}$, intent(in), target :: b(:), d(:)
+        !> Least-squares cost.
+        ${rt}$, intent(in) :: A(:, :), b(:)
+        !> Equality constraints.
+        ${rt}$, intent(in) :: C(:, :), d(:)
         integer(ilp), intent(out) :: lwork
         type(linalg_state_type), optional, intent(out) :: err
     end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space
diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp
index ce7b6c697..b29ce2732 100644
--- a/src/stdlib_linalg_least_squares.fypp
+++ b/src/stdlib_linalg_least_squares.fypp
@@ -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(:,:)
@@ -376,7 +376,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
             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(m, n) =', [mc, nc])
+            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
@@ -389,10 +389,15 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
     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, b, C, d, lwork, err)
-        ${rt}$, intent(in), target :: A(:, :), C(:, :)
-        ${rt}$, intent(in), target :: b(:), d(:)
+        !> Least-squares cost.
+        ${rt}$, intent(in) :: A(:, :), b(:)
+        !> Equality constraints.
+        ${rt}$, intent(in) :: C(:, :), d(:)
+        !> 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
@@ -402,18 +407,40 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
         !> Problem dimensions.
         m = size(A, 1) ; n = size(A, 2) ; p = size(C, 1)
         lwork = -1_ilp
-        !> Compute constrained lstsq solution.
+        !> 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, err)
         !> Optimal workspace size.
         lwork = ceiling(real(work(1), kind=${rk}$), kind=ilp)
     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)
-        !> Input matrices.
-        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
-        !> Right-hand side vectors.
-        ${rt}$, intent(inout), target :: b(:), d(:)
+        !!  ### 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.
@@ -494,10 +521,30 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
     end subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq
 
     module function stdlib_linalg_${ri}$_constrained_lstsq(A, b, C, d, overwrite_matrices, err) result(x)
-        !> Input matrices.
-        ${rt}$, intent(inout), target :: A(:, :), C(:, :)
-        !> Right-hand side vectors.
-        ${rt}$, intent(inout), target :: b(:), d(:)
+        !!  ### 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.

From ddc0f0da3e36a3f7a72c6eba25e84acbe006fbd8 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 12:07:06 +0200
Subject: [PATCH 07/17] Bug fix.

---
 src/lapack/stdlib_linalg_lapack_aux.fypp | 1 +
 src/stdlib_linalg_least_squares.fypp     | 5 ++++-
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp
index 5670bfc95..7142302db 100644
--- a/src/lapack/stdlib_linalg_lapack_aux.fypp
+++ b/src/lapack/stdlib_linalg_lapack_aux.fypp
@@ -1620,6 +1620,7 @@ module stdlib_linalg_lapack_aux
                 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)
diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp
index b29ce2732..ed12bcaa9 100644
--- a/src/stdlib_linalg_least_squares.fypp
+++ b/src/stdlib_linalg_least_squares.fypp
@@ -404,14 +404,17 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
         ${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, err)
+        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.

From b306f914088c9063c75ffa78c0d81684e5f8d2fa Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 12:07:18 +0200
Subject: [PATCH 08/17] Added examples.

---
 example/linalg/CMakeLists.txt                 |  2 +
 example/linalg/example_constrained_lstsq1.f90 | 40 ++++++++++++++++
 example/linalg/example_constrained_lstsq2.f90 | 48 +++++++++++++++++++
 3 files changed, 90 insertions(+)
 create mode 100644 example/linalg/example_constrained_lstsq1.f90
 create mode 100644 example/linalg/example_constrained_lstsq2.f90

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..cc8d074cb
--- /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, b, C, d, 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

From ed45fb55489281453da678963432719b24d19c24 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 16:41:29 +0200
Subject: [PATCH 09/17] Specs for constrained_lstsq

---
 doc/specs/stdlib_linalg.md | 49 ++++++++++++++++++++++++++++++++++++++
 src/stdlib_linalg.fypp     |  2 +-
 2 files changed, 50 insertions(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 7ed224e21..1423dde89 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -860,6 +860,55 @@ 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!}
+```
+
 ## `det` - Computes the determinant of a square matrix
 
 ### Status
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 7e9d3a5ea..3964238c4 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -590,7 +590,7 @@ module stdlib_linalg
     !!  subject to      Cx = d
     !!
     !!  where A is of size `m x n` and C of size `p x n`.
-    !!  ([Specification]())
+    !!  ([Specification](../page/specs/stdlib_linalg.html#constrained-lstsq))
     !!
     !!  ### Description
     !!

From 816d29eac177aec26102cd6613a5b00bd9b626e8 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 16:47:46 +0200
Subject: [PATCH 10/17] Specs for solve_constrained_lstsq

---
 doc/specs/stdlib_linalg.md | 46 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 1423dde89..d8af548ff 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -909,6 +909,52 @@ Exceptions trigger an `error stop`.
 {!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!}
+```
+
 ## `det` - Computes the determinant of a square matrix
 
 ### Status

From e9d10f0d47fbb6fc54b06b3e46108f28c6c228d8 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 16:55:31 +0200
Subject: [PATCH 11/17] All specs + minor modifs.

---
 doc/specs/stdlib_linalg.md                    | 26 +++++++++++++++++++
 example/linalg/example_constrained_lstsq2.f90 |  2 +-
 src/stdlib_linalg.fypp                        |  6 ++---
 src/stdlib_linalg_least_squares.fypp          |  8 +++---
 .../linalg/test_linalg_constrained_lstsq.fypp |  2 +-
 5 files changed, 35 insertions(+), 9 deletions(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index d8af548ff..09405cb42 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -955,6 +955,32 @@ where \(A\) is an \( m \times n \) matrix (with \(m \geq n\)) and \(C\) a \( p
 {!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
diff --git a/example/linalg/example_constrained_lstsq2.f90 b/example/linalg/example_constrained_lstsq2.f90
index cc8d074cb..997181f72 100644
--- a/example/linalg/example_constrained_lstsq2.f90
+++ b/example/linalg/example_constrained_lstsq2.f90
@@ -31,7 +31,7 @@ program example_constrained_lstsq2
    d = [1.0_dp, 3.0_dp, -1.0_dp]
 
    !> Optimal workspace size.
-   call constrained_lstsq_space(A, b, C, d, lwork)
+   call constrained_lstsq_space(A, C, lwork)
    allocate (work(lwork))
 
    ! Compute the solution.
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 3964238c4..40a385114 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -668,11 +668,11 @@ module stdlib_linalg
     !!  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, b, C, d, lwork, err)
+    module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, C, lwork, err)
         !> Least-squares cost.
-        ${rt}$, intent(in) :: A(:, :), b(:)
+        ${rt}$, intent(in) :: A(:, :)
         !> Equality constraints.
-        ${rt}$, intent(in) :: C(:, :), d(:)
+        ${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
diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp
index ed12bcaa9..bee76e40f 100644
--- a/src/stdlib_linalg_least_squares.fypp
+++ b/src/stdlib_linalg_least_squares.fypp
@@ -390,11 +390,11 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
 
     #: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, b, C, d, lwork, err)
+    module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, C, lwork, err)
         !> Least-squares cost.
-        ${rt}$, intent(in) :: A(:, :), b(:)
+        ${rt}$, intent(in) :: A(:, :)
         !> Equality constraints.
-        ${rt}$, intent(in) :: C(:, :), d(:)
+        ${rt}$, intent(in) :: C(:, :)
         !> Size of the workspace array.
         integer(ilp), intent(out)  :: lwork
         !> [optional] State return flag.
@@ -492,7 +492,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
         endif
 
         !> Retrieve workspace size.
-        call stdlib_linalg_${ri}$_constrained_lstsq_space(A, b, C, d, lwork, err0)
+        call stdlib_linalg_${ri}$_constrained_lstsq_space(A, C, lwork, err0)
 
         if (err0%ok()) then
             !> Workspace.
diff --git a/test/linalg/test_linalg_constrained_lstsq.fypp b/test/linalg/test_linalg_constrained_lstsq.fypp
index f20ed6aac..3b513f3d1 100644
--- a/test/linalg/test_linalg_constrained_lstsq.fypp
+++ b/test/linalg/test_linalg_constrained_lstsq.fypp
@@ -79,7 +79,7 @@ module test_linalg_constrained_least_squares
         if (allocated(error)) return
 
         !----- Pre-allocated storage -----
-        call constrained_lstsq_space(A, b, C, d, lwork, err=state)
+        call constrained_lstsq_space(A, C, lwork, err=state)
         call check(error, state%ok(), state%print())
         if (allocated(error)) return
         allocate(work(lwork))

From f0add5a996145994f452ed4fb7fefdc3333fac96 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 23 Oct 2025 17:01:47 +0200
Subject: [PATCH 12/17] Added links to the correct specs section.

---
 src/stdlib_linalg.fypp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 40a385114..8e8cec64d 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -626,7 +626,7 @@ module stdlib_linalg
     !!  subject to      Cx = d
     !!
     !!  where A is of size `m x n` and C of size `p x n`.
-    !!  ([Specification]())
+    !!  ([Specification](../page/specs/stdlib_linalg.html#solve-constrained-lstsq))
     !!
     !!  ### Description
     !!
@@ -658,7 +658,7 @@ module stdlib_linalg
     !!  version: experimental
     !!
     !!  Computes the size of the workspace required by the constrained least-squares solver.
-    !!  ([Specification]())
+    !!  ([Specification](../page/specs/stdlib_linalg.html#constrained-lstsq-space))
     !!
     !!  ### Description
     !!

From 5d53fe87bb90a18e7448be4afccc66a03b051038 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 4 Nov 2025 15:10:06 +0100
Subject: [PATCH 13/17] Update doc/specs/stdlib_linalg.md

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 doc/specs/stdlib_linalg.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 09405cb42..935c5c67a 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -979,7 +979,7 @@ call [stdlib_linalg(module):constrained_lstsq_space(interface)]`(a,  c,  lwork [
 
 `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.
+`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
 

From c43044a9c3900fc3f2fd17a8c249846c84dc5287 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 4 Nov 2025 15:10:25 +0100
Subject: [PATCH 14/17] Update src/stdlib_linalg.fypp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 src/stdlib_linalg.fypp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 8e8cec64d..26c60fd64 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -635,7 +635,7 @@ module stdlib_linalg
     !!  `complex`. If a pre-allocated workspace is provided, no internal memory allocation takes
     !!  place.
     !!
-    !!  @note The solution is based on LAPACK's `*GLLSE` methods.
+    !!  @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.

From ef6dee68dbcda1c55d9d320cd113832bd23072f7 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 4 Nov 2025 15:10:57 +0100
Subject: [PATCH 15/17] Update test/linalg/test_linalg_constrained_lstsq.fypp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 test/linalg/test_linalg_constrained_lstsq.fypp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/linalg/test_linalg_constrained_lstsq.fypp b/test/linalg/test_linalg_constrained_lstsq.fypp
index 3b513f3d1..215663362 100644
--- a/test/linalg/test_linalg_constrained_lstsq.fypp
+++ b/test/linalg/test_linalg_constrained_lstsq.fypp
@@ -91,7 +91,7 @@ module test_linalg_constrained_least_squares
         call check(error, all(abs(x-x_true) < 1.0e-4_${rk}$), 'Solver converged')
         if (allocated(error)) return
 
-        !----- Overwrite matrices (performances) -----
+        !----- 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())

From 5add96138a1994df9a0409fdd85d9a839b3f7d85 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 4 Nov 2025 15:11:06 +0100
Subject: [PATCH 16/17] Update src/stdlib_linalg.fypp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 src/stdlib_linalg.fypp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 26c60fd64..fe0fa4dca 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -598,7 +598,7 @@ module stdlib_linalg
     !!  least-squares problem using a function. Supported data types include `real` and
     !!  `complex`. 
     !!
-    !!  @note The solution is based on LAPACK's `*GLLSE` methods.
+    !!  @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

From fbb198e0144f7ea4f71f1030c5e22435db028802 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 4 Nov 2025 15:12:15 +0100
Subject: [PATCH 17/17] Update doc/specs/stdlib_linalg.md

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 doc/specs/stdlib_linalg.md | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 935c5c67a..1a9b852a9 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -973,12 +973,7 @@ call [stdlib_linalg(module):constrained_lstsq_space(interface)]`(a,  c,  lwork [
 
 `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 `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