diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 082efaefa..04f358fe4 100644 --- a/doc/specs/stdlib_linalg.md +++ b/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 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..638340daf --- /dev/null +++ b/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 diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 475424c91..1af3ab328 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,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 diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index 06169e071..bd8955f12 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 @@ -47,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 @@ -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 diff --git a/src/linalg/stdlib_linalg_least_squares.fypp b/src/linalg/stdlib_linalg_least_squares.fypp index bee76e40f..eaf10d55d 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, 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,172 @@ 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 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(:,:) + !> 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 stop + type(linalg_state_type), optional, intent(out) :: err + !> Result array x(n) + ${rt}$, allocatable :: 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 + 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) + n = size(a, 2, kind=ilp) + p = m ! For GLS, B is m×m + + ! 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:', shape(w, 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 + + 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) + copy_w = .not. optval(overwrite_w, .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/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) + 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) + if (copy_w) deallocate(lmat_alloc) + call linalg_error_handling(err0, err) + return + end if + 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, 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, lda, ldb, err0) + + ! Cleanup + if (copy_a) deallocate(amat_alloc) + if (copy_w) deallocate(lmat_alloc) + deallocate(d, y, work) + + call linalg_error_handling(err0, err) + + 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 7ec41a2e8..a0330493e 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 + 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) @@ -26,6 +27,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("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}$)) + call add_test(tests,new_unittest("solve_generalized_lstsq_${ri}$",test_solve_generalized_lstsq_${ri}$)) #:endfor end subroutine test_least_squares @@ -139,6 +145,206 @@ 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 + ${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 = 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) + + 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 + 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 = 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') + 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}$ + + !> 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 + 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 = 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) + 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^H + ! This is a DENSE matrix (not triangular like Cholesky factor) + ! 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) + 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 = 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) + 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}$ + + !> 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 = 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) + 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 subroutine add_test(tests,new_test) type(unittest_type), allocatable, intent(inout) :: tests(:)