From 9e16d374d658c4b03320c9c80c7c53ed75ae29cd Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Sat, 24 Jan 2026 18:48:57 +0000 Subject: [PATCH 01/12] feat(linalg): add generalized least squares solver --- doc/specs/stdlib_linalg.md | 44 ++++++ example/linalg/CMakeLists.txt | 1 + example/linalg/example_generalized_lstsq.f90 | 28 ++++ src/lapack/stdlib_linalg_lapack_aux.fypp | 27 +++- src/linalg/stdlib_linalg.fypp | 39 +++++ src/linalg/stdlib_linalg_least_squares.fypp | 151 ++++++++++++++++++- test/linalg/test_linalg_lstsq.fypp | 80 +++++++++- 7 files changed, 363 insertions(+), 7 deletions(-) create mode 100644 example/linalg/example_generalized_lstsq.f90 diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 082efaefa..93e0a6681 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -976,6 +976,50 @@ call [stdlib_linalg(module):constrained_lstsq_space(interface)]`(a, c, lwork [ `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. +## `generalized_lstsq` - Computes the generalized least squares solution {#generalized-lstsq} + +### Status + +Experimental + +### Description + +This function computes the generalized least-squares (GLS) solution to a linear matrix equation \( A \cdot x = b \) with correlated errors described by a covariance matrix \( W \). + +The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) where \( W \) is a symmetric positive definite covariance matrix. This is useful when observations have correlated errors. The solver is based on LAPACK's `*GGGLM` backends. + +### Syntax + +`x = ` [[stdlib_linalg(module):generalized_lstsq(interface)]] `(w, a, b [, prefactored_w, overwrite_a, err])` + +### Arguments + +`w`: Shall be a rank-2 `real` or `complex` symmetric positive definite array containing the covariance matrix, or its lower triangular Cholesky factor if `prefactored_w=.true.`. It is an `intent(in)` argument. + +`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. The matrix must have full column rank. It is an `intent(inout)` argument. + +`b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument. + +`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain the lower triangular Cholesky factor of the covariance matrix rather than the covariance matrix itself. Default: `.false.`. This is an `intent(in)` argument. + +`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns an array value of the same kind as `a`, containing the generalized least-squares solution. + +Raises `LINALG_VALUE_ERROR` if the covariance matrix is not square or has incompatible dimensions. +Raises `LINALG_ERROR` if the covariance matrix is not positive definite or if the design matrix does not have full column rank. +Exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_generalized_lstsq.f90!} +``` + ## `det` - Computes the determinant of a square matrix ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1c6c7e42c..c4e4a6238 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -37,6 +37,7 @@ ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(constrained_lstsq1) ADD_EXAMPLE(constrained_lstsq2) +ADD_EXAMPLE(generalized_lstsq) ADD_EXAMPLE(norm) ADD_EXAMPLE(mnorm) ADD_EXAMPLE(get_norm) diff --git a/example/linalg/example_generalized_lstsq.f90 b/example/linalg/example_generalized_lstsq.f90 new file mode 100644 index 000000000..c8f691b81 --- /dev/null +++ b/example/linalg/example_generalized_lstsq.f90 @@ -0,0 +1,28 @@ +! Generalized least-squares solver with correlated errors +program example_generalized_lstsq + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: generalized_lstsq + implicit none + + real(dp) :: A(3,2), b(3), W(3,3) + real(dp), allocatable :: x(:) + + ! Design matrix: intercept + slope + A(:,1) = 1.0_dp + A(:,2) = [1.0_dp, 2.0_dp, 3.0_dp] + + ! Observations + b = [1.0_dp, 2.1_dp, 2.9_dp] + + ! Covariance matrix (correlated errors) + W(1,:) = [1.0_dp, 0.5_dp, 0.25_dp] + W(2,:) = [0.5_dp, 1.0_dp, 0.5_dp] + W(3,:) = [0.25_dp, 0.5_dp, 1.0_dp] + + ! Solve generalized least-squares + x = generalized_lstsq(W, A, b) + + print '("GLS fit: intercept = ",f8.4,", slope = ",f8.4)', x(1), x(2) + ! GLS fit: intercept = 0.0500, slope = 0.9500 + +end program example_generalized_lstsq diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 475424c91..1a6e15eff 100644 --- a/src/lapack/stdlib_linalg_lapack_aux.fypp +++ b/src/lapack/stdlib_linalg_lapack_aux.fypp @@ -52,9 +52,10 @@ 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 + public :: handle_gglse_info + public :: handle_ggglm_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. ! An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if SELCTG is true, i.e., abstract interface @@ -1654,4 +1655,24 @@ module stdlib_linalg_lapack_aux end select end subroutine handle_gglse_info + elemental subroutine handle_ggglm_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 (0) + ! Success. + err%state = LINALG_SUCCESS + case (:-1) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid argument at position', -info) + case (1:) + ! From LAPACK: the upper triangular factor R of A does not have full rank + err = linalg_state_type(this, LINALG_ERROR, & + 'Design matrix A does not have full column rank (rank-deficient)') + case default + err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'catastrophic error.') + end select + end subroutine handle_ggglm_info + end module stdlib_linalg_lapack_aux diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index 06169e071..420c8e261 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -40,6 +40,7 @@ module stdlib_linalg public :: lstsq_space public :: constrained_lstsq public :: constrained_lstsq_space + public :: generalized_lstsq public :: norm public :: mnorm public :: get_norm @@ -679,6 +680,44 @@ module stdlib_linalg #:endfor end interface + interface generalized_lstsq + !! version: experimental + !! + !! Computes the generalized least-squares solution to \( \min_x (Ax-b)^T W^{-1} (Ax-b) \) + !! ([Specification](../page/specs/stdlib_linalg.html#generalized-lstsq)) + !! + !!### Summary + !! Function interface for computing generalized least-squares via GGGLM. + !! + !!### Description + !! + !! This interface provides methods for computing generalized least-squares + !! with a symmetric positive definite covariance matrix. + !! Supported data types include `real` and `complex`. + !! + !!@note The solution is based on LAPACK's `*GGGLM` routine. + !!@note The covariance matrix W is always preserved (copied internally). + !! + #:for rk,rt,ri in RC_KINDS_TYPES + module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,err) result(x) + !> Covariance matrix W[m,m] (SPD) or its lower triangular Cholesky factor + ${rt}$, intent(in) :: w(:,:) + !> Input matrix a[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector b[m] + ${rt}$, intent(in) :: b(:) + !> [optional] Is W already Cholesky-factored? Default: .false. + logical(lk), optional, intent(in) :: prefactored_w + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + !> Result array x[n] + ${rt}$, allocatable :: x(:) + end function stdlib_linalg_${ri}$_generalized_lstsq + #:endfor + end interface generalized_lstsq + ! QR factorization of rank-2 array A interface qr !! version: experimental diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index bee76e40f..a02234c24 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -7,10 +7,10 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !! Least-squares solution to Ax=b use stdlib_linalg_constants - use stdlib_linalg_lapack, only: gelsd, gglse, stdlib_ilaenv - use stdlib_linalg_lapack_aux, only: handle_gelsd_info, handle_gglse_info + use stdlib_linalg_lapack, only: gelsd, gglse, ggglm, potrf, stdlib_ilaenv + use stdlib_linalg_lapack_aux, only: handle_gelsd_info, handle_gglse_info, handle_ggglm_info use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & - LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS implicit none character(*), parameter :: this = 'lstsq' @@ -564,4 +564,149 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares end function stdlib_linalg_${ri}$_constrained_lstsq #:endfor + !------------------------------------------------------------- + !----- Generalized Least-Squares Solver ----- + !------------------------------------------------------------- + + #:for rk,rt,ri in RC_KINDS_TYPES + ! Generalized least-squares: minimize (Ax-b)^T W^{-1} (Ax-b) where W is SPD + module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,err) result(x) + !> Covariance matrix W[m,m] (SPD) or its lower triangular Cholesky factor + ${rt}$, intent(in) :: w(:,:) + !> Input matrix a[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector b[m] + ${rt}$, intent(in) :: b(:) + !> [optional] Is W already Cholesky-factored? Default: .false. + logical(lk), optional, intent(in) :: prefactored_w + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + !> Result array x[n] + ${rt}$, allocatable :: x(:) + + ! Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: m, n, p, lda, ldb, lwork, info, i, j + logical(lk) :: copy_a, is_prefactored + ${rt}$, pointer :: amat(:,:) + ${rt}$, allocatable, target :: amat_alloc(:,:) + ${rt}$, allocatable :: lmat(:,:), d(:), y(:), work(:) + character(*), parameter :: this = 'generalized_lstsq' + + m = size(a, 1, kind=ilp) + n = size(a, 2, kind=ilp) + p = m ! For GLS, B is m×m + + ! Allocate result early (prevents segfault on error return) + allocate(x(n)) + + ! Validate matrix dimensions + if (m < 1 .or. n < 1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid matrix size a(m, n) =', [m, n]) + call linalg_error_handling(err0, err) + return + end if + + ! Validate sizes + if (size(w, 1, kind=ilp) /= m .or. size(w, 2, kind=ilp) /= m) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, & + 'Covariance matrix must be square m×m:', [size(w, 1, kind=ilp), size(w, 2, kind=ilp)]) + call linalg_error_handling(err0, err) + return + end if + + if (size(b, kind=ilp) /= m) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, & + 'Right-hand side size must match rows of A:', size(b, kind=ilp), '/=', m) + call linalg_error_handling(err0, err) + return + end if + + if (m < n) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, & + 'GGGLM requires m >= n (overdetermined or square):', m, '<', n) + call linalg_error_handling(err0, err) + return + end if + + ! Process options + is_prefactored = .false._lk + if (present(prefactored_w)) is_prefactored = prefactored_w + + if (present(overwrite_a)) then + copy_a = .not.overwrite_a + else + copy_a = .true._lk + end if + + ! Handle A matrix + if (copy_a) then + allocate(amat_alloc(m, n), source=a) + amat => amat_alloc + else + amat => a + end if + + ! Handle covariance/Cholesky factor + ! ALWAYS copy W because GGGLM modifies it (protects user's data) + allocate(lmat(m, m), source=w) + + if (.not. is_prefactored) then + ! Compute Cholesky factorization: W = L * L^T + call potrf('L', m, lmat, m, info) + if (info /= 0) then + if (info > 0) then + err0 = linalg_state_type(this, LINALG_ERROR, & + 'Covariance matrix is not positive definite at position', info) + else + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, & + 'Invalid argument to POTRF at position', -info) + end if + ! Cleanup before early return + if (copy_a) deallocate(amat_alloc) + deallocate(lmat) + call linalg_error_handling(err0, err) + return + end if + end if + + ! Prepare for GGGLM + allocate(d(m), source=b) + allocate(y(p)) + + lda = m + ldb = m + + ! Zero out upper triangle of lmat (GGGLM reads full matrix, + ! but potrf only sets lower triangle) + do j = 1, m + do i = 1, j - 1 + lmat(i, j) = 0.0_${rk}$ + end do + end do + + ! Workspace query + allocate(work(1)) + call ggglm(m, n, p, amat, lda, lmat, ldb, d, x, y, work, -1_ilp, info) + lwork = ceiling(real(work(1), kind=${rk}$), kind=ilp) + deallocate(work) + allocate(work(lwork)) + + ! Solve GLS via GGGLM + call ggglm(m, n, p, amat, lda, lmat, ldb, d, x, y, work, lwork, info) + + ! Handle errors + call handle_ggglm_info(this, info, m, n, p, err0) + + ! Cleanup + if (copy_a) deallocate(amat_alloc) + deallocate(lmat, d, y, work) + + call linalg_error_handling(err0, err) + + end function stdlib_linalg_${ri}$_generalized_lstsq + #:endfor + end submodule stdlib_linalg_least_squares diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 7ec41a2e8..e759972c6 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -4,7 +4,7 @@ module test_linalg_least_squares use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg_constants - use stdlib_linalg, only: lstsq,solve_lstsq + use stdlib_linalg, only: lstsq, solve_lstsq, generalized_lstsq use stdlib_linalg_state, only: linalg_state_type implicit none (type,external) @@ -26,6 +26,8 @@ module test_linalg_least_squares #:for rk,rt,ri in REAL_KINDS_TYPES call add_test(tests,new_unittest("least_squares_${ri}$",test_lstsq_one_${ri}$)) call add_test(tests,new_unittest("least_squares_randm_${ri}$",test_lstsq_random_${ri}$)) + call add_test(tests,new_unittest("generalized_lstsq_${ri}$",test_generalized_lstsq_${ri}$)) + call add_test(tests,new_unittest("generalized_lstsq_identity_${ri}$",test_generalized_lstsq_identity_${ri}$)) #:endfor end subroutine test_least_squares @@ -139,6 +141,82 @@ module test_linalg_least_squares end subroutine test_issue_823 + #:for rk,rt,ri in REAL_KINDS_TYPES + !> Test basic generalized least-squares with correlated errors + subroutine test_generalized_lstsq_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp), parameter :: m = 3, n = 2 + real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: A(m,n), b(m), W(m,m) + ${rt}$, allocatable :: x(:) + + ! Design matrix + A(:,1) = 1.0_${rk}$ + A(:,2) = [1.0_${rk}$, 2.0_${rk}$, 3.0_${rk}$] + + ! Observations + b = [1.0_${rk}$, 2.1_${rk}$, 2.9_${rk}$] + + ! SPD covariance matrix (correlated errors) + W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] + W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] + W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + + x = generalized_lstsq(W, A, b, err=state) + + call check(error, state%ok(), 'generalized_lstsq failed: '//state%print()) + if (allocated(error)) return + + call check(error, size(x)==n, 'generalized_lstsq: wrong solution size') + if (allocated(error)) return + + ! Solution should be finite + call check(error, all(abs(x) < huge(0.0_${rk}$)), 'generalized_lstsq: solution not finite') + if (allocated(error)) return + + end subroutine test_generalized_lstsq_${ri}$ + + !> Test GLS with identity covariance = OLS + subroutine test_generalized_lstsq_identity_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp), parameter :: m = 4, n = 2 + integer(ilp) :: i + real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: A(m,n), b(m), W(m,m) + ${rt}$, allocatable :: x_gls(:), x_ols(:) + + ! Design matrix + A(:,1) = 1.0_${rk}$ + A(:,2) = [1.0_${rk}$, 2.0_${rk}$, 3.0_${rk}$, 4.0_${rk}$] + b = [2.0_${rk}$, 4.0_${rk}$, 5.0_${rk}$, 4.0_${rk}$] + + ! Identity covariance + W = 0.0_${rk}$ + do i = 1, m + W(i,i) = 1.0_${rk}$ + end do + + x_gls = generalized_lstsq(W, A, b, err=state) + call check(error, state%ok(), 'generalized_lstsq with I failed') + if (allocated(error)) return + + x_ols = lstsq(A, b, err=state) + call check(error, state%ok(), 'lstsq failed') + if (allocated(error)) return + + ! GLS with identity should equal OLS + call check(error, all(abs(x_gls - x_ols) < tol), & + 'generalized_lstsq with identity should equal lstsq') + if (allocated(error)) return + + end subroutine test_generalized_lstsq_identity_${ri}$ + + #:endfor + ! gcc-15 bugfix utility subroutine add_test(tests,new_test) type(unittest_type), allocatable, intent(inout) :: tests(:) From 82ab37151b4162a4a0db5ba5079c62efdd8602cb Mon Sep 17 00:00:00 2001 From: Amrinder Singh Date: Mon, 26 Jan 2026 02:54:26 +0530 Subject: [PATCH 02/12] Update src/linalg/stdlib_linalg_least_squares.fypp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/linalg/stdlib_linalg_least_squares.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index a02234c24..ff7263e94 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -569,9 +569,9 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !------------------------------------------------------------- #:for rk,rt,ri in RC_KINDS_TYPES - ! Generalized least-squares: minimize (Ax-b)^T W^{-1} (Ax-b) where W is SPD + ! Generalized least-squares: minimize (Ax-b)^T W^{-1} (Ax-b) where W is symmetric/Hermitian positive definite module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,err) result(x) - !> Covariance matrix W[m,m] (SPD) or its lower triangular Cholesky factor + !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its lower triangular Cholesky factor ${rt}$, intent(in) :: w(:,:) !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) From fbb6c69107681b7f943e8f50a6b45ab714f2837f Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Sun, 25 Jan 2026 21:34:52 +0000 Subject: [PATCH 03/12] docs: clarify SPD (real) / HPD (complex) for GLS, add zero parameter --- doc/specs/stdlib_linalg.md | 4 ++-- src/linalg/stdlib_linalg.fypp | 4 ++-- src/linalg/stdlib_linalg_least_squares.fypp | 5 +++-- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 93e0a6681..b7741bedc 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -986,7 +986,7 @@ Experimental This function computes the generalized least-squares (GLS) solution to a linear matrix equation \( A \cdot x = b \) with correlated errors described by a covariance matrix \( W \). -The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) where \( W \) is a symmetric positive definite covariance matrix. This is useful when observations have correlated errors. The solver is based on LAPACK's `*GGGLM` backends. +The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) where \( W \) is a symmetric (real) or Hermitian (complex) positive definite covariance matrix. This is useful when observations have correlated errors. The solver is based on LAPACK's `*GGGLM` backends. ### Syntax @@ -994,7 +994,7 @@ The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) w ### Arguments -`w`: Shall be a rank-2 `real` or `complex` symmetric positive definite array containing the covariance matrix, or its lower triangular Cholesky factor if `prefactored_w=.true.`. It is an `intent(in)` argument. +`w`: Shall be a rank-2 `real` or `complex` symmetric (real) or Hermitian (complex) positive definite array containing the covariance matrix, or its lower triangular Cholesky factor if `prefactored_w=.true.`. It is an `intent(in)` argument. `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. The matrix must have full column rank. It is an `intent(inout)` argument. diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index 420c8e261..d0a2e3c5d 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -692,7 +692,7 @@ module stdlib_linalg !!### Description !! !! This interface provides methods for computing generalized least-squares - !! with a symmetric positive definite covariance matrix. + !! with a symmetric (real) or Hermitian (complex) positive definite covariance matrix. !! Supported data types include `real` and `complex`. !! !!@note The solution is based on LAPACK's `*GGGLM` routine. @@ -700,7 +700,7 @@ module stdlib_linalg !! #:for rk,rt,ri in RC_KINDS_TYPES module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,err) result(x) - !> Covariance matrix W[m,m] (SPD) or its lower triangular Cholesky factor + !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its lower triangular Cholesky factor ${rt}$, intent(in) :: w(:,:) !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index ff7263e94..ab0e852b3 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -593,6 +593,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ${rt}$, pointer :: amat(:,:) ${rt}$, allocatable, target :: amat_alloc(:,:) ${rt}$, allocatable :: lmat(:,:), d(:), y(:), work(:) + ${rt}$, parameter :: zero = 0.0_${rk}$ character(*), parameter :: this = 'generalized_lstsq' m = size(a, 1, kind=ilp) @@ -654,7 +655,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares allocate(lmat(m, m), source=w) if (.not. is_prefactored) then - ! Compute Cholesky factorization: W = L * L^T + ! Compute Cholesky factorization: W = L * L^T (real) or W = L * L^H (complex) call potrf('L', m, lmat, m, info) if (info /= 0) then if (info > 0) then @@ -683,7 +684,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! but potrf only sets lower triangle) do j = 1, m do i = 1, j - 1 - lmat(i, j) = 0.0_${rk}$ + lmat(i, j) = zero end do end do From 703ae17eebd930320e8df630c1a4ae86f8ef95e2 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Sun, 25 Jan 2026 21:48:35 +0000 Subject: [PATCH 04/12] test: remove unused tol variable in generalized_lstsq test --- test/linalg/test_linalg_lstsq.fypp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index e759972c6..78c5a0d9b 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -148,7 +148,6 @@ module test_linalg_least_squares type(linalg_state_type) :: state integer(ilp), parameter :: m = 3, n = 2 - real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) ${rt}$ :: A(m,n), b(m), W(m,m) ${rt}$, allocatable :: x(:) From 4764aa99e653dcba955b02adb0d38a0eb263fc77 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Mon, 26 Jan 2026 14:30:11 +0000 Subject: [PATCH 05/12] refactor(generalized_lstsq): use optval, cholesky subroutine, and do concurrent --- src/linalg/stdlib_linalg_least_squares.fypp | 35 ++++++--------------- 1 file changed, 10 insertions(+), 25 deletions(-) diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index ab0e852b3..eb4d90a0a 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -7,7 +7,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !! Least-squares solution to Ax=b use stdlib_linalg_constants - use stdlib_linalg_lapack, only: gelsd, gglse, ggglm, potrf, stdlib_ilaenv + use stdlib_linalg_lapack, only: gelsd, gglse, ggglm, stdlib_ilaenv use stdlib_linalg_lapack_aux, only: handle_gelsd_info, handle_gglse_info, handle_ggglm_info use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS @@ -633,14 +633,9 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares end if ! Process options - is_prefactored = .false._lk - if (present(prefactored_w)) is_prefactored = prefactored_w + is_prefactored = optval(prefactored_w, .false._lk) - if (present(overwrite_a)) then - copy_a = .not.overwrite_a - else - copy_a = .true._lk - end if + copy_a = .not. optval(overwrite_a, .false._lk) ! Handle A matrix if (copy_a) then @@ -656,21 +651,19 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares if (.not. is_prefactored) then ! Compute Cholesky factorization: W = L * L^T (real) or W = L * L^H (complex) - call potrf('L', m, lmat, m, info) - if (info /= 0) then - if (info > 0) then - err0 = linalg_state_type(this, LINALG_ERROR, & - 'Covariance matrix is not positive definite at position', info) - else - err0 = linalg_state_type(this, LINALG_VALUE_ERROR, & - 'Invalid argument to POTRF at position', -info) - end if + call cholesky(lmat, lower=.true._lk, other_zeroed=.true._lk, err=err0) + if (err0%error()) then ! Cleanup before early return if (copy_a) deallocate(amat_alloc) deallocate(lmat) call linalg_error_handling(err0, err) return end if + else + ! User provided pre-factored L: zero out upper triangle for GGGLM + do concurrent(i=1:m, j=1:m, i < j) + lmat(i, j) = zero + end do end if ! Prepare for GGGLM @@ -680,14 +673,6 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares lda = m ldb = m - ! Zero out upper triangle of lmat (GGGLM reads full matrix, - ! but potrf only sets lower triangle) - do j = 1, m - do i = 1, j - 1 - lmat(i, j) = zero - end do - end do - ! Workspace query allocate(work(1)) call ggglm(m, n, p, amat, lda, lmat, ldb, d, x, y, work, -1_ilp, info) From 296727a40e94ea6ce612eabfd26be0f2389a7016 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Mon, 26 Jan 2026 15:08:52 +0000 Subject: [PATCH 06/12] fix(handle_ggglm_info): expand select case with specific error messages --- src/lapack/stdlib_linalg_lapack_aux.fypp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 1a6e15eff..7f054d0d2 100644 --- a/src/lapack/stdlib_linalg_lapack_aux.fypp +++ b/src/lapack/stdlib_linalg_lapack_aux.fypp @@ -1664,8 +1664,18 @@ module stdlib_linalg_lapack_aux case (0) ! Success. err%state = LINALG_SUCCESS - case (:-1) - err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid argument at position', -info) + case (-1) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of rows for A and B, m=', m) + case (-2) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of columns for A, n=', n) + case (-3) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of columns for B, p=', p) + case (-5) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda < m=', m) + case (-7) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for B, ldb < m=', m) + case (-12) + err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'Insufficient workspace size.') case (1:) ! From LAPACK: the upper triangular factor R of A does not have full rank err = linalg_state_type(this, LINALG_ERROR, & From 5a70184e075599f4883515b0d017f1d90022ec7f Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Fri, 30 Jan 2026 21:04:11 +0000 Subject: [PATCH 07/12] Add overwrite_w, SVD-based sqrt test, improve error messages --- doc/specs/stdlib_linalg.md | 10 ++- src/lapack/stdlib_linalg_lapack_aux.fypp | 12 +-- src/linalg/stdlib_linalg.fypp | 11 +-- src/linalg/stdlib_linalg_least_squares.fypp | 50 ++++++----- test/linalg/test_linalg_lstsq.fypp | 95 ++++++++++++++++++++- 5 files changed, 140 insertions(+), 38 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index b7741bedc..294115d3d 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -990,19 +990,21 @@ The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) w ### Syntax -`x = ` [[stdlib_linalg(module):generalized_lstsq(interface)]] `(w, a, b [, prefactored_w, overwrite_a, err])` +`x = ` [[stdlib_linalg(module):generalized_lstsq(interface)]] `(w, a, b [, prefactored_w, overwrite_a, overwrite_w, err])` ### Arguments -`w`: Shall be a rank-2 `real` or `complex` symmetric (real) or Hermitian (complex) positive definite array containing the covariance matrix, or its lower triangular Cholesky factor if `prefactored_w=.true.`. It is an `intent(in)` argument. +`w`: Shall be a rank-2 `real` or `complex` symmetric (real) or Hermitian (complex) positive definite array containing the covariance matrix, or a valid matrix square root (e.g., Cholesky factor, SVD-based) if `prefactored_w=.true.`. It is an `intent(inout)` argument. `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. The matrix must have full column rank. It is an `intent(inout)` argument. `b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument. -`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain the lower triangular Cholesky factor of the covariance matrix rather than the covariance matrix itself. Default: `.false.`. This is an `intent(in)` argument. +`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). Default: `.false.`. This is an `intent(in)` argument. -`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. +`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. Default: `.false.`. This is an `intent(in)` argument. + +`overwrite_w` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `w` will be used as temporary storage and overwritten. This avoids internal data allocation. Default: `.false.`. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 7f054d0d2..1af3ab328 100644 --- a/src/lapack/stdlib_linalg_lapack_aux.fypp +++ b/src/lapack/stdlib_linalg_lapack_aux.fypp @@ -1655,9 +1655,9 @@ module stdlib_linalg_lapack_aux end select end subroutine handle_gglse_info - elemental subroutine handle_ggglm_info(this, info, m, n, p, err) + elemental subroutine handle_ggglm_info(this, info, m, n, p, lda, ldb, err) character(len=*), intent(in) :: this - integer(ilp), intent(in) :: info, m, n, p + integer(ilp), intent(in) :: info, m, n, p, lda, ldb type(linalg_state_type), intent(out) :: err ! Process output. select case (info) @@ -1671,9 +1671,11 @@ module stdlib_linalg_lapack_aux case (-3) err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of columns for B, p=', p) case (-5) - err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda < m=', m) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda=', lda, & + ' is < m=', m) case (-7) - err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for B, ldb < m=', m) + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for B, ldb=', ldb, & + ' is < m=', m) case (-12) err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'Insufficient workspace size.') case (1:) @@ -1683,6 +1685,6 @@ module stdlib_linalg_lapack_aux case default err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'catastrophic error.') end select - end subroutine handle_ggglm_info + end subroutine handle_ggglm_info end module stdlib_linalg_lapack_aux diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index d0a2e3c5d..b6ffdf1be 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -696,20 +696,21 @@ module stdlib_linalg !! Supported data types include `real` and `complex`. !! !!@note The solution is based on LAPACK's `*GGGLM` routine. - !!@note The covariance matrix W is always preserved (copied internally). !! #:for rk,rt,ri in RC_KINDS_TYPES - module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,err) result(x) - !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its lower triangular Cholesky factor - ${rt}$, intent(in) :: w(:,:) + module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,overwrite_w,err) result(x) + !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + ${rt}$, intent(inout), target :: w(:,:) !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector b[m] ${rt}$, intent(in) :: b(:) - !> [optional] Is W already Cholesky-factored? Default: .false. + !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. logical(lk), optional, intent(in) :: prefactored_w !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Can W data be overwritten and destroyed? Default: .false. + logical(lk), optional, intent(in) :: overwrite_w !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err !> Result array x[n] diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index eb4d90a0a..40fadc3dc 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -570,17 +570,19 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:for rk,rt,ri in RC_KINDS_TYPES ! Generalized least-squares: minimize (Ax-b)^T W^{-1} (Ax-b) where W is symmetric/Hermitian positive definite - module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,err) result(x) - !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its lower triangular Cholesky factor - ${rt}$, intent(in) :: w(:,:) + module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,overwrite_w,err) result(x) + !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + ${rt}$, intent(inout), target :: w(:,:) !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector b[m] ${rt}$, intent(in) :: b(:) - !> [optional] Is W already Cholesky-factored? Default: .false. + !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. logical(lk), optional, intent(in) :: prefactored_w !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Can W data be overwritten and destroyed? Default: .false. + logical(lk), optional, intent(in) :: overwrite_w !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err !> Result array x[n] @@ -588,12 +590,11 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Local variables type(linalg_state_type) :: err0 - integer(ilp) :: m, n, p, lda, ldb, lwork, info, i, j - logical(lk) :: copy_a, is_prefactored - ${rt}$, pointer :: amat(:,:) - ${rt}$, allocatable, target :: amat_alloc(:,:) - ${rt}$, allocatable :: lmat(:,:), d(:), y(:), work(:) - ${rt}$, parameter :: zero = 0.0_${rk}$ + integer(ilp) :: m, n, p, lda, ldb, lwork, info + logical(lk) :: copy_a, copy_w, is_prefactored + ${rt}$, pointer :: amat(:,:), lmat(:,:) + ${rt}$, allocatable, target :: amat_alloc(:,:), lmat_alloc(:,:) + ${rt}$, allocatable :: d(:), y(:), work(:) character(*), parameter :: this = 'generalized_lstsq' m = size(a, 1, kind=ilp) @@ -613,7 +614,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Validate sizes if (size(w, 1, kind=ilp) /= m .or. size(w, 2, kind=ilp) /= m) then err0 = linalg_state_type(this, LINALG_VALUE_ERROR, & - 'Covariance matrix must be square m×m:', [size(w, 1, kind=ilp), size(w, 2, kind=ilp)]) + 'Covariance matrix must be square m×m:', shape(w, kind=ilp)) call linalg_error_handling(err0, err) return end if @@ -634,8 +635,8 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Process options is_prefactored = optval(prefactored_w, .false._lk) - copy_a = .not. optval(overwrite_a, .false._lk) + copy_w = .not. optval(overwrite_w, .false._lk) ! Handle A matrix if (copy_a) then @@ -645,9 +646,13 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares amat => a end if - ! Handle covariance/Cholesky factor - ! ALWAYS copy W because GGGLM modifies it (protects user's data) - allocate(lmat(m, m), source=w) + ! Handle covariance/matrix square root + if (copy_w) then + allocate(lmat_alloc(m, m), source=w) + lmat => lmat_alloc + else + lmat => w + end if if (.not. is_prefactored) then ! Compute Cholesky factorization: W = L * L^T (real) or W = L * L^H (complex) @@ -655,16 +660,14 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares if (err0%error()) then ! Cleanup before early return if (copy_a) deallocate(amat_alloc) - deallocate(lmat) + if (copy_w) deallocate(lmat_alloc) call linalg_error_handling(err0, err) return end if - else - ! User provided pre-factored L: zero out upper triangle for GGGLM - do concurrent(i=1:m, j=1:m, i < j) - lmat(i, j) = zero - end do end if + ! If prefactored_w=.true., user provides a valid matrix square root B where W = B*B^T. + ! This can be a Cholesky factor OR any other valid square root (e.g., SVD-based). + ! We do not modify the user's input. ! Prepare for GGGLM allocate(d(m), source=b) @@ -684,11 +687,12 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares call ggglm(m, n, p, amat, lda, lmat, ldb, d, x, y, work, lwork, info) ! Handle errors - call handle_ggglm_info(this, info, m, n, p, err0) + call handle_ggglm_info(this, info, m, n, p, lda, ldb, err0) ! Cleanup if (copy_a) deallocate(amat_alloc) - deallocate(lmat, d, y, work) + if (copy_w) deallocate(lmat_alloc) + deallocate(d, y, work) call linalg_error_handling(err0, err) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 78c5a0d9b..78f968cfd 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -4,7 +4,7 @@ module test_linalg_least_squares use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg_constants - use stdlib_linalg, only: lstsq, solve_lstsq, generalized_lstsq + use stdlib_linalg, only: lstsq, solve_lstsq, generalized_lstsq, eigh, diag use stdlib_linalg_state, only: linalg_state_type implicit none (type,external) @@ -28,6 +28,8 @@ module test_linalg_least_squares call add_test(tests,new_unittest("least_squares_randm_${ri}$",test_lstsq_random_${ri}$)) call add_test(tests,new_unittest("generalized_lstsq_${ri}$",test_generalized_lstsq_${ri}$)) call add_test(tests,new_unittest("generalized_lstsq_identity_${ri}$",test_generalized_lstsq_identity_${ri}$)) + call add_test(tests,new_unittest("generalized_lstsq_svd_sqrt_${ri}$",test_generalized_lstsq_svd_sqrt_${ri}$)) + call add_test(tests,new_unittest("generalized_lstsq_overwrite_${ri}$",test_generalized_lstsq_overwrite_${ri}$)) #:endfor end subroutine test_least_squares @@ -214,6 +216,97 @@ module test_linalg_least_squares end subroutine test_generalized_lstsq_identity_${ri}$ + !> Test GLS with eigendecomposition-based matrix square root (not Cholesky) + !> This verifies that prefactored_w works with any valid matrix square root + subroutine test_generalized_lstsq_svd_sqrt_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp), parameter :: m = 3, n = 2 + integer(ilp) :: i + real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: A(m,n), b(m), W(m,m), sqrt_W(m,m), U(m,m) + real(${rk}$) :: lambda(m) + ${rt}$, allocatable :: x_chol(:), x_sqrt(:) + + ! Design matrix + A(:,1) = 1.0_${rk}$ + A(:,2) = [1.0_${rk}$, 2.0_${rk}$, 3.0_${rk}$] + + ! Observations + b = [1.0_${rk}$, 2.1_${rk}$, 2.9_${rk}$] + + ! SPD covariance matrix + W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] + W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] + W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + + ! Solve using Cholesky (internal factorization) + x_chol = generalized_lstsq(W, A, b, err=state) + call check(error, state%ok(), 'generalized_lstsq (Cholesky) failed: '//state%print()) + if (allocated(error)) return + + ! Compute eigendecomposition: W = U * diag(lambda) * U^T + call eigh(W, lambda, vectors=U, err=state) + call check(error, state%ok(), 'eigh failed: '//state%print()) + if (allocated(error)) return + + ! Compute matrix square root: sqrt(W) = U * diag(sqrt(lambda)) * U^T + ! This is a DENSE matrix (not triangular like Cholesky factor) + sqrt_W = matmul(U, matmul(diag(sqrt(lambda)), transpose(U))) + + ! Solve using eigendecomposition-based square root (prefactored) + x_sqrt = generalized_lstsq(sqrt_W, A, b, prefactored_w=.true., err=state) + call check(error, state%ok(), 'generalized_lstsq (sqrt) failed: '//state%print()) + if (allocated(error)) return + + ! Both methods should give the same solution + call check(error, all(abs(x_chol - x_sqrt) < tol), & + 'generalized_lstsq: Cholesky and eigendecomposition-based sqrt should match') + if (allocated(error)) return + + end subroutine test_generalized_lstsq_svd_sqrt_${ri}$ + + !> Test GLS with overwrite_w=.true. (avoid internal allocation) + subroutine test_generalized_lstsq_overwrite_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp), parameter :: m = 3, n = 2 + real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: A(m,n), b(m), W(m,m), W_copy(m,m) + ${rt}$, allocatable :: x_copy(:), x_overwrite(:) + + ! Design matrix + A(:,1) = 1.0_${rk}$ + A(:,2) = [1.0_${rk}$, 2.0_${rk}$, 3.0_${rk}$] + + ! Observations + b = [1.0_${rk}$, 2.1_${rk}$, 2.9_${rk}$] + + ! SPD covariance matrix + W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] + W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] + W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + W_copy = W + + ! Solve with copy (default) + x_copy = generalized_lstsq(W_copy, A, b, err=state) + call check(error, state%ok(), 'generalized_lstsq (copy) failed: '//state%print()) + if (allocated(error)) return + + ! Solve with overwrite (no internal allocation for W) + x_overwrite = generalized_lstsq(W, A, b, overwrite_w=.true., err=state) + call check(error, state%ok(), 'generalized_lstsq (overwrite) failed: '//state%print()) + if (allocated(error)) return + + ! Both should give the same solution + call check(error, all(abs(x_copy - x_overwrite) < tol), & + 'generalized_lstsq: copy and overwrite solutions should match') + if (allocated(error)) return + + end subroutine test_generalized_lstsq_overwrite_${ri}$ + #:endfor ! gcc-15 bugfix utility From 0b5f9ad41c4f0c48957ab38824bd38d9fd19b777 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Sat, 7 Feb 2026 12:05:18 +0000 Subject: [PATCH 08/12] Use eye/hermitian in tests, document Cholesky zeroing, fix error propagation --- doc/specs/stdlib_linalg.md | 2 +- src/linalg/stdlib_linalg_least_squares.fypp | 2 +- test/linalg/test_linalg_lstsq.fypp | 14 +++++--------- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 294115d3d..9062c6b80 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1000,7 +1000,7 @@ The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) w `b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument. -`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). Default: `.false.`. This is an `intent(in)` argument. +`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). If a Cholesky factor is used, the upper triangular part must be zeroed out. Default: `.false.`. This is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. Default: `.false.`. This is an `intent(in)` argument. diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index 40fadc3dc..1b5db6ca7 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -661,7 +661,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Cleanup before early return if (copy_a) deallocate(amat_alloc) if (copy_w) deallocate(lmat_alloc) - call linalg_error_handling(err0, err) + call linalg_error_handling(err0, err, where_at=this) return end if end if diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 78f968cfd..614f91e72 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -4,7 +4,7 @@ module test_linalg_least_squares use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg_constants - use stdlib_linalg, only: lstsq, solve_lstsq, generalized_lstsq, eigh, diag + use stdlib_linalg, only: lstsq, solve_lstsq, generalized_lstsq, eigh, diag, eye, hermitian use stdlib_linalg_state, only: linalg_state_type implicit none (type,external) @@ -185,7 +185,6 @@ module test_linalg_least_squares type(linalg_state_type) :: state integer(ilp), parameter :: m = 4, n = 2 - integer(ilp) :: i real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) ${rt}$ :: A(m,n), b(m), W(m,m) ${rt}$, allocatable :: x_gls(:), x_ols(:) @@ -196,10 +195,7 @@ module test_linalg_least_squares b = [2.0_${rk}$, 4.0_${rk}$, 5.0_${rk}$, 4.0_${rk}$] ! Identity covariance - W = 0.0_${rk}$ - do i = 1, m - W(i,i) = 1.0_${rk}$ - end do + W = eye(m, mold=0.0_${rk}$) x_gls = generalized_lstsq(W, A, b, err=state) call check(error, state%ok(), 'generalized_lstsq with I failed') @@ -223,7 +219,6 @@ module test_linalg_least_squares type(linalg_state_type) :: state integer(ilp), parameter :: m = 3, n = 2 - integer(ilp) :: i real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) ${rt}$ :: A(m,n), b(m), W(m,m), sqrt_W(m,m), U(m,m) real(${rk}$) :: lambda(m) @@ -251,9 +246,10 @@ module test_linalg_least_squares call check(error, state%ok(), 'eigh failed: '//state%print()) if (allocated(error)) return - ! Compute matrix square root: sqrt(W) = U * diag(sqrt(lambda)) * U^T + ! Compute matrix square root: sqrt(W) = U * diag(sqrt(lambda)) * U^H ! This is a DENSE matrix (not triangular like Cholesky factor) - sqrt_W = matmul(U, matmul(diag(sqrt(lambda)), transpose(U))) + ! Using hermitian() instead of transpose() for complex-safety + sqrt_W = matmul(U, matmul(diag(sqrt(lambda)), hermitian(U))) ! Solve using eigendecomposition-based square root (prefactored) x_sqrt = generalized_lstsq(sqrt_W, A, b, prefactored_w=.true., err=state) From 07417f8d93e339c6fe5352351cff9d93c5cf49eb Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Tue, 10 Feb 2026 09:38:08 +0000 Subject: [PATCH 09/12] removed where_at error handling --- src/linalg/stdlib_linalg_least_squares.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index 1b5db6ca7..40fadc3dc 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -661,7 +661,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Cleanup before early return if (copy_a) deallocate(amat_alloc) if (copy_w) deallocate(lmat_alloc) - call linalg_error_handling(err0, err, where_at=this) + call linalg_error_handling(err0, err) return end if end if From 8c82fabf6c5f5cae6eea1d4e61da06747fb7846a Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Sat, 14 Feb 2026 18:02:44 +0000 Subject: [PATCH 10/12] add solve_generalized_lstsq subroutine, fix docs and allocate style --- doc/specs/stdlib_linalg.md | 50 ++++++++++++++++++- example/linalg/CMakeLists.txt | 1 + .../example_solve_generalized_lstsq.f90 | 27 ++++++++++ src/linalg/stdlib_linalg.fypp | 42 +++++++++++++++- src/linalg/stdlib_linalg_least_squares.fypp | 40 ++++++++++++--- test/linalg/test_linalg_lstsq.fypp | 42 +++++++++++++++- 6 files changed, 192 insertions(+), 10 deletions(-) create mode 100644 example/linalg/example_solve_generalized_lstsq.f90 diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 9062c6b80..ef191d7b5 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1010,7 +1010,7 @@ The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) w ### Return value -Returns an array value of the same kind as `a`, containing the generalized least-squares solution. +Returns a rank-1 array of the same kind as `a`, containing the generalized least-squares solution. Raises `LINALG_VALUE_ERROR` if the covariance matrix is not square or has incompatible dimensions. Raises `LINALG_ERROR` if the covariance matrix is not positive definite or if the design matrix does not have full column rank. @@ -1022,6 +1022,54 @@ Exceptions trigger an `error stop`. {!example/linalg/example_generalized_lstsq.f90!} ``` +## `solve_generalized_lstsq` - Solves the generalized least squares problem (subroutine) {#solve-generalized-lstsq} + +### Status + +Experimental + +### Description + +This subroutine computes the generalized least-squares (GLS) solution to a linear matrix equation \( A \cdot x = b \) with correlated errors described by a covariance matrix \( W \). + +The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) where \( W \) is a symmetric (real) or Hermitian (complex) positive definite covariance matrix. Unlike `generalized_lstsq`, this subroutine interface allows the user to provide a pre-allocated solution vector `x`. The solver is based on LAPACK's `*GGGLM` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):solve_generalized_lstsq(interface)]] `(w, a, b, x [, prefactored_w, overwrite_a, overwrite_w, err])` + +### Arguments + +`w`: Shall be a rank-2 `real` or `complex` symmetric (real) or Hermitian (complex) positive definite array containing the covariance matrix, or a valid matrix square root (e.g., Cholesky factor, SVD-based) if `prefactored_w=.true.`. It is an `intent(inout)` argument. + +`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. The matrix must have full column rank. It is an `intent(inout)` argument. + +`b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument. + +`x`: Shall be a rank-1 array of the same kind as `a`, containing the solution vector. It is an `intent(out)` argument. + +`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). Default: `.false.`. This is an `intent(in)` argument. + +`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. Default: `.false.`. This is an `intent(in)` argument. + +`overwrite_w` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `w` will be used as temporary storage and overwritten. This avoids internal data allocation. Default: `.false.`. This is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns the generalized least-squares solution in the pre-allocated rank-1 array `x`. + +Raises `LINALG_VALUE_ERROR` if the covariance matrix is not square or has incompatible dimensions. +Raises `LINALG_ERROR` if the covariance matrix is not positive definite or if the design matrix does not have full column rank. +Exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_solve_generalized_lstsq.f90!} +``` + ## `det` - Computes the determinant of a square matrix ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index c4e4a6238..338c62d36 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -38,6 +38,7 @@ ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(constrained_lstsq1) ADD_EXAMPLE(constrained_lstsq2) ADD_EXAMPLE(generalized_lstsq) +ADD_EXAMPLE(solve_generalized_lstsq) ADD_EXAMPLE(norm) ADD_EXAMPLE(mnorm) ADD_EXAMPLE(get_norm) diff --git a/example/linalg/example_solve_generalized_lstsq.f90 b/example/linalg/example_solve_generalized_lstsq.f90 new file mode 100644 index 000000000..cf90a6974 --- /dev/null +++ b/example/linalg/example_solve_generalized_lstsq.f90 @@ -0,0 +1,27 @@ +! Generalized least-squares solver (subroutine interface) +program example_solve_generalized_lstsq + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: solve_generalized_lstsq + implicit none + + real(dp) :: A(3,2), b(3), W(3,3), x(2) + + ! Design matrix: intercept + slope + A(:,1) = 1.0_dp + A(:,2) = [1.0_dp, 2.0_dp, 3.0_dp] + + ! Observations + b = [1.0_dp, 2.1_dp, 2.9_dp] + + ! Covariance matrix (correlated errors) + W(1,:) = [1.0_dp, 0.5_dp, 0.25_dp] + W(2,:) = [0.5_dp, 1.0_dp, 0.5_dp] + W(3,:) = [0.25_dp, 0.5_dp, 1.0_dp] + + ! Solve generalized least-squares + call solve_generalized_lstsq(W, A, b, x) + + print '("GLS fit: intercept = ",f8.4,", slope = ",f8.4)', x(1), x(2) + ! GLS fit: intercept = 0.0500, slope = 0.9500 + +end program example_solve_generalized_lstsq diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index b6ffdf1be..904a6d50f 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -41,6 +41,7 @@ module stdlib_linalg public :: constrained_lstsq public :: constrained_lstsq_space public :: generalized_lstsq + public :: solve_generalized_lstsq public :: norm public :: mnorm public :: get_norm @@ -714,11 +715,50 @@ module stdlib_linalg !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err !> Result array x[n] - ${rt}$, allocatable :: x(:) + ${rt}$, allocatable, target :: x(:) end function stdlib_linalg_${ri}$_generalized_lstsq #:endfor end interface generalized_lstsq + interface solve_generalized_lstsq + !! version: experimental + !! + !! Computes the generalized least-squares solution to \( \min_x (Ax-b)^T W^{-1} (Ax-b) \) + !! ([Specification](../page/specs/stdlib_linalg.html#solve-generalized-lstsq)) + !! + !!### Summary + !! Subroutine interface for computing generalized least-squares via GGGLM. + !! + !!### Description + !! + !! This interface provides methods for computing generalized least-squares + !! with a symmetric (real) or Hermitian (complex) positive definite covariance matrix. + !! Supported data types include `real` and `complex`. + !! + !!@note The solution is based on LAPACK's `*GGGLM` routine. + !! + #:for rk,rt,ri in RC_KINDS_TYPES + module subroutine stdlib_linalg_${ri}$_solve_generalized_lstsq(w,a,b,x,prefactored_w,overwrite_a,overwrite_w,err) + !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + ${rt}$, intent(inout), target :: w(:,:) + !> Input matrix a[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector b[m] + ${rt}$, intent(in) :: b(:) + !> Solution vector x[n] + ${rt}$, intent(out) :: x(:) + !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. + logical(lk), optional, intent(in) :: prefactored_w + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Can W data be overwritten and destroyed? Default: .false. + logical(lk), optional, intent(in) :: overwrite_w + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_${ri}$_solve_generalized_lstsq + #:endfor + end interface solve_generalized_lstsq + ! QR factorization of rank-2 array A interface qr !! version: experimental diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index 40fadc3dc..aaafd9372 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -569,7 +569,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !------------------------------------------------------------- #:for rk,rt,ri in RC_KINDS_TYPES - ! Generalized least-squares: minimize (Ax-b)^T W^{-1} (Ax-b) where W is symmetric/Hermitian positive definite + ! Generalized least-squares function: thin wrapper that allocates x and delegates to subroutine module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,overwrite_w,err) result(x) !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root ${rt}$, intent(inout), target :: w(:,:) @@ -586,7 +586,36 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err !> Result array x[n] - ${rt}$, allocatable :: x(:) + ${rt}$, allocatable, target :: x(:) + + integer(ilp) :: n + + n = size(a, 2, kind=ilp) + allocate(x(n)) + call stdlib_linalg_${ri}$_solve_generalized_lstsq(w, a, b, x, & + prefactored_w=prefactored_w, overwrite_a=overwrite_a, overwrite_w=overwrite_w, err=err) + end function stdlib_linalg_${ri}$_generalized_lstsq + #:endfor + + #:for rk,rt,ri in RC_KINDS_TYPES + ! Generalized least-squares subroutine: minimize (Ax-b)^T W^{-1} (Ax-b) + module subroutine stdlib_linalg_${ri}$_solve_generalized_lstsq(w,a,b,x,prefactored_w,overwrite_a,overwrite_w,err) + !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + ${rt}$, intent(inout), target :: w(:,:) + !> Input matrix a[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector b[m] + ${rt}$, intent(in) :: b(:) + !> Solution vector x[n] + ${rt}$, intent(out) :: x(:) + !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. + logical(lk), optional, intent(in) :: prefactored_w + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Can W data be overwritten and destroyed? Default: .false. + logical(lk), optional, intent(in) :: overwrite_w + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err ! Local variables type(linalg_state_type) :: err0 @@ -601,9 +630,6 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares n = size(a, 2, kind=ilp) p = m ! For GLS, B is m×m - ! Allocate result early (prevents segfault on error return) - allocate(x(n)) - ! Validate matrix dimensions if (m < 1 .or. n < 1) then err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid matrix size a(m, n) =', [m, n]) @@ -670,7 +696,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! We do not modify the user's input. ! Prepare for GGGLM - allocate(d(m), source=b) + allocate(d, source=b) allocate(y(p)) lda = m @@ -696,7 +722,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares call linalg_error_handling(err0, err) - end function stdlib_linalg_${ri}$_generalized_lstsq + end subroutine stdlib_linalg_${ri}$_solve_generalized_lstsq #:endfor end submodule stdlib_linalg_least_squares diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 614f91e72..65055b6f6 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -4,7 +4,8 @@ module test_linalg_least_squares use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg_constants - use stdlib_linalg, only: lstsq, solve_lstsq, generalized_lstsq, eigh, diag, eye, hermitian + use stdlib_linalg, only: lstsq, solve_lstsq, generalized_lstsq, solve_generalized_lstsq, & + eigh, diag, eye, hermitian use stdlib_linalg_state, only: linalg_state_type implicit none (type,external) @@ -30,6 +31,7 @@ module test_linalg_least_squares call add_test(tests,new_unittest("generalized_lstsq_identity_${ri}$",test_generalized_lstsq_identity_${ri}$)) call add_test(tests,new_unittest("generalized_lstsq_svd_sqrt_${ri}$",test_generalized_lstsq_svd_sqrt_${ri}$)) call add_test(tests,new_unittest("generalized_lstsq_overwrite_${ri}$",test_generalized_lstsq_overwrite_${ri}$)) + call add_test(tests,new_unittest("solve_generalized_lstsq_${ri}$",test_solve_generalized_lstsq_${ri}$)) #:endfor end subroutine test_least_squares @@ -303,6 +305,44 @@ module test_linalg_least_squares end subroutine test_generalized_lstsq_overwrite_${ri}$ + !> Test solve_generalized_lstsq matches generalized_lstsq + subroutine test_solve_generalized_lstsq_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp), parameter :: m = 4, n = 2 + real(${rk}$), parameter :: tol = 1000*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: A(m,n), b(m), W(m,m), x_sub(n) + ${rt}$, allocatable :: x_fun(:) + + ! Design matrix + A(:,1) = 1.0_${rk}$ + A(:,2) = [1.0_${rk}$, 2.0_${rk}$, 3.0_${rk}$, 4.0_${rk}$] + b = [2.0_${rk}$, 4.0_${rk}$, 5.0_${rk}$, 4.0_${rk}$] + + ! SPD covariance matrix + W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$, 0.25_${rk}$] + W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] + W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] + W(4,:) = [0.25_${rk}$, 0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + + ! Function form + x_fun = generalized_lstsq(W, A, b, err=state) + call check(error, state%ok(), 'generalized_lstsq (function) failed: '//state%print()) + if (allocated(error)) return + + ! Subroutine form + call solve_generalized_lstsq(W, A, b, x_sub, err=state) + call check(error, state%ok(), 'solve_generalized_lstsq failed: '//state%print()) + if (allocated(error)) return + + ! Both should give the same solution + call check(error, all(abs(x_fun - x_sub) < tol), & + 'solve_generalized_lstsq should match generalized_lstsq') + if (allocated(error)) return + + end subroutine test_solve_generalized_lstsq_${ri}$ + #:endfor ! gcc-15 bugfix utility From 34e840d0752589989109de17a92a34cdfc1ea746 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Mon, 9 Mar 2026 17:06:17 +0000 Subject: [PATCH 11/12] fix: address review comments for generalized_lstsq --- doc/specs/stdlib_linalg.md | 6 ++--- example/linalg/CMakeLists.txt | 1 - example/linalg/example_generalized_lstsq.f90 | 23 +++++++++------- .../example_solve_generalized_lstsq.f90 | 27 ------------------- src/linalg/stdlib_linalg.fypp | 20 +++++++------- src/linalg/stdlib_linalg_least_squares.fypp | 25 ++++++++++------- test/linalg/test_linalg_lstsq.fypp | 26 +++++++++--------- 7 files changed, 55 insertions(+), 73 deletions(-) delete mode 100644 example/linalg/example_solve_generalized_lstsq.f90 diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ef191d7b5..04f358fe4 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1000,7 +1000,7 @@ The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) w `b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument. -`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). If a Cholesky factor is used, the upper triangular part must be zeroed out. Default: `.false.`. This is an `intent(in)` argument. +`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). If a Cholesky factor is used, the upper triangular part must be zeroed out. It is the user's responsibility to ensure that the prefactored matrix is valid. Default: `.false.`. This is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. Default: `.false.`. This is an `intent(in)` argument. @@ -1048,7 +1048,7 @@ The solver minimizes the Mahalanobis distance \( (Ax - b)^T W^{-1} (Ax - b) \) w `x`: Shall be a rank-1 array of the same kind as `a`, containing the solution vector. It is an `intent(out)` argument. -`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). Default: `.false.`. This is an `intent(in)` argument. +`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). It is the user's responsibility to ensure that the prefactored matrix is valid. Default: `.false.`. This is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. Default: `.false.`. This is an `intent(in)` argument. @@ -1067,7 +1067,7 @@ Exceptions trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_solve_generalized_lstsq.f90!} +{!example/linalg/example_generalized_lstsq.f90!} ``` ## `det` - Computes the determinant of a square matrix diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 338c62d36..c4e4a6238 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -38,7 +38,6 @@ ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(constrained_lstsq1) ADD_EXAMPLE(constrained_lstsq2) ADD_EXAMPLE(generalized_lstsq) -ADD_EXAMPLE(solve_generalized_lstsq) ADD_EXAMPLE(norm) ADD_EXAMPLE(mnorm) ADD_EXAMPLE(get_norm) diff --git a/example/linalg/example_generalized_lstsq.f90 b/example/linalg/example_generalized_lstsq.f90 index c8f691b81..638340daf 100644 --- a/example/linalg/example_generalized_lstsq.f90 +++ b/example/linalg/example_generalized_lstsq.f90 @@ -1,11 +1,12 @@ ! Generalized least-squares solver with correlated errors program example_generalized_lstsq use stdlib_linalg_constants, only: dp - use stdlib_linalg, only: generalized_lstsq + use stdlib_linalg, only: generalized_lstsq, solve_generalized_lstsq implicit none - real(dp) :: A(3,2), b(3), W(3,3) - real(dp), allocatable :: x(:) + integer, parameter :: m = 3 + real(dp) :: A(m,2), b(m), W(m,m), x(2) + real(dp), allocatable :: x_fun(:) ! Design matrix: intercept + slope A(:,1) = 1.0_dp @@ -15,14 +16,16 @@ program example_generalized_lstsq b = [1.0_dp, 2.1_dp, 2.9_dp] ! Covariance matrix (correlated errors) - W(1,:) = [1.0_dp, 0.5_dp, 0.25_dp] - W(2,:) = [0.5_dp, 1.0_dp, 0.5_dp] - W(3,:) = [0.25_dp, 0.5_dp, 1.0_dp] + W = reshape([1.0_dp, 0.5_dp, 0.25_dp, & + 0.5_dp, 1.0_dp, 0.5_dp, & + 0.25_dp, 0.5_dp, 1.0_dp], [m, m]) - ! Solve generalized least-squares - x = generalized_lstsq(W, A, b) + ! Function interface: allocates solution + x_fun = generalized_lstsq(W, A, b) + print '("GLS (function): intercept = ",f8.4,", slope = ",f8.4)', x_fun(1), x_fun(2) - print '("GLS fit: intercept = ",f8.4,", slope = ",f8.4)', x(1), x(2) - ! GLS fit: intercept = 0.0500, slope = 0.9500 + ! Subroutine interface: user-provided solution vector + call solve_generalized_lstsq(W, A, b, x) + print '("GLS (subroutine): intercept = ",f8.4,", slope = ",f8.4)', x(1), x(2) end program example_generalized_lstsq diff --git a/example/linalg/example_solve_generalized_lstsq.f90 b/example/linalg/example_solve_generalized_lstsq.f90 deleted file mode 100644 index cf90a6974..000000000 --- a/example/linalg/example_solve_generalized_lstsq.f90 +++ /dev/null @@ -1,27 +0,0 @@ -! Generalized least-squares solver (subroutine interface) -program example_solve_generalized_lstsq - use stdlib_linalg_constants, only: dp - use stdlib_linalg, only: solve_generalized_lstsq - implicit none - - real(dp) :: A(3,2), b(3), W(3,3), x(2) - - ! Design matrix: intercept + slope - A(:,1) = 1.0_dp - A(:,2) = [1.0_dp, 2.0_dp, 3.0_dp] - - ! Observations - b = [1.0_dp, 2.1_dp, 2.9_dp] - - ! Covariance matrix (correlated errors) - W(1,:) = [1.0_dp, 0.5_dp, 0.25_dp] - W(2,:) = [0.5_dp, 1.0_dp, 0.5_dp] - W(3,:) = [0.25_dp, 0.5_dp, 1.0_dp] - - ! Solve generalized least-squares - call solve_generalized_lstsq(W, A, b, x) - - print '("GLS fit: intercept = ",f8.4,", slope = ",f8.4)', x(1), x(2) - ! GLS fit: intercept = 0.0500, slope = 0.9500 - -end program example_solve_generalized_lstsq diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index 904a6d50f..84612500a 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -41,7 +41,6 @@ module stdlib_linalg public :: constrained_lstsq public :: constrained_lstsq_space public :: generalized_lstsq - public :: solve_generalized_lstsq public :: norm public :: mnorm public :: get_norm @@ -49,6 +48,7 @@ module stdlib_linalg public :: solve_lu public :: solve_lstsq public :: solve_constrained_lstsq + public :: solve_generalized_lstsq public :: trace public :: svd public :: svdvals @@ -700,11 +700,11 @@ module stdlib_linalg !! #:for rk,rt,ri in RC_KINDS_TYPES module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,overwrite_w,err) result(x) - !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + !> Covariance matrix W(m,m) (symmetric/Hermitian positive definite) or its matrix square root ${rt}$, intent(inout), target :: w(:,:) - !> Input matrix a[m,n] + !> Input matrix a(m,n) ${rt}$, intent(inout), target :: a(:,:) - !> Right hand side vector b[m] + !> Right hand side vector b(m) ${rt}$, intent(in) :: b(:) !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. logical(lk), optional, intent(in) :: prefactored_w @@ -714,8 +714,8 @@ module stdlib_linalg logical(lk), optional, intent(in) :: overwrite_w !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err - !> Result array x[n] - ${rt}$, allocatable, target :: x(:) + !> Result array x(n) + ${rt}$, allocatable :: x(:) end function stdlib_linalg_${ri}$_generalized_lstsq #:endfor end interface generalized_lstsq @@ -739,13 +739,13 @@ module stdlib_linalg !! #:for rk,rt,ri in RC_KINDS_TYPES module subroutine stdlib_linalg_${ri}$_solve_generalized_lstsq(w,a,b,x,prefactored_w,overwrite_a,overwrite_w,err) - !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + !> Covariance matrix W(m,m) (symmetric/Hermitian positive definite) or its matrix square root ${rt}$, intent(inout), target :: w(:,:) - !> Input matrix a[m,n] + !> Input matrix a(m,n) ${rt}$, intent(inout), target :: a(:,:) - !> Right hand side vector b[m] + !> Right hand side vector b(m) ${rt}$, intent(in) :: b(:) - !> Solution vector x[n] + !> Solution vector x(n) ${rt}$, intent(out) :: x(:) !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. logical(lk), optional, intent(in) :: prefactored_w diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index aaafd9372..eaf10d55d 100644 --- a/src/linalg/stdlib_linalg_least_squares.fypp +++ b/src/linalg/stdlib_linalg_least_squares.fypp @@ -571,11 +571,11 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:for rk,rt,ri in RC_KINDS_TYPES ! Generalized least-squares function: thin wrapper that allocates x and delegates to subroutine module function stdlib_linalg_${ri}$_generalized_lstsq(w,a,b,prefactored_w,overwrite_a,overwrite_w,err) result(x) - !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + !> Covariance matrix W(m,m) (symmetric/Hermitian positive definite) or its matrix square root ${rt}$, intent(inout), target :: w(:,:) - !> Input matrix a[m,n] + !> Input matrix a(m,n) ${rt}$, intent(inout), target :: a(:,:) - !> Right hand side vector b[m] + !> Right hand side vector b(m) ${rt}$, intent(in) :: b(:) !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. logical(lk), optional, intent(in) :: prefactored_w @@ -585,8 +585,8 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares logical(lk), optional, intent(in) :: overwrite_w !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err - !> Result array x[n] - ${rt}$, allocatable, target :: x(:) + !> Result array x(n) + ${rt}$, allocatable :: x(:) integer(ilp) :: n @@ -600,13 +600,13 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:for rk,rt,ri in RC_KINDS_TYPES ! Generalized least-squares subroutine: minimize (Ax-b)^T W^{-1} (Ax-b) module subroutine stdlib_linalg_${ri}$_solve_generalized_lstsq(w,a,b,x,prefactored_w,overwrite_a,overwrite_w,err) - !> Covariance matrix W[m,m] (symmetric/Hermitian positive definite) or its matrix square root + !> Covariance matrix W(m,m) (symmetric/Hermitian positive definite) or its matrix square root ${rt}$, intent(inout), target :: w(:,:) - !> Input matrix a[m,n] + !> Input matrix a(m,n) ${rt}$, intent(inout), target :: a(:,:) - !> Right hand side vector b[m] + !> Right hand side vector b(m) ${rt}$, intent(in) :: b(:) - !> Solution vector x[n] + !> Solution vector x(n) ${rt}$, intent(out) :: x(:) !> [optional] Is W already a matrix square root (e.g., Cholesky factor)? Default: .false. logical(lk), optional, intent(in) :: prefactored_w @@ -659,6 +659,13 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares return end if + if (size(x, kind=ilp) /= n) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, & + 'Solution vector size must match columns of A:', size(x, kind=ilp), '/=', n) + call linalg_error_handling(err0, err) + return + end if + ! Process options is_prefactored = optval(prefactored_w, .false._lk) copy_a = .not. optval(overwrite_a, .false._lk) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 65055b6f6..a0330493e 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -163,9 +163,9 @@ module test_linalg_least_squares b = [1.0_${rk}$, 2.1_${rk}$, 2.9_${rk}$] ! SPD covariance matrix (correlated errors) - W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] - W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] - W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + W = reshape([2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$, & + 1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$, & + 0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$], [m, m]) x = generalized_lstsq(W, A, b, err=state) @@ -234,9 +234,9 @@ module test_linalg_least_squares b = [1.0_${rk}$, 2.1_${rk}$, 2.9_${rk}$] ! SPD covariance matrix - W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] - W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] - W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + W = reshape([2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$, & + 1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$, & + 0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$], [m, m]) ! Solve using Cholesky (internal factorization) x_chol = generalized_lstsq(W, A, b, err=state) @@ -283,9 +283,9 @@ module test_linalg_least_squares b = [1.0_${rk}$, 2.1_${rk}$, 2.9_${rk}$] ! SPD covariance matrix - W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] - W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] - W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + W = reshape([2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$, & + 1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$, & + 0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$], [m, m]) W_copy = W ! Solve with copy (default) @@ -321,10 +321,10 @@ module test_linalg_least_squares b = [2.0_${rk}$, 4.0_${rk}$, 5.0_${rk}$, 4.0_${rk}$] ! SPD covariance matrix - W(1,:) = [2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$, 0.25_${rk}$] - W(2,:) = [1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$] - W(3,:) = [0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$] - W(4,:) = [0.25_${rk}$, 0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$] + W = reshape([2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$, 0.25_${rk}$, & + 1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$, 0.5_${rk}$, & + 0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$, 1.0_${rk}$, & + 0.25_${rk}$, 0.5_${rk}$, 1.0_${rk}$, 2.0_${rk}$], [m, m]) ! Function form x_fun = generalized_lstsq(W, A, b, err=state) From 3640093406e3c20e03dfd88d265197e489e57f9a Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Mar 2026 08:43:23 +0100 Subject: [PATCH 12/12] Update src/linalg/stdlib_linalg.fypp Co-authored-by: Jeremie Vandenplas --- src/linalg/stdlib_linalg.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index 84612500a..bd8955f12 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -712,7 +712,7 @@ module stdlib_linalg logical(lk), optional, intent(in) :: overwrite_a !> [optional] Can W data be overwritten and destroyed? Default: .false. logical(lk), optional, intent(in) :: overwrite_w - !> [optional] state return flag. On error if not requested, the code will stop + !> [optional] state return flag. On error if not requested, the code will trigger an error stop type(linalg_state_type), optional, intent(out) :: err !> Result array x(n) ${rt}$, allocatable :: x(:)