diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 082efaefa..10f42e716 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/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
diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt
index c6c41f78f..c5c0b4dd7 100644
--- a/example/linalg/CMakeLists.txt
+++ b/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)
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/example/linalg/example_weighted_lstsq.f90 b/example/linalg/example_weighted_lstsq.f90
new file mode 100644
index 000000000..849c1bde5
--- /dev/null
+++ b/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
diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp
index 475424c91..887aebaf5 100644
--- a/src/lapack/stdlib_linalg_lapack_aux.fypp
+++ b/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
diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp
index 06169e071..1104a84ab 100644
--- a/src/linalg/stdlib_linalg.fypp
+++ b/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 
diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp
index bee76e40f..a8db2b05f 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,259 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
     end function stdlib_linalg_${ri}$_constrained_lstsq
     #:endfor
 
+    !-------------------------------------------------------------
+    !-----     Weighted Least-Squares Solver                 -----
+    !-------------------------------------------------------------
+
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    ! Weighted least-squares: minimize ||D(Ax - b)||^2 where D = diag(sqrt(w))
+    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(:)
+
+        ! Local variables
+        type(linalg_state_type) :: err0
+        integer(ilp) :: m, n, j
+        logical(lk) :: copy_a
+        ${rt}$, pointer :: amat(:,:)
+        ${rt}$, allocatable, target :: amat_alloc(:,:)
+        ${rt}$, allocatable :: b_scaled(:)
+        real(${rk}$), allocatable :: sqrt_w(:)
+        character(*), parameter :: this = 'weighted_lstsq'
+
+        m = size(a, 1, kind=ilp)
+        n = size(a, 2, kind=ilp)
+
+        ! Allocate result (even on error, to prevent segfault on 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)
+            if (present(rank)) rank = 0
+            return
+        end if
+
+        ! Validate inputs
+        if (size(w, kind=ilp) /= m) then
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, &
+                   'Weight vector size must match number of rows:', size(w, kind=ilp), '/=', m)
+            call linalg_error_handling(err0, err)
+            if (present(rank)) rank = 0
+            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)
+            if (present(rank)) rank = 0
+            return
+        end if
+
+        if (any(w <= 0.0_${rk}$)) then
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Weights must be positive')
+            call linalg_error_handling(err0, err)
+            if (present(rank)) rank = 0
+            return
+        end if
+
+        ! Can A be overwritten? By default, do not overwrite
+        if (present(overwrite_a)) then
+            copy_a = .not.overwrite_a
+        else
+            copy_a = .true._lk
+        end if
+
+        ! Compute sqrt of weights
+        allocate(sqrt_w(m))
+        sqrt_w = sqrt(w)
+
+        ! Handle A matrix: either copy or use original
+        if (copy_a) then
+            allocate(amat_alloc(m, n))
+            amat => amat_alloc
+        else
+            amat => a
+        end if
+
+        ! Scale A column-wise (cache-friendly: column-major order)
+        do j = 1, n
+            amat(:, j) = sqrt_w(:) * a(:, j)
+        end do
+
+        ! Scale b
+        allocate(b_scaled(m))
+        b_scaled = sqrt_w * b
+
+        ! Solve transformed OLS problem
+        x = stdlib_linalg_${ri}$_lstsq_one(amat, b_scaled, cond=cond, overwrite_a=.true., rank=rank, err=err)
+
+        ! Cleanup
+        if (copy_a) deallocate(amat_alloc)
+        deallocate(b_scaled, sqrt_w)
+
+    end function stdlib_linalg_${ri}$_weighted_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..30b46a1c7 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, weighted_lstsq, generalized_lstsq
     use stdlib_linalg_state, only: linalg_state_type
 
     implicit none (type,external)
