From fbfe4dcd323984389e3c2bfe799862365dba9c1b Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Fri, 23 Jan 2026 08:30:17 +0000 Subject: [PATCH 1/8] fix: rebase onto v0.8.0, resolve merge conflicts --- doc/specs/stdlib_linalg.md | 96 +++++++++ example/linalg/CMakeLists.txt | 2 + example/linalg/example_cholesky_solve.f90 | 27 +++ example/linalg/example_solve_chol.f90 | 27 +++ src/lapack/stdlib_linalg_lapack_aux.fypp | 61 ++++++ src/linalg/stdlib_linalg.fypp | 88 +++++++- src/linalg/stdlib_linalg_solve.fypp | 157 +++++++++++++- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_solve_chol.fypp | 244 ++++++++++++++++++++++ 9 files changed, 698 insertions(+), 6 deletions(-) create mode 100644 example/linalg/example_cholesky_solve.f90 create mode 100644 example/linalg/example_solve_chol.f90 create mode 100644 test/linalg/test_linalg_solve_chol.fypp diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 082efaefa..51933242e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -736,6 +736,102 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_solve3.f90!} ``` +## `solve_chol` - Solves a linear system using pre-computed Cholesky factors (subroutine interface). + +### Status + +Experimental + +### Description + +This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix that has been **previously factorized** using the Cholesky decomposition (via `cholesky`). + +Result vector or array `x` returns the exact solution to within numerical precision, provided that the factorization is correct. +An error is returned if the matrix and right-hand-side have incompatible sizes. +The solver is based on LAPACK's `*POTRS` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):solve_chol(interface)]] `(a, b, x, lower [, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` square array containing the Cholesky-factorized matrix (output of `cholesky`). It is an `intent(in)` argument. + +`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. + +`x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. + +`lower`: Shall be an input `logical` flag. If `.true.`, the lower triangular Cholesky factor (`L`) is stored in `a`. If `.false.`, the upper triangular factor (`U`) is stored. This must match the `lower` flag used during the Cholesky factorization. It is a **required** `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +For a correctly factorized matrix, returns an array value that represents the solution to the linear system of equations. + +Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_solve_chol.f90!} +``` + +## `cholesky_solve` - Solves a linear matrix equation using Cholesky factorization (subroutine interface). + +### Status + +Experimental + +### Description + +This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix. It combines Cholesky factorization and the solve step in a single call. + +Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is positive definite. +An error is returned if the matrix is not positive definite or has incompatible sizes with the right-hand-side. +Use this routine for one-time solves. For repeated solves with the same matrix but different right-hand sides, use `cholesky` followed by `solve_chol` for better performance. +The solver is based on LAPACK's `*POSV` backends. + +### Syntax + +Simple (`Pure`) interface: + +`call ` [[stdlib_linalg(module):cholesky_solve(interface)]] `(a, b, x)` + +Expert (`Pure`) interface: + +`call ` [[stdlib_linalg(module):cholesky_solve(interface)]] `(a, b, x [, lower, overwrite_a, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. + +`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. + +`x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. + +`lower` (optional): Shall be an input `logical` flag. If `.true.` (default), the lower triangular Cholesky factorization is computed. If `.false.`, the upper triangular factorization is computed. It 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 + +For a positive definite matrix, returns an array value that represents the solution to the linear system of equations. + +Raises `LINALG_ERROR` if the matrix is not positive definite. +Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_cholesky_solve.f90!} +``` + ## `lstsq` - Computes the least squares solution to a linear matrix equation. ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1c6c7e42c..1e313cf6f 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -63,5 +63,7 @@ ADD_EXAMPLE(qr_space) ADD_EXAMPLE(pivoting_qr_space) ADD_EXAMPLE(cholesky) ADD_EXAMPLE(chol) +ADD_EXAMPLE(solve_chol) +ADD_EXAMPLE(cholesky_solve) ADD_EXAMPLE(expm) ADD_EXAMPLE(matrix_exp) diff --git a/example/linalg/example_cholesky_solve.f90 b/example/linalg/example_cholesky_solve.f90 new file mode 100644 index 000000000..6e25eafeb --- /dev/null +++ b/example/linalg/example_cholesky_solve.f90 @@ -0,0 +1,27 @@ +program example_cholesky_solve + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: cholesky_solve, linalg_state_type + implicit none + + real(dp) :: A(3,3), b(3), x(3) + type(linalg_state_type) :: state + + ! Symmetric positive definite matrix + A(1,:) = [4.0_dp, 2.0_dp, 2.0_dp] + A(2,:) = [2.0_dp, 5.0_dp, 1.0_dp] + A(3,:) = [2.0_dp, 1.0_dp, 6.0_dp] + + ! Right-hand side + b = [1.0_dp, 2.0_dp, 3.0_dp] + + ! One-shot factorization and solve (A is preserved by default) + call cholesky_solve(A, b, x, lower=.true., err=state) + if (state%error()) error stop state%print() + + print '("Solution: ",*(f8.4,1x))', x + + ! For performance-critical code, use overwrite_a=.true. + ! to avoid internal allocation (but A will be destroyed) + ! call cholesky_solve(A, b, x, lower=.true., overwrite_a=.true., err=state) + +end program example_cholesky_solve diff --git a/example/linalg/example_solve_chol.f90 b/example/linalg/example_solve_chol.f90 new file mode 100644 index 000000000..7d5be06d0 --- /dev/null +++ b/example/linalg/example_solve_chol.f90 @@ -0,0 +1,27 @@ +program example_solve_chol + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: cholesky, solve_chol, linalg_state_type + implicit none + + real(dp) :: A(3,3), L(3,3), b(3), x(3) + type(linalg_state_type) :: state + + ! Symmetric positive definite matrix + A(1,:) = [4.0_dp, 2.0_dp, 2.0_dp] + A(2,:) = [2.0_dp, 5.0_dp, 1.0_dp] + A(3,:) = [2.0_dp, 1.0_dp, 6.0_dp] + + ! Right-hand side + b = [1.0_dp, 2.0_dp, 3.0_dp] + + ! Compute Cholesky factorization: A = L * L^T + call cholesky(A, L, lower=.true., err=state) + if (state%error()) error stop state%print() + + ! Solve using pre-computed Cholesky factors + call solve_chol(L, b, x, lower=.true., err=state) + if (state%error()) error stop state%print() + + print '("Solution: ",*(f8.4,1x))', x + +end program example_solve_chol diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 475424c91..cc75d9cc3 100644 --- a/src/lapack/stdlib_linalg_lapack_aux.fypp +++ b/src/lapack/stdlib_linalg_lapack_aux.fypp @@ -41,6 +41,8 @@ module stdlib_linalg_lapack_aux public :: stdlib_selctg_${ri}$ #:endfor public :: handle_potrf_info + public :: handle_potrs_info + public :: handle_posv_info public :: handle_getri_info public :: handle_gesdd_info public :: handle_gesv_info @@ -1323,6 +1325,65 @@ module stdlib_linalg_lapack_aux end subroutine handle_potrf_info + ! Cholesky solve (triangular solve with pre-computed factors) + elemental subroutine handle_potrs_info(this,info,triangle,n,nrhs,lda,ldb,err) + character(len=*), intent(in) :: this + character, intent(in) :: triangle + integer(ilp), intent(in) :: info,n,nrhs,lda,ldb + type(linalg_state_type), intent(out) :: err + + ! Process output + select case (info) + case (0) + ! Success + case (-1) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', & + triangle,'. should be U/L') + case (-2) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n) + case (-3) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size nrhs=',nrhs) + case (-5) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid lda=',lda,': should be >=',n) + case (-7) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid ldb=',ldb,': should be >=',n) + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select + + end subroutine handle_potrs_info + + ! Cholesky factorization and solve (combined) + elemental subroutine handle_posv_info(this,info,triangle,n,nrhs,lda,ldb,err) + character(len=*), intent(in) :: this + character, intent(in) :: triangle + integer(ilp), intent(in) :: info,n,nrhs,lda,ldb + type(linalg_state_type), intent(out) :: err + + ! Process output + select case (info) + case (0) + ! Success + case (-1) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', & + triangle,'. should be U/L') + case (-2) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n) + case (-3) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size nrhs=',nrhs) + case (-5) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid lda=',lda,': should be >=',n) + case (-7) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid ldb=',ldb,': should be >=',n) + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'matrix is not positive definite: ', & + 'leading minor of order',info,' is not positive definite') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select + + end subroutine handle_posv_info + elemental subroutine handle_getri_info(this,info,lda,n,err) character(len=*), intent(in) :: this integer(ilp), intent(in) :: info,lda,n diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index 06169e071..dfe376f35 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -44,7 +44,9 @@ module stdlib_linalg public :: mnorm public :: get_norm public :: solve - public :: solve_lu + public :: solve_lu + public :: solve_chol + public :: cholesky_solve public :: solve_lstsq public :: solve_constrained_lstsq public :: trace @@ -457,9 +459,91 @@ module stdlib_linalg !> [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_lu_${ndsuf}$ + #:endfor + #:endfor + end interface solve_lu + + ! Solve linear system Ax = b using pre-computed Cholesky decomposition (subroutine interface) + interface solve_chol + !! version: experimental + !! + !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a + !! symmetric positive definite matrix \( A \) that has been pre-factorized using Cholesky. + !! ([Specification](../page/specs/stdlib_linalg.html#solve_chol-solve-a-linear-system-using-cholesky-factors)) + !! + !!### Summary + !! Subroutine interface for solving a linear system using pre-computed Cholesky factors. + !! + !!### Description + !! + !! This interface provides methods for computing the solution of a linear matrix system using + !! Cholesky factors. Supported data types include `real` and `complex`. Preallocated space + !! for the solution vector `x` is user-provided. The `lower` argument is REQUIRED and must + !! match the `lower` used during the Cholesky factorization. + !! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`) + !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems. + !! + !!@note The solution is based on LAPACK's `*POTRS` routines. + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + pure module subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$(a,b,x,lower,err) + !> Input matrix a[n,n] containing Cholesky factors from cholesky + ${rt}$, intent(in) :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> Is the lower triangular factor stored? (REQUIRED) + logical(lk), intent(in) :: lower + !> [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_chol_${ndsuf}$ #:endfor #:endfor - end interface solve_lu + end interface solve_chol + + ! One-shot Cholesky factorization and solve (convenience wrapper using POSV) + interface cholesky_solve + !! version: experimental + !! + !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a + !! symmetric positive definite matrix \( A \). Combines factorization and solve in one call. + !! ([Specification](../page/specs/stdlib_linalg.html#cholesky_solve-one-shot-cholesky-solve)) + !! + !!### Summary + !! One-shot factorization and solve for SPD systems (wraps LAPACK POSV). + !! + !!### Description + !! + !! This interface computes both the Cholesky factorization and solves the linear system + !! in a single call. Use this for one-time solves. For repeated solves with the same + !! matrix but different RHS, use `cholesky` + `solve_chol` for better performance. + !! Supported data types include `real` and `complex`. + !! By default, A is not overwritten. Set `overwrite_a=.true.` to allow in-place + !! factorization for better performance. + !! + !!@note The solution is based on LAPACK's `*POSV` routines. + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + pure module subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$(a,b,x,lower,overwrite_a,err) + !> Input SPD matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] Use lower triangular factorization? Default = .true. + logical(lk), optional, intent(in) :: lower + !> [optional] Can A data be overwritten and destroyed? Default = .false. + 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 + end subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$ + #:endfor + #:endfor + end interface cholesky_solve ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. interface lstsq diff --git a/src/linalg/stdlib_linalg_solve.fypp b/src/linalg/stdlib_linalg_solve.fypp index b4df33133..a25cf034e 100644 --- a/src/linalg/stdlib_linalg_solve.fypp +++ b/src/linalg/stdlib_linalg_solve.fypp @@ -7,8 +7,8 @@ submodule (stdlib_linalg) stdlib_linalg_solve !! Solve linear system Ax=b use stdlib_linalg_constants - use stdlib_linalg_lapack, only: gesv - use stdlib_linalg_lapack_aux, only: handle_gesv_info + use stdlib_linalg_lapack, only: gesv, potrs, posv + use stdlib_linalg_lapack_aux, only: handle_gesv_info, handle_potrs_info, handle_posv_info use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR implicit none @@ -140,7 +140,156 @@ submodule (stdlib_linalg) stdlib_linalg_solve end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$ - #:endfor - #:endfor + #:endfor + #:endfor + + !--------------------------------------------------------------------------- + !> solve_chol: Solve using PRE-COMPUTED Cholesky factors (POTRS) + !--------------------------------------------------------------------------- + + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + !> Solve the linear system A*x = b using pre-computed Cholesky factorization + pure module subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$(a,b,x,lower,err) + !> Cholesky-factorized matrix a[n,n] + ${rt}$, intent(in) :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> Is the lower triangular factor stored? (REQUIRED) + logical(lk), intent(in) :: lower + !> [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) :: lda,n,ldb,ldx,nrhs,nrhsx,info + character :: triangle + ${rt}$, pointer :: xmat(:,:) + + ! Problem sizes + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + ldb = size(b,1,kind=ilp) + nrhs = size(b,kind=ilp)/ldb + ldx = size(x,1,kind=ilp) + nrhsx = size(x,kind=ilp)/ldx + + ! Set triangle based on lower flag (REQUIRED argument) + triangle = merge('L','U',lower) + + ! Validate dimensions + if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + 'b=',[ldb,nrhs],' x=',[ldx,nrhsx]) + call linalg_error_handling(err0,err) + return + end if + + ! Copy RHS to solution array (POTRS overwrites with solution) + x = b + + ! Create 2D pointer for LAPACK call + xmat(1:n,1:nrhs) => x + + ! Solve the system using LAPACK POTRS + call potrs(triangle,n,nrhs,a,lda,xmat,n,info) + + ! Handle errors using standard handler + call handle_potrs_info(this,info,triangle,n,nrhs,lda,n,err0) + + ! Process output and return + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$ + + #:endfor + #:endfor + + !--------------------------------------------------------------------------- + !> cholesky_solve: One-shot factorize + solve (POSV) + !--------------------------------------------------------------------------- + + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + !> Factorize and solve A*x = b in one call (uses LAPACK POSV) + pure module subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$(a,b,x,lower,overwrite_a,err) + !> Input SPD matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] Use lower triangular factorization? Default = .true. + logical(lk), optional, intent(in) :: lower + !> [optional] Can A data be overwritten and destroyed? Default = .false. + 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 + + ! Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: lda,n,ldb,ldx,nrhs,nrhsx,info + logical(lk) :: lower_,copy_a + character :: triangle + ${rt}$, pointer :: xmat(:,:),amat(:,:) + + ! Problem sizes + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + ldb = size(b,1,kind=ilp) + nrhs = size(b,kind=ilp)/ldb + ldx = size(x,1,kind=ilp) + nrhsx = size(x,kind=ilp)/ldx + + ! Default: use lower triangular + lower_ = .true._lk + if (present(lower)) lower_ = lower + triangle = merge('L','U',lower_) + + ! Can A be overwritten? By default, do not overwrite + if (present(overwrite_a)) then + copy_a = .not.overwrite_a + else + copy_a = .true._lk + endif + + ! Validate dimensions + if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + 'b=',[ldb,nrhs],' x=',[ldx,nrhsx]) + call linalg_error_handling(err0,err) + return + end if + + ! Initialize a matrix temporary + if (copy_a) then + allocate(amat(lda,n),source=a) + else + amat => a + endif + + ! Copy RHS to solution array (POSV overwrites with solution) + x = b + + ! Create 2D pointer for LAPACK call + xmat(1:n,1:nrhs) => x + + ! Factorize AND solve using LAPACK POSV + call posv(triangle,n,nrhs,amat,lda,xmat,n,info) + + ! Handle errors using standard handler + call handle_posv_info(this,info,triangle,n,nrhs,lda,n,err0) + + if (copy_a) deallocate(amat) + + ! Process output and return + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$ + + #:endfor + #:endfor end submodule stdlib_linalg_solve diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index f27dedd4f..097591893 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -19,6 +19,7 @@ set( "test_linalg_sparse.fypp" "test_linalg_specialmatrices.fypp" "test_linalg_cholesky.fypp" + "test_linalg_solve_chol.fypp" "test_linalg_expm.fypp" ) @@ -42,6 +43,7 @@ ADDTEST(linalg_pseudoinverse) ADDTEST(linalg_norm) ADDTEST(linalg_mnorm) ADDTEST(linalg_solve) +ADDTEST(linalg_solve_chol) ADDTEST(linalg_lstsq) ADDTEST(linalg_constrained_lstsq) ADDTEST(linalg_qr) diff --git a/test/linalg/test_linalg_solve_chol.fypp b/test/linalg/test_linalg_solve_chol.fypp new file mode 100644 index 000000000..2bcccf27a --- /dev/null +++ b/test/linalg/test_linalg_solve_chol.fypp @@ -0,0 +1,244 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test solve_chol and cholesky_solve +module test_linalg_solve_chol + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: cholesky, solve_chol, cholesky_solve + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + private + + public :: test_solve_chol_factorization + + contains + + !> Cholesky solve tests + subroutine test_solve_chol_factorization(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in RC_KINDS_TYPES + call add_test(tests,new_unittest("solve_chol_${ri}$",test_solve_chol_${ri}$)) + call add_test(tests,new_unittest("cholesky_solve_${ri}$",test_cholesky_solve_${ri}$)) + call add_test(tests,new_unittest("cholesky_solve_overwrite_${ri}$",test_cholesky_solve_overwrite_${ri}$)) + call add_test(tests,new_unittest("solve_chol_multi_rhs_${ri}$",test_solve_chol_multi_rhs_${ri}$)) + #:endfor + + end subroutine test_solve_chol_factorization + + !> Test solve_chol with pre-computed Cholesky factors + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_chol_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n,n), l(n,n), b(n), x(n), x_expected(n) + type(linalg_state_type) :: state + + ! Set symmetric positive definite matrix + a(1,:) = [4, 2, 2] + a(2,:) = [2, 5, 1] + a(3,:) = [2, 1, 6] + + ! Known solution + x_expected = [1, 2, 3] + + ! Compute RHS: b = A * x_expected + b = matmul(a, x_expected) + + ! Compute Cholesky factorization + call cholesky(a, l, lower=.true., err=state) + + call check(error, state%ok(), 'cholesky factorization failed: '//state%print()) + if (allocated(error)) return + + ! Solve using Cholesky factors + call solve_chol(l, b, x, lower=.true., err=state) + + call check(error, state%ok(), 'solve_chol failed: '//state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(x - x_expected) < tol), 'solve_chol: solution mismatch') + if (allocated(error)) return + + end subroutine test_solve_chol_${ri}$ + + #:endfor + + !> Test cholesky_solve (one-shot) - default preserves A + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_cholesky_solve_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n,n), a_copy(n,n), b(n), x(n), x_expected(n) + type(linalg_state_type) :: state + + ! Set symmetric positive definite matrix + a(1,:) = [4, 2, 2] + a(2,:) = [2, 5, 1] + a(3,:) = [2, 1, 6] + a_copy = a + + ! Known solution + x_expected = [1, 2, 3] + + ! Compute RHS: b = A * x_expected + b = matmul(a, x_expected) + + ! One-shot solve with default overwrite_a=.false. (A should be preserved) + call cholesky_solve(a, b, x, lower=.true., err=state) + + call check(error, state%ok(), 'cholesky_solve failed: '//state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(x - x_expected) < tol), 'cholesky_solve: solution mismatch') + if (allocated(error)) return + + ! Check that A is preserved (default behavior) + call check(error, all(abs(a - a_copy) < tol), 'cholesky_solve: A should be preserved by default') + if (allocated(error)) return + + end subroutine test_cholesky_solve_${ri}$ + + #:endfor + + !> Test cholesky_solve with overwrite_a=.true. + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_cholesky_solve_overwrite_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n,n), a_copy(n,n), b(n), x(n), x_expected(n) + type(linalg_state_type) :: state + + ! Set symmetric positive definite matrix + a(1,:) = [4, 2, 2] + a(2,:) = [2, 5, 1] + a(3,:) = [2, 1, 6] + a_copy = a + + ! Known solution + x_expected = [1, 2, 3] + + ! Compute RHS: b = A * x_expected + b = matmul(a, x_expected) + + ! One-shot solve with overwrite_a=.true. (A will be destroyed) + call cholesky_solve(a, b, x, lower=.true., overwrite_a=.true., err=state) + + call check(error, state%ok(), 'cholesky_solve overwrite failed: '//state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(x - x_expected) < tol), 'cholesky_solve overwrite: solution mismatch') + if (allocated(error)) return + + ! Check that A was overwritten (not equal to original) + call check(error, any(abs(a - a_copy) > tol), 'cholesky_solve: A should be overwritten with overwrite_a=.true.') + if (allocated(error)) return + + end subroutine test_cholesky_solve_overwrite_${ri}$ + + #:endfor + + !> Test solve_chol with multiple RHS + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_chol_multi_rhs_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + integer(ilp), parameter :: nrhs = 2_ilp + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n,n), l(n,n), b(n,nrhs), x(n,nrhs), x_expected(n,nrhs) + type(linalg_state_type) :: state + + ! Set symmetric positive definite matrix + a(1,:) = [4, 2, 2] + a(2,:) = [2, 5, 1] + a(3,:) = [2, 1, 6] + + ! Known solutions (two RHS) + x_expected(:,1) = [1, 2, 3] + x_expected(:,2) = [4, 5, 6] + + ! Compute RHS: B = A * X_expected + b = matmul(a, x_expected) + + ! Compute Cholesky factorization + call cholesky(a, l, lower=.true., err=state) + + call check(error, state%ok(), 'cholesky factorization failed: '//state%print()) + if (allocated(error)) return + + ! Solve using Cholesky factors + call solve_chol(l, b, x, lower=.true., err=state) + + call check(error, state%ok(), 'solve_chol multi-rhs failed: '//state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(x - x_expected) < tol), 'solve_chol multi-rhs: solution mismatch') + if (allocated(error)) return + + end subroutine test_solve_chol_multi_rhs_${ri}$ + + #:endfor + + ! gcc-15 bugfix utility + subroutine add_test(tests,new_test) + type(unittest_type), allocatable, intent(inout) :: tests(:) + type(unittest_type), intent(in) :: new_test + + integer :: n + type(unittest_type), allocatable :: new_tests(:) + + if (allocated(tests)) then + n = size(tests) + else + n = 0 + end if + + allocate(new_tests(n+1)) + if (n>0) new_tests(1:n) = tests(1:n) + new_tests(1+n) = new_test + call move_alloc(from=new_tests,to=tests) + + end subroutine add_test + +end module test_linalg_solve_chol + +program test_solve_chol + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_solve_chol, only : test_solve_chol_factorization + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_solve_chol", test_solve_chol_factorization) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_solve_chol From 840d0eebd3a279e4ec7e8ab6d68cef5b14792d1c Mon Sep 17 00:00:00 2001 From: Amrinder Singh Date: Fri, 23 Jan 2026 22:07:50 +0530 Subject: [PATCH 2/8] Update doc/specs/stdlib_linalg.md Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 51933242e..f737e761e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -806,7 +806,7 @@ Expert (`Pure`) interface: ### Arguments -`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. +`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is an `intent(inout)` argument. By default (`overwrite_a=.false.`) its contents are preserved; if `overwrite_a=.true.`, it is used as temporary storage and its contents are destroyed by the call. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. From a90185ba697130c8730c94ba37b10dbc9ec92c6a Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Fri, 23 Jan 2026 16:41:01 +0000 Subject: [PATCH 3/8] test: add upper triangular tests for solve_chol --- test/linalg/test_linalg_solve_chol.fypp | 83 +++++++++++++++++++++++++ 1 file changed, 83 insertions(+) diff --git a/test/linalg/test_linalg_solve_chol.fypp b/test/linalg/test_linalg_solve_chol.fypp index 2bcccf27a..550f837b5 100644 --- a/test/linalg/test_linalg_solve_chol.fypp +++ b/test/linalg/test_linalg_solve_chol.fypp @@ -23,7 +23,9 @@ module test_linalg_solve_chol #:for rk,rt,ri in RC_KINDS_TYPES call add_test(tests,new_unittest("solve_chol_${ri}$",test_solve_chol_${ri}$)) + call add_test(tests,new_unittest("solve_chol_upper_${ri}$",test_solve_chol_upper_${ri}$)) call add_test(tests,new_unittest("cholesky_solve_${ri}$",test_cholesky_solve_${ri}$)) + call add_test(tests,new_unittest("cholesky_solve_upper_${ri}$",test_cholesky_solve_upper_${ri}$)) call add_test(tests,new_unittest("cholesky_solve_overwrite_${ri}$",test_cholesky_solve_overwrite_${ri}$)) call add_test(tests,new_unittest("solve_chol_multi_rhs_${ri}$",test_solve_chol_multi_rhs_${ri}$)) #:endfor @@ -71,6 +73,47 @@ module test_linalg_solve_chol #:endfor + !> Test solve_chol with upper triangular factors (lower=.false.) + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_chol_upper_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n,n), u(n,n), b(n), x(n), x_expected(n) + type(linalg_state_type) :: state + + ! Set symmetric positive definite matrix + a(1,:) = [4, 2, 2] + a(2,:) = [2, 5, 1] + a(3,:) = [2, 1, 6] + + ! Known solution + x_expected = [1, 2, 3] + + ! Compute RHS: b = A * x_expected + b = matmul(a, x_expected) + + ! Compute Cholesky factorization (upper triangular: A = U^T * U) + call cholesky(a, u, lower=.false., err=state) + + call check(error, state%ok(), 'cholesky factorization (upper) failed: '//state%print()) + if (allocated(error)) return + + ! Solve using upper Cholesky factors + call solve_chol(u, b, x, lower=.false., err=state) + + call check(error, state%ok(), 'solve_chol (upper) failed: '//state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(x - x_expected) < tol), 'solve_chol (upper): solution mismatch') + if (allocated(error)) return + + end subroutine test_solve_chol_upper_${ri}$ + + #:endfor + !> Test cholesky_solve (one-shot) - default preserves A #:for rk,rt,ri in RC_KINDS_TYPES subroutine test_cholesky_solve_${ri}$(error) @@ -111,6 +154,46 @@ module test_linalg_solve_chol #:endfor + !> Test cholesky_solve with upper triangular (lower=.false.) + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_cholesky_solve_upper_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n,n), a_copy(n,n), b(n), x(n), x_expected(n) + type(linalg_state_type) :: state + + ! Set symmetric positive definite matrix + a(1,:) = [4, 2, 2] + a(2,:) = [2, 5, 1] + a(3,:) = [2, 1, 6] + a_copy = a + + ! Known solution + x_expected = [1, 2, 3] + + ! Compute RHS: b = A * x_expected + b = matmul(a, x_expected) + + ! One-shot solve with upper triangular (A should be preserved by default) + call cholesky_solve(a, b, x, lower=.false., err=state) + + call check(error, state%ok(), 'cholesky_solve (upper) failed: '//state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(x - x_expected) < tol), 'cholesky_solve (upper): solution mismatch') + if (allocated(error)) return + + ! Check that A is preserved (default behavior) + call check(error, all(abs(a - a_copy) < tol), 'cholesky_solve (upper): A should be preserved') + if (allocated(error)) return + + end subroutine test_cholesky_solve_upper_${ri}$ + + #:endfor + !> Test cholesky_solve with overwrite_a=.true. #:for rk,rt,ri in RC_KINDS_TYPES subroutine test_cholesky_solve_overwrite_${ri}$(error) From 921ba4b19d6f6754c3bb4cf30ac4ef498f4ed8fa Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Tue, 27 Jan 2026 16:55:24 +0000 Subject: [PATCH 4/8] refactor: use uplo/optval, add error tests --- src/linalg/stdlib_linalg_solve.fypp | 27 +++++------- test/linalg/test_linalg_solve_chol.fypp | 56 +++++++++++++++++++++++++ 2 files changed, 67 insertions(+), 16 deletions(-) diff --git a/src/linalg/stdlib_linalg_solve.fypp b/src/linalg/stdlib_linalg_solve.fypp index a25cf034e..76a6baa70 100644 --- a/src/linalg/stdlib_linalg_solve.fypp +++ b/src/linalg/stdlib_linalg_solve.fypp @@ -165,7 +165,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve ! Local variables type(linalg_state_type) :: err0 integer(ilp) :: lda,n,ldb,ldx,nrhs,nrhsx,info - character :: triangle + character :: uplo ${rt}$, pointer :: xmat(:,:) ! Problem sizes @@ -176,8 +176,8 @@ submodule (stdlib_linalg) stdlib_linalg_solve ldx = size(x,1,kind=ilp) nrhsx = size(x,kind=ilp)/ldx - ! Set triangle based on lower flag (REQUIRED argument) - triangle = merge('L','U',lower) + ! Set uplo based on lower flag (REQUIRED argument) + uplo = merge('L','U',lower) ! Validate dimensions if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then @@ -194,10 +194,10 @@ submodule (stdlib_linalg) stdlib_linalg_solve xmat(1:n,1:nrhs) => x ! Solve the system using LAPACK POTRS - call potrs(triangle,n,nrhs,a,lda,xmat,n,info) + call potrs(uplo,n,nrhs,a,lda,xmat,n,info) ! Handle errors using standard handler - call handle_potrs_info(this,info,triangle,n,nrhs,lda,n,err0) + call handle_potrs_info(this,info,uplo,n,nrhs,lda,n,err0) ! Process output and return call linalg_error_handling(err0,err) @@ -232,7 +232,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve type(linalg_state_type) :: err0 integer(ilp) :: lda,n,ldb,ldx,nrhs,nrhsx,info logical(lk) :: lower_,copy_a - character :: triangle + character :: uplo ${rt}$, pointer :: xmat(:,:),amat(:,:) ! Problem sizes @@ -244,16 +244,11 @@ submodule (stdlib_linalg) stdlib_linalg_solve nrhsx = size(x,kind=ilp)/ldx ! Default: use lower triangular - lower_ = .true._lk - if (present(lower)) lower_ = lower - triangle = merge('L','U',lower_) + lower_ = optval(lower, .true._lk) + uplo = merge('L','U',lower_) ! Can A be overwritten? By default, do not overwrite - if (present(overwrite_a)) then - copy_a = .not.overwrite_a - else - copy_a = .true._lk - endif + copy_a = .not. optval(overwrite_a, .false._lk) ! Validate dimensions if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then @@ -277,10 +272,10 @@ submodule (stdlib_linalg) stdlib_linalg_solve xmat(1:n,1:nrhs) => x ! Factorize AND solve using LAPACK POSV - call posv(triangle,n,nrhs,amat,lda,xmat,n,info) + call posv(uplo,n,nrhs,amat,lda,xmat,n,info) ! Handle errors using standard handler - call handle_posv_info(this,info,triangle,n,nrhs,lda,n,err0) + call handle_posv_info(this,info,uplo,n,nrhs,lda,n,err0) if (copy_a) deallocate(amat) diff --git a/test/linalg/test_linalg_solve_chol.fypp b/test/linalg/test_linalg_solve_chol.fypp index 550f837b5..bbdd17bda 100644 --- a/test/linalg/test_linalg_solve_chol.fypp +++ b/test/linalg/test_linalg_solve_chol.fypp @@ -28,6 +28,8 @@ module test_linalg_solve_chol call add_test(tests,new_unittest("cholesky_solve_upper_${ri}$",test_cholesky_solve_upper_${ri}$)) call add_test(tests,new_unittest("cholesky_solve_overwrite_${ri}$",test_cholesky_solve_overwrite_${ri}$)) call add_test(tests,new_unittest("solve_chol_multi_rhs_${ri}$",test_solve_chol_multi_rhs_${ri}$)) + call add_test(tests,new_unittest("cholesky_solve_indefinite_${ri}$",test_cholesky_solve_indefinite_${ri}$)) + call add_test(tests,new_unittest("cholesky_solve_semidefinite_${ri}$",test_cholesky_solve_semidefinite_${ri}$)) #:endfor end subroutine test_solve_chol_factorization @@ -277,6 +279,60 @@ module test_linalg_solve_chol #:endfor + !> Test cholesky_solve with symmetric indefinite matrix (should fail) + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_cholesky_solve_indefinite_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 2_ilp + ${rt}$ :: a(n,n), b(n), x(n) + type(linalg_state_type) :: state + + ! Set symmetric INDEFINITE matrix (eigenvalues: 3, -1) + a(1,:) = [1, 2] + a(2,:) = [2, 1] + + ! Arbitrary RHS + b = [1, 1] + + ! cholesky_solve should fail for indefinite matrix + call cholesky_solve(a, b, x, lower=.true., err=state) + + ! Check that it failed (not positive definite) + call check(error, state%error(), 'cholesky_solve should fail for indefinite matrix') + if (allocated(error)) return + + end subroutine test_cholesky_solve_indefinite_${ri}$ + + #:endfor + + !> Test cholesky_solve with symmetric positive semi-definite matrix (should fail) + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_cholesky_solve_semidefinite_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 2_ilp + ${rt}$ :: a(n,n), b(n), x(n) + type(linalg_state_type) :: state + + ! Set symmetric positive SEMI-DEFINITE matrix (eigenvalues: 2, 0 - singular) + a(1,:) = [1, 1] + a(2,:) = [1, 1] + + ! Arbitrary RHS + b = [1, 1] + + ! cholesky_solve should fail for semi-definite (singular) matrix + call cholesky_solve(a, b, x, lower=.true., err=state) + + ! Check that it failed (not strictly positive definite) + call check(error, state%error(), 'cholesky_solve should fail for positive semi-definite matrix') + if (allocated(error)) return + + end subroutine test_cholesky_solve_semidefinite_${ri}$ + + #:endfor + ! gcc-15 bugfix utility subroutine add_test(tests,new_test) type(unittest_type), allocatable, intent(inout) :: tests(:) From f006e138d567c58233813ac7cfe004f19f8542e8 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Wed, 28 Jan 2026 10:07:14 +0000 Subject: [PATCH 5/8] Set LINALG_SUCCESS in handle_potrs_info and handle_posv_info --- src/lapack/stdlib_linalg_lapack_aux.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index cc75d9cc3..63e6ae0f7 100644 --- a/src/lapack/stdlib_linalg_lapack_aux.fypp +++ b/src/lapack/stdlib_linalg_lapack_aux.fypp @@ -1335,7 +1335,7 @@ module stdlib_linalg_lapack_aux ! Process output select case (info) case (0) - ! Success + err%state = LINALG_SUCCESS case (-1) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', & triangle,'. should be U/L') @@ -1363,7 +1363,7 @@ module stdlib_linalg_lapack_aux ! Process output select case (info) case (0) - ! Success + err%state = LINALG_SUCCESS case (-1) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', & triangle,'. should be U/L') From ecb81de6b3122a66e6d4be4fc1f8d1a6343ee1fd Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Wed, 28 Jan 2026 15:34:49 +0000 Subject: [PATCH 6/8] Refactor Cholesky API: solve_chol (one-shot), solve_lower_chol, solve_upper_chol --- doc/specs/stdlib_linalg.md | 91 +++++++++----- example/linalg/example_cholesky_solve.f90 | 28 +++-- example/linalg/example_solve_chol.f90 | 16 +-- src/linalg/stdlib_linalg.fypp | 108 ++++++++++------ src/linalg/stdlib_linalg_solve.fypp | 144 +++++++++++++++------- test/linalg/test_linalg_solve_chol.fypp | 122 +++++++++--------- 6 files changed, 318 insertions(+), 191 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index f737e761e..6ec1a1770 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -736,7 +736,7 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_solve3.f90!} ``` -## `solve_chol` - Solves a linear system using pre-computed Cholesky factors (subroutine interface). +## `solve_chol` - Solves a linear system using Cholesky factorization (one-shot interface). ### Status @@ -744,32 +744,42 @@ Experimental ### Description -This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix that has been **previously factorized** using the Cholesky decomposition (via `cholesky`). +This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix. It combines Cholesky factorization and the solve step in a single call. -Result vector or array `x` returns the exact solution to within numerical precision, provided that the factorization is correct. -An error is returned if the matrix and right-hand-side have incompatible sizes. -The solver is based on LAPACK's `*POTRS` backends. +Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is positive definite. +An error is returned if the matrix is not positive definite or has incompatible sizes with the right-hand-side. +Use this routine for one-time solves. For repeated solves with the same matrix but different right-hand sides, use `cholesky` followed by `solve_lower_chol`/`solve_upper_chol` for better performance. +The solver is based on LAPACK's `*POSV` backends. ### Syntax -`call ` [[stdlib_linalg(module):solve_chol(interface)]] `(a, b, x, lower [, err])` +Simple (`Pure`) interface: + +`call ` [[stdlib_linalg(module):solve_chol(interface)]] `(a, b, x)` + +Expert (`Pure`) interface: + +`call ` [[stdlib_linalg(module):solve_chol(interface)]] `(a, b, x [, lower, overwrite_a, err])` ### Arguments -`a`: Shall be a rank-2 `real` or `complex` square array containing the Cholesky-factorized matrix (output of `cholesky`). It is an `intent(in)` argument. +`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is an `intent(inout)` argument. By default (`overwrite_a=.false.`) its contents are preserved; if `overwrite_a=.true.`, it is used as temporary storage and its contents are destroyed by the call. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. -`lower`: Shall be an input `logical` flag. If `.true.`, the lower triangular Cholesky factor (`L`) is stored in `a`. If `.false.`, the upper triangular factor (`U`) is stored. This must match the `lower` flag used during the Cholesky factorization. It is a **required** `intent(in)` argument. +`lower` (optional): Shall be an input `logical` flag. If `.true.` (default), the lower triangular Cholesky factorization is computed. If `.false.`, the upper triangular factorization is computed. It 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 -For a correctly factorized matrix, returns an array value that represents the solution to the linear system of equations. +For a positive definite matrix, returns an array value that represents the solution to the linear system of equations. +Raises `LINALG_ERROR` if the matrix is not positive definite. Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. @@ -779,7 +789,7 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_solve_chol.f90!} ``` -## `cholesky_solve` - Solves a linear matrix equation using Cholesky factorization (subroutine interface). +## `solve_lower_chol` - Solves a linear system using pre-computed lower Cholesky factor. ### Status @@ -787,42 +797,30 @@ Experimental ### Description -This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix. It combines Cholesky factorization and the solve step in a single call. +This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix that has been **previously factorized** using the Cholesky decomposition (via `cholesky` with `lower=.true.`). -Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is positive definite. -An error is returned if the matrix is not positive definite or has incompatible sizes with the right-hand-side. -Use this routine for one-time solves. For repeated solves with the same matrix but different right-hand sides, use `cholesky` followed by `solve_chol` for better performance. -The solver is based on LAPACK's `*POSV` backends. +Result vector or array `x` returns the exact solution to within numerical precision, provided that the factorization is correct. +An error is returned if the matrix and right-hand-side have incompatible sizes. +The solver is based on LAPACK's `*POTRS` backends. ### Syntax -Simple (`Pure`) interface: - -`call ` [[stdlib_linalg(module):cholesky_solve(interface)]] `(a, b, x)` - -Expert (`Pure`) interface: - -`call ` [[stdlib_linalg(module):cholesky_solve(interface)]] `(a, b, x [, lower, overwrite_a, err])` +`call ` [[stdlib_linalg(module):solve_lower_chol(interface)]] `(a, b, x [, err])` ### Arguments -`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is an `intent(inout)` argument. By default (`overwrite_a=.false.`) its contents are preserved; if `overwrite_a=.true.`, it is used as temporary storage and its contents are destroyed by the call. +`a`: Shall be a rank-2 `real` or `complex` square array containing the **lower** Cholesky factor `L` (output of `cholesky(..., lower=.true.)`). It is an `intent(in)` argument. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. -`lower` (optional): Shall be an input `logical` flag. If `.true.` (default), the lower triangular Cholesky factorization is computed. If `.false.`, the upper triangular factorization is computed. It 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 -For a positive definite matrix, returns an array value that represents the solution to the linear system of equations. +For a correctly factorized matrix, returns an array value that represents the solution to the linear system of equations. -Raises `LINALG_ERROR` if the matrix is not positive definite. Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. @@ -832,6 +830,41 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_cholesky_solve.f90!} ``` +## `solve_upper_chol` - Solves a linear system using pre-computed upper Cholesky factor. + +### Status + +Experimental + +### Description + +This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix that has been **previously factorized** using the Cholesky decomposition (via `cholesky` with `lower=.false.`). + +Result vector or array `x` returns the exact solution to within numerical precision, provided that the factorization is correct. +An error is returned if the matrix and right-hand-side have incompatible sizes. +The solver is based on LAPACK's `*POTRS` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):solve_upper_chol(interface)]] `(a, b, x [, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` square array containing the **upper** Cholesky factor `U` (output of `cholesky(..., lower=.false.)`). It is an `intent(in)` argument. + +`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. + +`x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +For a correctly factorized matrix, returns an array value that represents the solution to the linear system of equations. + +Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + ## `lstsq` - Computes the least squares solution to a linear matrix equation. ### Status diff --git a/example/linalg/example_cholesky_solve.f90 b/example/linalg/example_cholesky_solve.f90 index 6e25eafeb..2b4fdf952 100644 --- a/example/linalg/example_cholesky_solve.f90 +++ b/example/linalg/example_cholesky_solve.f90 @@ -1,9 +1,12 @@ +! Example: solve_lower_chol - Solve using pre-computed Cholesky factors +! For repeated solves with the same matrix, pre-compute the factorization +! once and reuse it for better performance. program example_cholesky_solve use stdlib_linalg_constants, only: dp - use stdlib_linalg, only: cholesky_solve, linalg_state_type + use stdlib_linalg, only: cholesky, solve_lower_chol, linalg_state_type implicit none - real(dp) :: A(3,3), b(3), x(3) + real(dp) :: A(3,3), L(3,3), b1(3), b2(3), x(3) type(linalg_state_type) :: state ! Symmetric positive definite matrix @@ -11,17 +14,20 @@ program example_cholesky_solve A(2,:) = [2.0_dp, 5.0_dp, 1.0_dp] A(3,:) = [2.0_dp, 1.0_dp, 6.0_dp] - ! Right-hand side - b = [1.0_dp, 2.0_dp, 3.0_dp] - - ! One-shot factorization and solve (A is preserved by default) - call cholesky_solve(A, b, x, lower=.true., err=state) + ! Compute Cholesky factorization once: A = L * L^T + call cholesky(A, L, lower=.true., err=state) if (state%error()) error stop state%print() - print '("Solution: ",*(f8.4,1x))', x + ! First right-hand side + b1 = [1.0_dp, 2.0_dp, 3.0_dp] + call solve_lower_chol(L, b1, x, err=state) + if (state%error()) error stop state%print() + print '("Solution 1: ",*(f8.4,1x))', x - ! For performance-critical code, use overwrite_a=.true. - ! to avoid internal allocation (but A will be destroyed) - ! call cholesky_solve(A, b, x, lower=.true., overwrite_a=.true., err=state) + ! Second right-hand side (reusing the same factorization) + b2 = [4.0_dp, 5.0_dp, 6.0_dp] + call solve_lower_chol(L, b2, x, err=state) + if (state%error()) error stop state%print() + print '("Solution 2: ",*(f8.4,1x))', x end program example_cholesky_solve diff --git a/example/linalg/example_solve_chol.f90 b/example/linalg/example_solve_chol.f90 index 7d5be06d0..00a418e0a 100644 --- a/example/linalg/example_solve_chol.f90 +++ b/example/linalg/example_solve_chol.f90 @@ -1,9 +1,9 @@ program example_solve_chol use stdlib_linalg_constants, only: dp - use stdlib_linalg, only: cholesky, solve_chol, linalg_state_type + use stdlib_linalg, only: solve_chol, linalg_state_type implicit none - real(dp) :: A(3,3), L(3,3), b(3), x(3) + real(dp) :: A(3,3), b(3), x(3) type(linalg_state_type) :: state ! Symmetric positive definite matrix @@ -14,14 +14,14 @@ program example_solve_chol ! Right-hand side b = [1.0_dp, 2.0_dp, 3.0_dp] - ! Compute Cholesky factorization: A = L * L^T - call cholesky(A, L, lower=.true., err=state) - if (state%error()) error stop state%print() - - ! Solve using pre-computed Cholesky factors - call solve_chol(L, b, x, lower=.true., err=state) + ! One-shot Cholesky factorization and solve (A is preserved by default) + call solve_chol(A, b, x, lower=.true., err=state) if (state%error()) error stop state%print() print '("Solution: ",*(f8.4,1x))', x + ! For performance-critical code, use overwrite_a=.true. + ! to avoid internal allocation (but A will be destroyed) + ! call solve_chol(A, b, x, lower=.true., overwrite_a=.true., err=state) + end program example_solve_chol diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index dfe376f35..58b633883 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -46,7 +46,8 @@ module stdlib_linalg public :: solve public :: solve_lu public :: solve_chol - public :: cholesky_solve + public :: solve_lower_chol + public :: solve_upper_chol public :: solve_lstsq public :: solve_constrained_lstsq public :: trace @@ -463,39 +464,42 @@ module stdlib_linalg #:endfor end interface solve_lu - ! Solve linear system Ax = b using pre-computed Cholesky decomposition (subroutine interface) + ! One-shot Cholesky factorization and solve (uses POSV) interface solve_chol !! version: experimental !! !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a - !! symmetric positive definite matrix \( A \) that has been pre-factorized using Cholesky. - !! ([Specification](../page/specs/stdlib_linalg.html#solve_chol-solve-a-linear-system-using-cholesky-factors)) + !! symmetric positive definite matrix \( A \). Combines factorization and solve in one call. + !! ([Specification](../page/specs/stdlib_linalg.html#solve_chol-solve-spd-system-with-cholesky-factorization)) !! !!### Summary - !! Subroutine interface for solving a linear system using pre-computed Cholesky factors. + !! One-shot factorization and solve for SPD systems (wraps LAPACK POSV). !! !!### Description !! - !! This interface provides methods for computing the solution of a linear matrix system using - !! Cholesky factors. Supported data types include `real` and `complex`. Preallocated space - !! for the solution vector `x` is user-provided. The `lower` argument is REQUIRED and must - !! match the `lower` used during the Cholesky factorization. - !! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`) - !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems. + !! This interface computes both the Cholesky factorization and solves the linear system + !! in a single call. Use this for one-time solves. For repeated solves with the same + !! matrix but different RHS, use `cholesky` + `solve_lower_chol`/`solve_upper_chol` for + !! better performance. + !! Supported data types include `real` and `complex`. + !! By default, A is not overwritten. Set `overwrite_a=.true.` to allow in-place + !! factorization for better performance. !! - !!@note The solution is based on LAPACK's `*POTRS` routines. + !!@note The solution is based on LAPACK's `*POSV` routines. !! #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - pure module subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$(a,b,x,lower,err) - !> Input matrix a[n,n] containing Cholesky factors from cholesky - ${rt}$, intent(in) :: a(:,:) + pure module subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$(a,b,x,lower,overwrite_a,err) + !> Input SPD matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, intent(inout), contiguous, target :: x${nd}$ - !> Is the lower triangular factor stored? (REQUIRED) - logical(lk), intent(in) :: lower + !> [optional] Use lower triangular factorization? Default = .true. + logical(lk), optional, intent(in) :: lower + !> [optional] Can A data be overwritten and destroyed? Default = .false. + 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 end subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$ @@ -503,47 +507,77 @@ module stdlib_linalg #:endfor end interface solve_chol - ! One-shot Cholesky factorization and solve (convenience wrapper using POSV) - interface cholesky_solve + ! Solve linear system using pre-computed LOWER Cholesky factor (subroutine interface) + interface solve_lower_chol !! version: experimental !! !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a - !! symmetric positive definite matrix \( A \). Combines factorization and solve in one call. - !! ([Specification](../page/specs/stdlib_linalg.html#cholesky_solve-one-shot-cholesky-solve)) + !! symmetric positive definite matrix \( A \) using pre-computed LOWER Cholesky factor \( L \). + !! ([Specification](../page/specs/stdlib_linalg.html#solve_lower_chol-solve-using-lower-cholesky-factor)) !! !!### Summary - !! One-shot factorization and solve for SPD systems (wraps LAPACK POSV). + !! Subroutine interface for solving a linear system using pre-computed lower Cholesky factor. !! !!### Description !! - !! This interface computes both the Cholesky factorization and solves the linear system - !! in a single call. Use this for one-time solves. For repeated solves with the same - !! matrix but different RHS, use `cholesky` + `solve_chol` for better performance. + !! This interface solves a linear system using a pre-computed lower triangular Cholesky + !! factor \( L \) where \( A = L \cdot L^T \). The input matrix must come from a prior + !! call to `cholesky` with `lower=.true.`. !! Supported data types include `real` and `complex`. - !! By default, A is not overwritten. Set `overwrite_a=.true.` to allow in-place - !! factorization for better performance. !! - !!@note The solution is based on LAPACK's `*POSV` routines. + !!@note The solution is based on LAPACK's `*POTRS` routines. !! #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - pure module subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$(a,b,x,lower,overwrite_a,err) - !> Input SPD matrix a[n,n] - ${rt}$, intent(inout), target :: a(:,:) + pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(a,b,x,err) + !> Input matrix a[n,n] containing lower Cholesky factor L from cholesky(...,lower=.true.) + ${rt}$, intent(in) :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [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_lower_chol_${ndsuf}$ + #:endfor + #:endfor + end interface solve_lower_chol + + ! Solve linear system using pre-computed UPPER Cholesky factor (subroutine interface) + interface solve_upper_chol + !! version: experimental + !! + !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a + !! symmetric positive definite matrix \( A \) using pre-computed UPPER Cholesky factor \( U \). + !! ([Specification](../page/specs/stdlib_linalg.html#solve_upper_chol-solve-using-upper-cholesky-factor)) + !! + !!### Summary + !! Subroutine interface for solving a linear system using pre-computed upper Cholesky factor. + !! + !!### Description + !! + !! This interface solves a linear system using a pre-computed upper triangular Cholesky + !! factor \( U \) where \( A = U^T \cdot U \). The input matrix must come from a prior + !! call to `cholesky` with `lower=.false.`. + !! Supported data types include `real` and `complex`. + !! + !!@note The solution is based on LAPACK's `*POTRS` routines. + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + pure module subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$(a,b,x,err) + !> Input matrix a[n,n] containing upper Cholesky factor U from cholesky(...,lower=.false.) + ${rt}$, intent(in) :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, intent(inout), contiguous, target :: x${nd}$ - !> [optional] Use lower triangular factorization? Default = .true. - logical(lk), optional, intent(in) :: lower - !> [optional] Can A data be overwritten and destroyed? Default = .false. - 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 - end subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$ + end subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$ #:endfor #:endfor - end interface cholesky_solve + end interface solve_upper_chol ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. interface lstsq diff --git a/src/linalg/stdlib_linalg_solve.fypp b/src/linalg/stdlib_linalg_solve.fypp index 76a6baa70..b2b4174ac 100644 --- a/src/linalg/stdlib_linalg_solve.fypp +++ b/src/linalg/stdlib_linalg_solve.fypp @@ -144,29 +144,32 @@ submodule (stdlib_linalg) stdlib_linalg_solve #:endfor !--------------------------------------------------------------------------- - !> solve_chol: Solve using PRE-COMPUTED Cholesky factors (POTRS) + !> solve_chol: One-shot factorize + solve (POSV) !--------------------------------------------------------------------------- #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - !> Solve the linear system A*x = b using pre-computed Cholesky factorization - pure module subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$(a,b,x,lower,err) - !> Cholesky-factorized matrix a[n,n] - ${rt}$, intent(in) :: a(:,:) + !> Factorize and solve A*x = b in one call (uses LAPACK POSV) + pure module subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$(a,b,x,lower,overwrite_a,err) + !> Input SPD matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, intent(inout), contiguous, target :: x${nd}$ - !> Is the lower triangular factor stored? (REQUIRED) - logical(lk), intent(in) :: lower + !> [optional] Use lower triangular factorization? Default = .true. + logical(lk), optional, intent(in) :: lower + !> [optional] Can A data be overwritten and destroyed? Default = .false. + 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 ! Local variables type(linalg_state_type) :: err0 integer(ilp) :: lda,n,ldb,ldx,nrhs,nrhsx,info + logical(lk) :: lower_,copy_a character :: uplo - ${rt}$, pointer :: xmat(:,:) + ${rt}$, pointer :: xmat(:,:),amat(:,:) ! Problem sizes lda = size(a,1,kind=ilp) @@ -176,8 +179,12 @@ submodule (stdlib_linalg) stdlib_linalg_solve ldx = size(x,1,kind=ilp) nrhsx = size(x,kind=ilp)/ldx - ! Set uplo based on lower flag (REQUIRED argument) - uplo = merge('L','U',lower) + ! Default: use lower triangular + lower_ = optval(lower, .true._lk) + uplo = merge('L','U',lower_) + + ! Can A be overwritten? By default, do not overwrite + copy_a = .not. optval(overwrite_a, .false._lk) ! Validate dimensions if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then @@ -187,17 +194,26 @@ submodule (stdlib_linalg) stdlib_linalg_solve return end if - ! Copy RHS to solution array (POTRS overwrites with solution) + ! Initialize a matrix temporary + if (copy_a) then + allocate(amat(lda,n),source=a) + else + amat => a + endif + + ! Copy RHS to solution array (POSV overwrites with solution) x = b ! Create 2D pointer for LAPACK call xmat(1:n,1:nrhs) => x - ! Solve the system using LAPACK POTRS - call potrs(uplo,n,nrhs,a,lda,xmat,n,info) + ! Factorize AND solve using LAPACK POSV + call posv(uplo,n,nrhs,amat,lda,xmat,n,info) ! Handle errors using standard handler - call handle_potrs_info(this,info,uplo,n,nrhs,lda,n,err0) + call handle_posv_info(this,info,uplo,n,nrhs,lda,n,err0) + + if (copy_a) deallocate(amat) ! Process output and return call linalg_error_handling(err0,err) @@ -208,32 +224,27 @@ submodule (stdlib_linalg) stdlib_linalg_solve #:endfor !--------------------------------------------------------------------------- - !> cholesky_solve: One-shot factorize + solve (POSV) + !> solve_lower_chol: Solve using PRE-COMPUTED LOWER Cholesky factor (POTRS) !--------------------------------------------------------------------------- #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - !> Factorize and solve A*x = b in one call (uses LAPACK POSV) - pure module subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$(a,b,x,lower,overwrite_a,err) - !> Input SPD matrix a[n,n] - ${rt}$, intent(inout), target :: a(:,:) + !> Solve the linear system A*x = b using pre-computed lower Cholesky factor + pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(a,b,x,err) + !> Lower Cholesky factor L[n,n] from cholesky(...,lower=.true.) + ${rt}$, intent(in) :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, intent(inout), contiguous, target :: x${nd}$ - !> [optional] Use lower triangular factorization? Default = .true. - logical(lk), optional, intent(in) :: lower - !> [optional] Can A data be overwritten and destroyed? Default = .false. - 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 ! Local variables type(linalg_state_type) :: err0 integer(ilp) :: lda,n,ldb,ldx,nrhs,nrhsx,info - logical(lk) :: lower_,copy_a - character :: uplo - ${rt}$, pointer :: xmat(:,:),amat(:,:) + character, parameter :: uplo = 'L' + ${rt}$, pointer :: xmat(:,:) ! Problem sizes lda = size(a,1,kind=ilp) @@ -243,12 +254,64 @@ submodule (stdlib_linalg) stdlib_linalg_solve ldx = size(x,1,kind=ilp) nrhsx = size(x,kind=ilp)/ldx - ! Default: use lower triangular - lower_ = optval(lower, .true._lk) - uplo = merge('L','U',lower_) + ! Validate dimensions + if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + 'b=',[ldb,nrhs],' x=',[ldx,nrhsx]) + call linalg_error_handling(err0,err) + return + end if - ! Can A be overwritten? By default, do not overwrite - copy_a = .not. optval(overwrite_a, .false._lk) + ! Copy RHS to solution array (POTRS overwrites with solution) + x = b + + ! Create 2D pointer for LAPACK call + xmat(1:n,1:nrhs) => x + + ! Solve the system using LAPACK POTRS with lower triangular factor + call potrs(uplo,n,nrhs,a,lda,xmat,n,info) + + ! Handle errors using standard handler + call handle_potrs_info(this,info,uplo,n,nrhs,lda,n,err0) + + ! Process output and return + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$ + + #:endfor + #:endfor + + !--------------------------------------------------------------------------- + !> solve_upper_chol: Solve using PRE-COMPUTED UPPER Cholesky factor (POTRS) + !--------------------------------------------------------------------------- + + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + !> Solve the linear system A*x = b using pre-computed upper Cholesky factor + pure module subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$(a,b,x,err) + !> Upper Cholesky factor U[n,n] from cholesky(...,lower=.false.) + ${rt}$, intent(in) :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [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) :: lda,n,ldb,ldx,nrhs,nrhsx,info + character, parameter :: uplo = 'U' + ${rt}$, pointer :: xmat(:,:) + + ! Problem sizes + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + ldb = size(b,1,kind=ilp) + nrhs = size(b,kind=ilp)/ldb + ldx = size(x,1,kind=ilp) + nrhsx = size(x,kind=ilp)/ldx ! Validate dimensions if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then @@ -258,31 +321,22 @@ submodule (stdlib_linalg) stdlib_linalg_solve return end if - ! Initialize a matrix temporary - if (copy_a) then - allocate(amat(lda,n),source=a) - else - amat => a - endif - - ! Copy RHS to solution array (POSV overwrites with solution) + ! Copy RHS to solution array (POTRS overwrites with solution) x = b ! Create 2D pointer for LAPACK call xmat(1:n,1:nrhs) => x - ! Factorize AND solve using LAPACK POSV - call posv(uplo,n,nrhs,amat,lda,xmat,n,info) + ! Solve the system using LAPACK POTRS with upper triangular factor + call potrs(uplo,n,nrhs,a,lda,xmat,n,info) ! Handle errors using standard handler - call handle_posv_info(this,info,uplo,n,nrhs,lda,n,err0) - - if (copy_a) deallocate(amat) + call handle_potrs_info(this,info,uplo,n,nrhs,lda,n,err0) ! Process output and return call linalg_error_handling(err0,err) - end subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$ + end subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$ #:endfor #:endfor diff --git a/test/linalg/test_linalg_solve_chol.fypp b/test/linalg/test_linalg_solve_chol.fypp index bbdd17bda..a4719935a 100644 --- a/test/linalg/test_linalg_solve_chol.fypp +++ b/test/linalg/test_linalg_solve_chol.fypp @@ -1,10 +1,10 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -! Test solve_chol and cholesky_solve +! Test solve_chol, solve_lower_chol, solve_upper_chol module test_linalg_solve_chol use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg_constants - use stdlib_linalg, only: cholesky, solve_chol, cholesky_solve + use stdlib_linalg, only: cholesky, solve_chol, solve_lower_chol, solve_upper_chol use stdlib_linalg_state, only: linalg_state_type implicit none (type,external) @@ -22,21 +22,21 @@ module test_linalg_solve_chol allocate(tests(0)) #:for rk,rt,ri in RC_KINDS_TYPES + call add_test(tests,new_unittest("solve_lower_chol_${ri}$",test_solve_lower_chol_${ri}$)) + call add_test(tests,new_unittest("solve_upper_chol_${ri}$",test_solve_upper_chol_${ri}$)) call add_test(tests,new_unittest("solve_chol_${ri}$",test_solve_chol_${ri}$)) call add_test(tests,new_unittest("solve_chol_upper_${ri}$",test_solve_chol_upper_${ri}$)) - call add_test(tests,new_unittest("cholesky_solve_${ri}$",test_cholesky_solve_${ri}$)) - call add_test(tests,new_unittest("cholesky_solve_upper_${ri}$",test_cholesky_solve_upper_${ri}$)) - call add_test(tests,new_unittest("cholesky_solve_overwrite_${ri}$",test_cholesky_solve_overwrite_${ri}$)) - call add_test(tests,new_unittest("solve_chol_multi_rhs_${ri}$",test_solve_chol_multi_rhs_${ri}$)) - call add_test(tests,new_unittest("cholesky_solve_indefinite_${ri}$",test_cholesky_solve_indefinite_${ri}$)) - call add_test(tests,new_unittest("cholesky_solve_semidefinite_${ri}$",test_cholesky_solve_semidefinite_${ri}$)) + call add_test(tests,new_unittest("solve_chol_overwrite_${ri}$",test_solve_chol_overwrite_${ri}$)) + call add_test(tests,new_unittest("solve_lower_chol_multi_rhs_${ri}$",test_solve_lower_chol_multi_rhs_${ri}$)) + call add_test(tests,new_unittest("solve_chol_indefinite_${ri}$",test_solve_chol_indefinite_${ri}$)) + call add_test(tests,new_unittest("solve_chol_semidefinite_${ri}$",test_solve_chol_semidefinite_${ri}$)) #:endfor end subroutine test_solve_chol_factorization - !> Test solve_chol with pre-computed Cholesky factors + !> Test solve_lower_chol with pre-computed lower Cholesky factors #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_solve_chol_${ri}$(error) + subroutine test_solve_lower_chol_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 3_ilp @@ -61,23 +61,23 @@ module test_linalg_solve_chol call check(error, state%ok(), 'cholesky factorization failed: '//state%print()) if (allocated(error)) return - ! Solve using Cholesky factors - call solve_chol(l, b, x, lower=.true., err=state) + ! Solve using lower Cholesky factors + call solve_lower_chol(l, b, x, err=state) - call check(error, state%ok(), 'solve_chol failed: '//state%print()) + call check(error, state%ok(), 'solve_lower_chol failed: '//state%print()) if (allocated(error)) return ! Check solution - call check(error, all(abs(x - x_expected) < tol), 'solve_chol: solution mismatch') + call check(error, all(abs(x - x_expected) < tol), 'solve_lower_chol: solution mismatch') if (allocated(error)) return - end subroutine test_solve_chol_${ri}$ + end subroutine test_solve_lower_chol_${ri}$ #:endfor - !> Test solve_chol with upper triangular factors (lower=.false.) + !> Test solve_upper_chol with pre-computed upper Cholesky factors #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_solve_chol_upper_${ri}$(error) + subroutine test_solve_upper_chol_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 3_ilp @@ -103,22 +103,22 @@ module test_linalg_solve_chol if (allocated(error)) return ! Solve using upper Cholesky factors - call solve_chol(u, b, x, lower=.false., err=state) + call solve_upper_chol(u, b, x, err=state) - call check(error, state%ok(), 'solve_chol (upper) failed: '//state%print()) + call check(error, state%ok(), 'solve_upper_chol failed: '//state%print()) if (allocated(error)) return ! Check solution - call check(error, all(abs(x - x_expected) < tol), 'solve_chol (upper): solution mismatch') + call check(error, all(abs(x - x_expected) < tol), 'solve_upper_chol: solution mismatch') if (allocated(error)) return - end subroutine test_solve_chol_upper_${ri}$ + end subroutine test_solve_upper_chol_${ri}$ #:endfor - !> Test cholesky_solve (one-shot) - default preserves A + !> Test solve_chol (one-shot) - default preserves A #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_cholesky_solve_${ri}$(error) + subroutine test_solve_chol_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 3_ilp @@ -139,26 +139,26 @@ module test_linalg_solve_chol b = matmul(a, x_expected) ! One-shot solve with default overwrite_a=.false. (A should be preserved) - call cholesky_solve(a, b, x, lower=.true., err=state) + call solve_chol(a, b, x, lower=.true., err=state) - call check(error, state%ok(), 'cholesky_solve failed: '//state%print()) + call check(error, state%ok(), 'solve_chol failed: '//state%print()) if (allocated(error)) return ! Check solution - call check(error, all(abs(x - x_expected) < tol), 'cholesky_solve: solution mismatch') + call check(error, all(abs(x - x_expected) < tol), 'solve_chol: solution mismatch') if (allocated(error)) return ! Check that A is preserved (default behavior) - call check(error, all(abs(a - a_copy) < tol), 'cholesky_solve: A should be preserved by default') + call check(error, all(abs(a - a_copy) < tol), 'solve_chol: A should be preserved by default') if (allocated(error)) return - end subroutine test_cholesky_solve_${ri}$ + end subroutine test_solve_chol_${ri}$ #:endfor - !> Test cholesky_solve with upper triangular (lower=.false.) + !> Test solve_chol with upper triangular (lower=.false.) #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_cholesky_solve_upper_${ri}$(error) + subroutine test_solve_chol_upper_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 3_ilp @@ -179,26 +179,26 @@ module test_linalg_solve_chol b = matmul(a, x_expected) ! One-shot solve with upper triangular (A should be preserved by default) - call cholesky_solve(a, b, x, lower=.false., err=state) + call solve_chol(a, b, x, lower=.false., err=state) - call check(error, state%ok(), 'cholesky_solve (upper) failed: '//state%print()) + call check(error, state%ok(), 'solve_chol (upper) failed: '//state%print()) if (allocated(error)) return ! Check solution - call check(error, all(abs(x - x_expected) < tol), 'cholesky_solve (upper): solution mismatch') + call check(error, all(abs(x - x_expected) < tol), 'solve_chol (upper): solution mismatch') if (allocated(error)) return ! Check that A is preserved (default behavior) - call check(error, all(abs(a - a_copy) < tol), 'cholesky_solve (upper): A should be preserved') + call check(error, all(abs(a - a_copy) < tol), 'solve_chol (upper): A should be preserved') if (allocated(error)) return - end subroutine test_cholesky_solve_upper_${ri}$ + end subroutine test_solve_chol_upper_${ri}$ #:endfor - !> Test cholesky_solve with overwrite_a=.true. + !> Test solve_chol with overwrite_a=.true. #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_cholesky_solve_overwrite_${ri}$(error) + subroutine test_solve_chol_overwrite_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 3_ilp @@ -219,26 +219,26 @@ module test_linalg_solve_chol b = matmul(a, x_expected) ! One-shot solve with overwrite_a=.true. (A will be destroyed) - call cholesky_solve(a, b, x, lower=.true., overwrite_a=.true., err=state) + call solve_chol(a, b, x, lower=.true., overwrite_a=.true., err=state) - call check(error, state%ok(), 'cholesky_solve overwrite failed: '//state%print()) + call check(error, state%ok(), 'solve_chol overwrite failed: '//state%print()) if (allocated(error)) return ! Check solution - call check(error, all(abs(x - x_expected) < tol), 'cholesky_solve overwrite: solution mismatch') + call check(error, all(abs(x - x_expected) < tol), 'solve_chol overwrite: solution mismatch') if (allocated(error)) return ! Check that A was overwritten (not equal to original) - call check(error, any(abs(a - a_copy) > tol), 'cholesky_solve: A should be overwritten with overwrite_a=.true.') + call check(error, any(abs(a - a_copy) > tol), 'solve_chol: A should be overwritten with overwrite_a=.true.') if (allocated(error)) return - end subroutine test_cholesky_solve_overwrite_${ri}$ + end subroutine test_solve_chol_overwrite_${ri}$ #:endfor - !> Test solve_chol with multiple RHS + !> Test solve_lower_chol with multiple RHS #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_solve_chol_multi_rhs_${ri}$(error) + subroutine test_solve_lower_chol_multi_rhs_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 3_ilp @@ -265,23 +265,23 @@ module test_linalg_solve_chol call check(error, state%ok(), 'cholesky factorization failed: '//state%print()) if (allocated(error)) return - ! Solve using Cholesky factors - call solve_chol(l, b, x, lower=.true., err=state) + ! Solve using lower Cholesky factors + call solve_lower_chol(l, b, x, err=state) - call check(error, state%ok(), 'solve_chol multi-rhs failed: '//state%print()) + call check(error, state%ok(), 'solve_lower_chol multi-rhs failed: '//state%print()) if (allocated(error)) return ! Check solution - call check(error, all(abs(x - x_expected) < tol), 'solve_chol multi-rhs: solution mismatch') + call check(error, all(abs(x - x_expected) < tol), 'solve_lower_chol multi-rhs: solution mismatch') if (allocated(error)) return - end subroutine test_solve_chol_multi_rhs_${ri}$ + end subroutine test_solve_lower_chol_multi_rhs_${ri}$ #:endfor - !> Test cholesky_solve with symmetric indefinite matrix (should fail) + !> Test solve_chol with symmetric indefinite matrix (should fail) #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_cholesky_solve_indefinite_${ri}$(error) + subroutine test_solve_chol_indefinite_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 2_ilp @@ -295,20 +295,20 @@ module test_linalg_solve_chol ! Arbitrary RHS b = [1, 1] - ! cholesky_solve should fail for indefinite matrix - call cholesky_solve(a, b, x, lower=.true., err=state) + ! solve_chol should fail for indefinite matrix + call solve_chol(a, b, x, lower=.true., err=state) ! Check that it failed (not positive definite) - call check(error, state%error(), 'cholesky_solve should fail for indefinite matrix') + call check(error, state%error(), 'solve_chol should fail for indefinite matrix') if (allocated(error)) return - end subroutine test_cholesky_solve_indefinite_${ri}$ + end subroutine test_solve_chol_indefinite_${ri}$ #:endfor - !> Test cholesky_solve with symmetric positive semi-definite matrix (should fail) + !> Test solve_chol with symmetric positive semi-definite matrix (should fail) #:for rk,rt,ri in RC_KINDS_TYPES - subroutine test_cholesky_solve_semidefinite_${ri}$(error) + subroutine test_solve_chol_semidefinite_${ri}$(error) type(error_type), allocatable, intent(out) :: error integer(ilp), parameter :: n = 2_ilp @@ -322,14 +322,14 @@ module test_linalg_solve_chol ! Arbitrary RHS b = [1, 1] - ! cholesky_solve should fail for semi-definite (singular) matrix - call cholesky_solve(a, b, x, lower=.true., err=state) + ! solve_chol should fail for semi-definite (singular) matrix + call solve_chol(a, b, x, lower=.true., err=state) ! Check that it failed (not strictly positive definite) - call check(error, state%error(), 'cholesky_solve should fail for positive semi-definite matrix') + call check(error, state%error(), 'solve_chol should fail for positive semi-definite matrix') if (allocated(error)) return - end subroutine test_cholesky_solve_semidefinite_${ri}$ + end subroutine test_solve_chol_semidefinite_${ri}$ #:endfor From da8c1df71c15fc41b1035bafff6f0a3cc08bbbf9 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Fri, 30 Jan 2026 20:17:52 +0000 Subject: [PATCH 7/8] Refactor: DRY driver for solve_lower/upper_chol, rename examples to match API --- doc/specs/stdlib_linalg.md | 8 +- example/linalg/CMakeLists.txt | 3 +- ...solve.f90 => example_solve_lower_chol.f90} | 9 +-- example/linalg/example_solve_upper_chol.f90 | 30 ++++++++ src/linalg/stdlib_linalg_solve.fypp | 75 ++++++++----------- 5 files changed, 75 insertions(+), 50 deletions(-) rename example/linalg/{example_cholesky_solve.f90 => example_solve_lower_chol.f90} (73%) create mode 100644 example/linalg/example_solve_upper_chol.f90 diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 6ec1a1770..728b39ed4 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -827,7 +827,7 @@ If `err` is not present, exceptions trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_cholesky_solve.f90!} +{!example/linalg/example_solve_lower_chol.f90!} ``` ## `solve_upper_chol` - Solves a linear system using pre-computed upper Cholesky factor. @@ -865,6 +865,12 @@ For a correctly factorized matrix, returns an array value that represents the so Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. +### Example + +```fortran +{!example/linalg/example_solve_upper_chol.f90!} +``` + ## `lstsq` - Computes the least squares solution to a linear matrix equation. ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1e313cf6f..76dc35963 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -64,6 +64,7 @@ ADD_EXAMPLE(pivoting_qr_space) ADD_EXAMPLE(cholesky) ADD_EXAMPLE(chol) ADD_EXAMPLE(solve_chol) -ADD_EXAMPLE(cholesky_solve) +ADD_EXAMPLE(solve_lower_chol) +ADD_EXAMPLE(solve_upper_chol) ADD_EXAMPLE(expm) ADD_EXAMPLE(matrix_exp) diff --git a/example/linalg/example_cholesky_solve.f90 b/example/linalg/example_solve_lower_chol.f90 similarity index 73% rename from example/linalg/example_cholesky_solve.f90 rename to example/linalg/example_solve_lower_chol.f90 index 2b4fdf952..98c2856a6 100644 --- a/example/linalg/example_cholesky_solve.f90 +++ b/example/linalg/example_solve_lower_chol.f90 @@ -1,7 +1,4 @@ -! Example: solve_lower_chol - Solve using pre-computed Cholesky factors -! For repeated solves with the same matrix, pre-compute the factorization -! once and reuse it for better performance. -program example_cholesky_solve +program example_solve_lower_chol use stdlib_linalg_constants, only: dp use stdlib_linalg, only: cholesky, solve_lower_chol, linalg_state_type implicit none @@ -14,7 +11,7 @@ program example_cholesky_solve A(2,:) = [2.0_dp, 5.0_dp, 1.0_dp] A(3,:) = [2.0_dp, 1.0_dp, 6.0_dp] - ! Compute Cholesky factorization once: A = L * L^T + ! Compute lower Cholesky factorization once: A = L * L^T call cholesky(A, L, lower=.true., err=state) if (state%error()) error stop state%print() @@ -30,4 +27,4 @@ program example_cholesky_solve if (state%error()) error stop state%print() print '("Solution 2: ",*(f8.4,1x))', x -end program example_cholesky_solve +end program example_solve_lower_chol diff --git a/example/linalg/example_solve_upper_chol.f90 b/example/linalg/example_solve_upper_chol.f90 new file mode 100644 index 000000000..7c3935041 --- /dev/null +++ b/example/linalg/example_solve_upper_chol.f90 @@ -0,0 +1,30 @@ +program example_solve_upper_chol + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: cholesky, solve_upper_chol, linalg_state_type + implicit none + + real(dp) :: A(3,3), U(3,3), b1(3), b2(3), x(3) + type(linalg_state_type) :: state + + ! Symmetric positive definite matrix + A(1,:) = [4.0_dp, 2.0_dp, 2.0_dp] + A(2,:) = [2.0_dp, 5.0_dp, 1.0_dp] + A(3,:) = [2.0_dp, 1.0_dp, 6.0_dp] + + ! Compute upper Cholesky factorization once: A = U^T * U + call cholesky(A, U, lower=.false., err=state) + if (state%error()) error stop state%print() + + ! First right-hand side + b1 = [1.0_dp, 2.0_dp, 3.0_dp] + call solve_upper_chol(U, b1, x, err=state) + if (state%error()) error stop state%print() + print '("Solution 1: ",*(f8.4,1x))', x + + ! Second right-hand side (reusing the same factorization) + b2 = [4.0_dp, 5.0_dp, 6.0_dp] + call solve_upper_chol(U, b2, x, err=state) + if (state%error()) error stop state%print() + print '("Solution 2: ",*(f8.4,1x))', x + +end program example_solve_upper_chol diff --git a/src/linalg/stdlib_linalg_solve.fypp b/src/linalg/stdlib_linalg_solve.fypp index b2b4174ac..ee9fe3d95 100644 --- a/src/linalg/stdlib_linalg_solve.fypp +++ b/src/linalg/stdlib_linalg_solve.fypp @@ -224,26 +224,28 @@ submodule (stdlib_linalg) stdlib_linalg_solve #:endfor !--------------------------------------------------------------------------- - !> solve_lower_chol: Solve using PRE-COMPUTED LOWER Cholesky factor (POTRS) + !> Private driver: Solve using pre-computed Cholesky factor (POTRS) + !> Not exported - used internally by solve_lower_chol and solve_upper_chol !--------------------------------------------------------------------------- #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - !> Solve the linear system A*x = b using pre-computed lower Cholesky factor - pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(a,b,x,err) - !> Lower Cholesky factor L[n,n] from cholesky(...,lower=.true.) + !> Low-level driver for solving A*x = b using pre-computed Cholesky factor + pure subroutine solve_chol_${ri}$_${ndsuf}$_driver(a,b,x,uplo,err) + !> Cholesky factor (L or U)[n,n] from cholesky(...) ${rt}$, intent(in) :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> Triangle selector: 'L' for lower, 'U' for upper + character, intent(in) :: uplo !> [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) :: lda,n,ldb,ldx,nrhs,nrhsx,info - character, parameter :: uplo = 'L' ${rt}$, pointer :: xmat(:,:) ! Problem sizes @@ -268,7 +270,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve ! Create 2D pointer for LAPACK call xmat(1:n,1:nrhs) => x - ! Solve the system using LAPACK POTRS with lower triangular factor + ! Solve the system using LAPACK POTRS call potrs(uplo,n,nrhs,a,lda,xmat,n,info) ! Handle errors using standard handler @@ -277,6 +279,30 @@ submodule (stdlib_linalg) stdlib_linalg_solve ! Process output and return call linalg_error_handling(err0,err) + end subroutine solve_chol_${ri}$_${ndsuf}$_driver + + #:endfor + #:endfor + + !--------------------------------------------------------------------------- + !> solve_lower_chol: Solve using PRE-COMPUTED LOWER Cholesky factor (POTRS) + !--------------------------------------------------------------------------- + + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + !> Solve the linear system A*x = b using pre-computed lower Cholesky factor + pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(a,b,x,err) + !> Lower Cholesky factor L[n,n] from cholesky(...,lower=.true.) + ${rt}$, intent(in) :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] State return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + call solve_chol_${ri}$_${ndsuf}$_driver(a,b,x,'L',err) + end subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$ #:endfor @@ -299,42 +325,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve !> [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) :: lda,n,ldb,ldx,nrhs,nrhsx,info - character, parameter :: uplo = 'U' - ${rt}$, pointer :: xmat(:,:) - - ! Problem sizes - lda = size(a,1,kind=ilp) - n = size(a,2,kind=ilp) - ldb = size(b,1,kind=ilp) - nrhs = size(b,kind=ilp)/ldb - ldx = size(x,1,kind=ilp) - nrhsx = size(x,kind=ilp)/ldx - - ! Validate dimensions - if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & - 'b=',[ldb,nrhs],' x=',[ldx,nrhsx]) - call linalg_error_handling(err0,err) - return - end if - - ! Copy RHS to solution array (POTRS overwrites with solution) - x = b - - ! Create 2D pointer for LAPACK call - xmat(1:n,1:nrhs) => x - - ! Solve the system using LAPACK POTRS with upper triangular factor - call potrs(uplo,n,nrhs,a,lda,xmat,n,info) - - ! Handle errors using standard handler - call handle_potrs_info(this,info,uplo,n,nrhs,lda,n,err0) - - ! Process output and return - call linalg_error_handling(err0,err) + call solve_chol_${ri}$_${ndsuf}$_driver(a,b,x,'U',err) end subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$ From f429a98cada96c1a8d9458cc49d4fb7f4854c729 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Mon, 2 Feb 2026 19:45:59 +0000 Subject: [PATCH 8/8] Rename Cholesky solver params: a -> l/u for clarity --- doc/specs/stdlib_linalg.md | 12 ++++++------ src/linalg/stdlib_linalg.fypp | 12 ++++++------ src/linalg/stdlib_linalg_solve.fypp | 16 ++++++++-------- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 728b39ed4..be970de9f 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -805,13 +805,13 @@ The solver is based on LAPACK's `*POTRS` backends. ### Syntax -`call ` [[stdlib_linalg(module):solve_lower_chol(interface)]] `(a, b, x [, err])` +`call ` [[stdlib_linalg(module):solve_lower_chol(interface)]] `(l, b, x [, err])` ### Arguments -`a`: Shall be a rank-2 `real` or `complex` square array containing the **lower** Cholesky factor `L` (output of `cholesky(..., lower=.true.)`). It is an `intent(in)` argument. +`l`: Shall be a rank-2 `real` or `complex` square array containing the **lower** Cholesky factor `L` (output of `cholesky(..., lower=.true.)`). It is an `intent(in)` argument. -`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. +`b`: Shall be a rank-1 or rank-2 array of the same kind as `l`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. @@ -846,13 +846,13 @@ The solver is based on LAPACK's `*POTRS` backends. ### Syntax -`call ` [[stdlib_linalg(module):solve_upper_chol(interface)]] `(a, b, x [, err])` +`call ` [[stdlib_linalg(module):solve_upper_chol(interface)]] `(u, b, x [, err])` ### Arguments -`a`: Shall be a rank-2 `real` or `complex` square array containing the **upper** Cholesky factor `U` (output of `cholesky(..., lower=.false.)`). It is an `intent(in)` argument. +`u`: Shall be a rank-2 `real` or `complex` square array containing the **upper** Cholesky factor `U` (output of `cholesky(..., lower=.false.)`). It is an `intent(in)` argument. -`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. +`b`: Shall be a rank-1 or rank-2 array of the same kind as `u`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp index 58b633883..0af97b84e 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -529,9 +529,9 @@ module stdlib_linalg !! #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(a,b,x,err) - !> Input matrix a[n,n] containing lower Cholesky factor L from cholesky(...,lower=.true.) - ${rt}$, intent(in) :: a(:,:) + pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(l,b,x,err) + !> Input matrix l[n,n] containing lower Cholesky factor L from cholesky(...,lower=.true.) + ${rt}$, intent(in) :: l(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] @@ -565,9 +565,9 @@ module stdlib_linalg !! #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES - pure module subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$(a,b,x,err) - !> Input matrix a[n,n] containing upper Cholesky factor U from cholesky(...,lower=.false.) - ${rt}$, intent(in) :: a(:,:) + pure module subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$(u,b,x,err) + !> Input matrix u[n,n] containing upper Cholesky factor U from cholesky(...,lower=.false.) + ${rt}$, intent(in) :: u(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] diff --git a/src/linalg/stdlib_linalg_solve.fypp b/src/linalg/stdlib_linalg_solve.fypp index ee9fe3d95..6cf1f9ed5 100644 --- a/src/linalg/stdlib_linalg_solve.fypp +++ b/src/linalg/stdlib_linalg_solve.fypp @@ -291,9 +291,9 @@ submodule (stdlib_linalg) stdlib_linalg_solve #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES !> Solve the linear system A*x = b using pre-computed lower Cholesky factor - pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(a,b,x,err) - !> Lower Cholesky factor L[n,n] from cholesky(...,lower=.true.) - ${rt}$, intent(in) :: a(:,:) + pure module subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$(l,b,x,err) + !> Lower Cholesky factor l[n,n] from cholesky(...,lower=.true.) + ${rt}$, intent(in) :: l(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] @@ -301,7 +301,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve !> [optional] State return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err - call solve_chol_${ri}$_${ndsuf}$_driver(a,b,x,'L',err) + call solve_chol_${ri}$_${ndsuf}$_driver(l,b,x,'L',err) end subroutine stdlib_linalg_${ri}$_solve_lower_chol_${ndsuf}$ @@ -315,9 +315,9 @@ submodule (stdlib_linalg) stdlib_linalg_solve #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES !> Solve the linear system A*x = b using pre-computed upper Cholesky factor - pure module subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$(a,b,x,err) - !> Upper Cholesky factor U[n,n] from cholesky(...,lower=.false.) - ${rt}$, intent(in) :: a(:,:) + pure module subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$(u,b,x,err) + !> Upper Cholesky factor u[n,n] from cholesky(...,lower=.false.) + ${rt}$, intent(in) :: u(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Result array/matrix x[n] or x[n,nrhs] @@ -325,7 +325,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve !> [optional] State return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err - call solve_chol_${ri}$_${ndsuf}$_driver(a,b,x,'U',err) + call solve_chol_${ri}$_${ndsuf}$_driver(u,b,x,'U',err) end subroutine stdlib_linalg_${ri}$_solve_upper_chol_${ndsuf}$