feat(linalg): move to solve submodule, add overwrite_a

Author: aamrindersingh

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

example/linalg/CMakeLists.txt

@@ -61,5 +61,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/CMakeLists.txt

@@ -50,7 +50,6 @@ set(fppFiles
     stdlib_linalg_state.fypp
     stdlib_linalg_svd.fypp
     stdlib_linalg_cholesky.fypp
-    stdlib_linalg_solve_chol.fypp
     stdlib_linalg_schur.fypp
     stdlib_linalg_matrix_functions.fypp
     stdlib_optval.fypp

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/stdlib_linalg.fypp

@@ -493,7 +493,7 @@ module stdlib_linalg
          !> 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(out), contiguous, target :: x${nd}$
+         ${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
@@ -520,20 +520,24 @@ module stdlib_linalg
      !! 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,err)     
-         !> Input SPD matrix a[n,n] (overwritten with Cholesky factors)
-         ${rt}$, intent(inout) :: a(:,:)
+     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(out), contiguous, target :: x${nd}$
+         ${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}$