feat(linalg): add generalized least squares solver

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 (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
+
+`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 (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.
+
+`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

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)

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

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,34 @@ 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 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, &
+                      '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

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 (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.
+    !!@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(:,:)
+        !> 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 

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, 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,135 @@ 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 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(:,:)
+        !> 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(:)
+        ${rt}$, parameter :: zero = 0.0_${rk}$
+        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 = optval(prefactored_w, .false._lk)
+
+        copy_a = .not. optval(overwrite_a, .false._lk)
+
+        ! 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 (real) or W = L * L^H (complex)
+            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
+        allocate(d(m), source=b)
+        allocate(y(p))
+
+        lda = m
+        ldb = m
+
+        ! 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