feat(linalg): Add weighted and generalized least squares solvers

doc/specs/stdlib_linalg.md

@@ -976,6 +976,96 @@ 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.
 
+## `weighted_lstsq` - Computes the weighted least squares solution to a linear matrix equation {#weighted-lstsq}
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes the weighted least-squares solution to a linear matrix equation \( A \cdot x = b \) where each observation has a different weight.
+
+The solver minimizes the weighted 2-norm \( \| D(Ax - b) \|_2^2 \) where \( D = \mathrm{diag}(\sqrt{w}) \) is a diagonal matrix formed from the square roots of the weight vector. This is equivalent to transforming the problem to ordinary least-squares through row scaling. The solver is based on LAPACK's `*GELSD` backends.
+
+### Syntax
+
+`x = ` [[stdlib_linalg(module):weighted_lstsq(interface)]] `(w, a, b [, cond, overwrite_a, rank, err])`
+
+### Arguments
+
+`w`: Shall be a rank-1 `real` or `complex` array containing the weight vector. All weights must be positive. It is an `intent(in)` argument.
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. 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.
+
+`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `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.
+
+`rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `a`. This is an `intent(out)` 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 weighted least-squares solution.
+
+Raises `LINALG_VALUE_ERROR` if any weight is non-positive or if the matrix and weight/right-hand-side have incompatible sizes.
+Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_weighted_lstsq.f90!}
+```
+
+## `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

example/linalg/CMakeLists.txt

@@ -37,6 +37,8 @@ ADD_EXAMPLE(lstsq1)
 ADD_EXAMPLE(lstsq2)
 ADD_EXAMPLE(constrained_lstsq1)
 ADD_EXAMPLE(constrained_lstsq2)
+ADD_EXAMPLE(weighted_lstsq)
+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

example/linalg/example_weighted_lstsq.f90

@@ -0,0 +1,26 @@
+! Weighted least-squares solver
+program example_weighted_lstsq
+  use stdlib_linalg_constants, only: dp
+  use stdlib_linalg, only: weighted_lstsq
+  implicit none
+
+  real(dp) :: A(4,2), b(4), w(4)
+  real(dp), allocatable :: x(:)
+
+  ! Design matrix: intercept + slope
+  A(:,1) = 1.0_dp
+  A(:,2) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp]
+
+  ! Observations (one outlier at position 3)
+  b = [2.1_dp, 4.0_dp, 10.0_dp, 7.9_dp]
+
+  ! Weights: downweight the outlier
+  w = [1.0_dp, 1.0_dp, 0.1_dp, 1.0_dp]
+
+  ! Solve weighted least-squares
+  x = weighted_lstsq(w, A, b)
+
+  print '("Weighted fit: intercept = ",f8.4,", slope = ",f8.4)', x(1), x(2)
+  ! Weighted fit: intercept =   0.0667, slope =   1.9556
+
+end program example_weighted_lstsq

src/lapack/stdlib_linalg_lapack_aux.fypp

@@ -53,6 +53,7 @@ module stdlib_linalg_lapack_aux
      public :: handle_ggev_info
      public :: handle_heev_info
      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. 
@@ -1654,4 +1655,25 @@ module stdlib_linalg_lapack_aux
         end select
     end subroutine handle_gglse_info
 
+    !> Process GGGLM output flags
+    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

src/linalg/stdlib_linalg.fypp

@@ -40,6 +40,8 @@ module stdlib_linalg
   public :: lstsq_space
   public :: constrained_lstsq
   public :: constrained_lstsq_space
+  public :: weighted_lstsq
+  public :: generalized_lstsq
   public :: norm
   public :: mnorm
   public :: get_norm
@@ -679,6 +681,87 @@ module stdlib_linalg
     #:endfor
   end interface
 
+  ! Weighted least-squares: minimize ||D(Ax - b)||^2 where D = diag(sqrt(w))
+  interface weighted_lstsq
+    !! version: experimental
+    !!
+    !! Computes the weighted least-squares solution to \( \min_x \|D(Ax - b)\|_2^2 \)
+    !! ([Specification](../page/specs/stdlib_linalg.html#weighted-lstsq))
+    !!
+    !!### Summary
+    !! Function interface for computing weighted least-squares via row scaling.
+    !!
+    !!### Description
+    !!
+    !! This interface provides methods for computing weighted least-squares by
+    !! transforming to ordinary least-squares through row scaling.
+    !! Supported data types include `real` and `complex`.
+    !!
+    !!@note The solution is based on LAPACK's `*GELSD` after applying diagonal weights.
+    !!@warning Avoid extreme weight ratios (e.g., max(w)/min(w) > 1e6) as this may 
+    !!         cause loss of precision in the SVD-based solver.
+    !!
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    module function stdlib_linalg_${ri}$_weighted_lstsq(w,a,b,cond,overwrite_a,rank,err) result(x)
+        !> Weight vector (must be positive, always real)
+        real(${rk}$), intent(in) :: w(:)
+        !> Input matrix a[m,n]
+        ${rt}$, intent(inout), target :: a(:,:)
+        !> Right hand side vector b[m]
+        ${rt}$, intent(in) :: b(:)
+        !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
+        real(${rk}$), optional, intent(in) :: cond
+        !> [optional] Can A data be overwritten and destroyed?
+        logical(lk), optional, intent(in) :: overwrite_a
+        !> [optional] Return rank of A
+        integer(ilp), optional, intent(out) :: rank
+        !> [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}$_weighted_lstsq
+    #:endfor
+  end interface weighted_lstsq
+
+  ! Generalized least-squares: minimize (Ax-b)^T W^{-1} (Ax-b) where W is SPD
+  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