@@ -26,6 +26,11 @@ 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("weighted_lstsq_${ri}$",test_weighted_lstsq_${ri}$))
+        call add_test(tests,new_unittest("weighted_lstsq_effect_${ri}$",test_weighted_lstsq_effect_${ri}$))
+        call add_test(tests,new_unittest("weighted_lstsq_negative_${ri}$",test_weighted_lstsq_negative_${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 +144,178 @@ module test_linalg_least_squares
 
     end subroutine test_issue_823
 
+    !-------------------------------------------------------------
+    !-----     Weighted Least-Squares Tests                  -----
+    !-------------------------------------------------------------
+
+    #:for rk,rt,ri in REAL_KINDS_TYPES
+    !> Test basic weighted least-squares
+    subroutine test_weighted_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 = 100*sqrt(epsilon(0.0_${rk}$))
+        ${rt}$ :: A(m,n), b(m)
+        real(${rk}$) :: w(m)
+        ${rt}$, allocatable :: x(:)
+
+        ! Simple test case
+        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}$]
+        w = [1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$]  ! Uniform weights = OLS
+
+        x = weighted_lstsq(w, A, b, err=state)
+
+        call check(error, state%ok(), 'weighted_lstsq failed: '//state%print())
+        if (allocated(error)) return
+        
+        call check(error, size(x)==n, 'weighted_lstsq: wrong solution size')
+        if (allocated(error)) return
+
+        ! Verify residual is small
+        call check(error, norm2(matmul(A,x) - b) < 2.0_${rk}$, 'weighted_lstsq: residual too large')
+        if (allocated(error)) return
+
+    end subroutine test_weighted_lstsq_${ri}$
+
+    !> Test that non-uniform weights change the solution
+    subroutine test_weighted_lstsq_effect_${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 = 100*sqrt(epsilon(0.0_${rk}$))
+        ${rt}$ :: A(m,n), b(m)
+        real(${rk}$) :: w_uniform(m), w_nonuniform(m)
+        ${rt}$, allocatable :: x_uniform(:), x_weighted(:)
+
+        ! Setup problem
+        A(:,1) = 1.0_${rk}$
+        A(:,2) = [1.0_${rk}$, 2.0_${rk}$, 3.0_${rk}$, 4.0_${rk}$]
+        b = [1.0_${rk}$, 3.0_${rk}$, 2.0_${rk}$, 5.0_${rk}$]
+        
+        w_uniform = 1.0_${rk}$
+        w_nonuniform = [10.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$, 10.0_${rk}$]  ! Weight first and last more
+
+        x_uniform = weighted_lstsq(w_uniform, A, b, err=state)
+        call check(error, state%ok(), 'uniform weighted_lstsq failed')
+        if (allocated(error)) return
+
+        x_weighted = weighted_lstsq(w_nonuniform, A, b, err=state)
+        call check(error, state%ok(), 'non-uniform weighted_lstsq failed')
+        if (allocated(error)) return
+
+        ! Solutions should be different
+        call check(error, any(abs(x_uniform - x_weighted) > tol), &
+                   'weighted_lstsq: non-uniform weights should change solution')
+        if (allocated(error)) return
+
+    end subroutine test_weighted_lstsq_effect_${ri}$
+
+    !> Test error on negative weights
+    subroutine test_weighted_lstsq_negative_${ri}$(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        type(linalg_state_type) :: state
+        ${rt}$ :: A(3,2), b(3)
+        real(${rk}$) :: w(3)
+        ${rt}$, allocatable :: x(:)
+
+        A = 1.0_${rk}$
+        b = 1.0_${rk}$
+        w = [-1.0_${rk}$, 1.0_${rk}$, 1.0_${rk}$]  ! Invalid: negative weight!
+
+        x = weighted_lstsq(w, A, b, err=state)
+
+        call check(error, state%error(), 'weighted_lstsq should fail on negative weights')
+        if (allocated(error)) return
+
+    end subroutine test_weighted_lstsq_negative_${ri}$
+
+    #:endfor
+
+    !-------------------------------------------------------------
+    !-----     Generalized Least-Squares Tests               -----
+    !-------------------------------------------------------------
+
+    #: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(:)