feat(linalg): add generalized least squares solver

doc/specs/stdlib_linalg.md

@@ -976,6 +976,100 @@ 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, 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.
+
+`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.
+
+`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 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.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!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). 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.
+
+`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_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,31 @@
+! Generalized least-squares solver with correlated errors
+program example_generalized_lstsq
+  use stdlib_linalg_constants, only: dp
+  use stdlib_linalg, only: generalized_lstsq, solve_generalized_lstsq
+  implicit none
+
+  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
+  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 = 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])
+
+  ! 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)
+
+  ! 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

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,36 @@ module stdlib_linalg_lapack_aux
         end select
     end subroutine handle_gglse_info
 
+   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, lda, ldb
+        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=', lda, &
+                      ' is < m=', m)
+            case (-7)
+                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:)
+                ! 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,13 +40,15 @@ module stdlib_linalg
   public :: lstsq_space
   public :: constrained_lstsq
   public :: constrained_lstsq_space
+  public :: generalized_lstsq
   public :: norm
   public :: mnorm
   public :: get_norm
   public :: solve
   public :: solve_lu  
   public :: solve_lstsq
   public :: solve_constrained_lstsq
+  public :: solve_generalized_lstsq
   public :: trace
   public :: svd
   public :: svdvals
@@ -679,6 +681,84 @@ 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.
+    !!
+    #: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
+        ${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 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 trigger an error 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
+
+  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