Solving a linear system with a symmetric positive definite matrix (#1093)

Author: perazz

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

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)

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

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

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

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