Solving a linear system with a symmetric positive definite matrix

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 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_cholesky_solve.f90!}
+```
+
 ## `lstsq` - Computes the least squares solution to a linear matrix equation. 
 
 ### Status

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)

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

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

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

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