diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 082efaefa..be970de9f 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -736,6 +736,141 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_solve3.f90!} ``` +## `solve_chol` - Solves a linear system using Cholesky factorization (one-shot 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_lower_chol`/`solve_upper_chol` for better performance. +The solver is based on LAPACK's `*POSV` backends. + +### Syntax + +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 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` (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_solve_chol.f90!} +``` + +## `solve_lower_chol` - Solves a linear system using pre-computed lower 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=.true.`). + +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_lower_chol(interface)]] `(l, b, x [, err])` + +### Arguments + +`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 `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. + +`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_lower_chol.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)]] `(u, b, x [, err])` + +### Arguments + +`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 `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. + +`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_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 1c6c7e42c..76dc35963 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -63,5 +63,8 @@ ADD_EXAMPLE(qr_space) ADD_EXAMPLE(pivoting_qr_space) ADD_EXAMPLE(cholesky) ADD_EXAMPLE(chol) +ADD_EXAMPLE(solve_chol) +ADD_EXAMPLE(solve_lower_chol) +ADD_EXAMPLE(solve_upper_chol) ADD_EXAMPLE(expm) ADD_EXAMPLE(matrix_exp) diff --git a/example/linalg/example_solve_chol.f90 b/example/linalg/example_solve_chol.f90 new file mode 100644 index 000000000..00a418e0a --- /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: solve_chol, 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 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/example/linalg/example_solve_lower_chol.f90 b/example/linalg/example_solve_lower_chol.f90 new file mode 100644 index 000000000..98c2856a6 --- /dev/null +++ b/example/linalg/example_solve_lower_chol.f90 @@ -0,0 +1,30 @@ +program example_solve_lower_chol + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: cholesky, solve_lower_chol, linalg_state_type + implicit none + + real(dp) :: A(3,3), L(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 lower Cholesky factorization once: A = L * L^T + call cholesky(A, L, lower=.true., 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_lower_chol(L, 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_lower_chol(L, b2, x, err=state) + if (state%error()) error stop state%print() + print '("Solution 2: ",*(f8.4,1x))', x + +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/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 475424c91..63e6ae0f7 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) + err%state = LINALG_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) + err%state = LINALG_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..0af97b84e 100644 --- a/src/linalg/stdlib_linalg.fypp +++ b/src/linalg/stdlib_linalg.fypp @@ -44,7 +44,10 @@ module stdlib_linalg public :: mnorm public :: get_norm public :: solve - public :: solve_lu + public :: solve_lu + public :: solve_chol + public :: solve_lower_chol + public :: solve_upper_chol public :: solve_lstsq public :: solve_constrained_lstsq public :: trace @@ -457,9 +460,124 @@ 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 + + ! 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 \). Combines factorization and solve in one call. + !! ([Specification](../page/specs/stdlib_linalg.html#solve_chol-solve-spd-system-with-cholesky-factorization)) + !! + !!### 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_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 `*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,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}$_solve_chol_${ndsuf}$ #:endfor #:endfor - end interface solve_lu + end interface solve_chol + + ! 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 \) using pre-computed LOWER Cholesky factor \( L \). + !! ([Specification](../page/specs/stdlib_linalg.html#solve_lower_chol-solve-using-lower-cholesky-factor)) + !! + !!### Summary + !! Subroutine interface for solving a linear system using pre-computed lower Cholesky factor. + !! + !!### Description + !! + !! 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`. + !! + !!@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_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] + ${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}$(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] + ${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_upper_chol_${ndsuf}$ + #:endfor + #:endfor + 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 b4df33133..6cf1f9ed5 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,196 @@ submodule (stdlib_linalg) stdlib_linalg_solve end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$ - #:endfor - #:endfor + #:endfor + #:endfor + + !--------------------------------------------------------------------------- + !> solve_chol: 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}$_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}$ + !> [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(:,:) + + ! 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_ = 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 + 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(uplo,n,nrhs,amat,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) + + ! Process output and return + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$ + + #:endfor + #:endfor + + !--------------------------------------------------------------------------- + !> 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 + !> 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 + ${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 + 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 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}$(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] + ${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(l,b,x,'L',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}$(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] + ${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(u,b,x,'U',err) + + end subroutine stdlib_linalg_${ri}$_solve_upper_chol_${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..a4719935a --- /dev/null +++ b/test/linalg/test_linalg_solve_chol.fypp @@ -0,0 +1,383 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! 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, solve_lower_chol, solve_upper_chol + 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_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("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_lower_chol with pre-computed lower Cholesky factors + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_lower_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 lower Cholesky factors + call solve_lower_chol(l, b, x, err=state) + + 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_lower_chol: solution mismatch') + if (allocated(error)) return + + end subroutine test_solve_lower_chol_${ri}$ + + #:endfor + + !> Test solve_upper_chol with pre-computed upper Cholesky factors + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_upper_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), 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_upper_chol(u, b, x, err=state) + + 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_upper_chol: solution mismatch') + if (allocated(error)) return + + end subroutine test_solve_upper_chol_${ri}$ + + #:endfor + + !> Test solve_chol (one-shot) - default preserves A + #: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), 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 solve_chol(a, 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 + + ! Check that A is preserved (default behavior) + call check(error, all(abs(a - a_copy) < tol), 'solve_chol: A should be preserved by default') + if (allocated(error)) return + + end subroutine test_solve_chol_${ri}$ + + #:endfor + + !> Test solve_chol with upper triangular (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), 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 solve_chol(a, 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 + + ! Check that A is preserved (default behavior) + call check(error, all(abs(a - a_copy) < tol), 'solve_chol (upper): A should be preserved') + if (allocated(error)) return + + end subroutine test_solve_chol_upper_${ri}$ + + #:endfor + + !> Test solve_chol with overwrite_a=.true. + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_chol_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 solve_chol(a, b, x, lower=.true., overwrite_a=.true., err=state) + + 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), '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), 'solve_chol: A should be overwritten with overwrite_a=.true.') + if (allocated(error)) return + + end subroutine test_solve_chol_overwrite_${ri}$ + + #:endfor + + !> Test solve_lower_chol with multiple RHS + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_lower_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 lower Cholesky factors + call solve_lower_chol(l, b, x, err=state) + + 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_lower_chol multi-rhs: solution mismatch') + if (allocated(error)) return + + end subroutine test_solve_lower_chol_multi_rhs_${ri}$ + + #:endfor + + !> Test solve_chol with symmetric indefinite matrix (should fail) + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_chol_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] + + ! 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(), 'solve_chol should fail for indefinite matrix') + if (allocated(error)) return + + end subroutine test_solve_chol_indefinite_${ri}$ + + #:endfor + + !> Test solve_chol with symmetric positive semi-definite matrix (should fail) + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_solve_chol_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] + + ! 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(), 'solve_chol should fail for positive semi-definite matrix') + if (allocated(error)) return + + end subroutine test_solve_chol_semidefinite_${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