Skip to content

Commit bfd9b12

Browse files
aamrindersinghaamrindersingh
committed
fix: rebase onto v0.8.0, resolve merge conflicts
1 parent b32b42a commit bfd9b12

9 files changed

Lines changed: 698 additions & 6 deletions

File tree

‎doc/specs/stdlib_linalg.md‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -736,6 +736,102 @@ If `err` is not present, exceptions trigger an `error stop`.
736736
{!example/linalg/example_solve3.f90!}
737737
```
738738

739+
## `solve_chol` - Solves a linear system using pre-computed Cholesky factors (subroutine interface).
740+
741+
### Status
742+
743+
Experimental
744+
745+
### Description
746+
747+
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`).
748+
749+
Result vector or array `x` returns the exact solution to within numerical precision, provided that the factorization is correct.
750+
An error is returned if the matrix and right-hand-side have incompatible sizes.
751+
The solver is based on LAPACK's `*POTRS` backends.
752+
753+
### Syntax
754+
755+
`call ` [[stdlib_linalg(module):solve_chol(interface)]] `(a, b, x, lower [, err])`
756+
757+
### Arguments
758+
759+
`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.
760+
761+
`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.
762+
763+
`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.
764+
765+
`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.
766+
767+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
768+
769+
### Return value
770+
771+
For a correctly factorized matrix, returns an array value that represents the solution to the linear system of equations.
772+
773+
Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes.
774+
If `err` is not present, exceptions trigger an `error stop`.
775+
776+
### Example
777+
778+
```fortran
779+
{!example/linalg/example_solve_chol.f90!}
780+
```
781+
782+
## `cholesky_solve` - Solves a linear matrix equation using Cholesky factorization (subroutine interface).
783+
784+
### Status
785+
786+
Experimental
787+
788+
### Description
789+
790+
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.
791+
792+
Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is positive definite.
793+
An error is returned if the matrix is not positive definite or has incompatible sizes with the right-hand-side.
794+
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.
795+
The solver is based on LAPACK's `*POSV` backends.
796+
797+
### Syntax
798+
799+
Simple (`Pure`) interface:
800+
801+
`call ` [[stdlib_linalg(module):cholesky_solve(interface)]] `(a, b, x)`
802+
803+
Expert (`Pure`) interface:
804+
805+
`call ` [[stdlib_linalg(module):cholesky_solve(interface)]] `(a, b, x [, lower, overwrite_a, err])`
806+
807+
### Arguments
808+
809+
`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.
810+
811+
`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.
812+
813+
`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.
814+
815+
`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.
816+
817+
`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.
818+
819+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
820+
821+
### Return value
822+
823+
For a positive definite matrix, returns an array value that represents the solution to the linear system of equations.
824+
825+
Raises `LINALG_ERROR` if the matrix is not positive definite.
826+
Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes.
827+
If `err` is not present, exceptions trigger an `error stop`.
828+
829+
### Example
830+
831+
```fortran
832+
{!example/linalg/example_cholesky_solve.f90!}
833+
```
834+
739835
## `lstsq` - Computes the least squares solution to a linear matrix equation.
740836

741837
### Status

‎example/linalg/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,5 +61,7 @@ ADD_EXAMPLE(qr_space)
6161
ADD_EXAMPLE(pivoting_qr_space)
6262
ADD_EXAMPLE(cholesky)
6363
ADD_EXAMPLE(chol)
64+
ADD_EXAMPLE(solve_chol)
65+
ADD_EXAMPLE(cholesky_solve)
6466
ADD_EXAMPLE(expm)
6567
ADD_EXAMPLE(matrix_exp)

‎example/linalg/example_cholesky_solve.f90‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
program example_cholesky_solve
2+
use stdlib_linalg_constants, only: dp
3+
use stdlib_linalg, only: cholesky_solve, linalg_state_type
4+
implicit none
5+
6+
real(dp) :: A(3,3), b(3), x(3)
7+
type(linalg_state_type) :: state
8+
9+
! Symmetric positive definite matrix
10+
A(1,:) = [4.0_dp, 2.0_dp, 2.0_dp]
11+
A(2,:) = [2.0_dp, 5.0_dp, 1.0_dp]
12+
A(3,:) = [2.0_dp, 1.0_dp, 6.0_dp]
13+
14+
! Right-hand side
15+
b = [1.0_dp, 2.0_dp, 3.0_dp]
16+
17+
! One-shot factorization and solve (A is preserved by default)
18+
call cholesky_solve(A, b, x, lower=.true., err=state)
19+
if (state%error()) error stop state%print()
20+
21+
print '("Solution: ",*(f8.4,1x))', x
22+
23+
! For performance-critical code, use overwrite_a=.true.
24+
! to avoid internal allocation (but A will be destroyed)
25+
! call cholesky_solve(A, b, x, lower=.true., overwrite_a=.true., err=state)
26+
27+
end program example_cholesky_solve

‎example/linalg/example_solve_chol.f90‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
program example_solve_chol
2+
use stdlib_linalg_constants, only: dp
3+
use stdlib_linalg, only: cholesky, solve_chol, linalg_state_type
4+
implicit none
5+
6+
real(dp) :: A(3,3), L(3,3), b(3), x(3)
7+
type(linalg_state_type) :: state
8+
9+
! Symmetric positive definite matrix
10+
A(1,:) = [4.0_dp, 2.0_dp, 2.0_dp]
11+
A(2,:) = [2.0_dp, 5.0_dp, 1.0_dp]
12+
A(3,:) = [2.0_dp, 1.0_dp, 6.0_dp]
13+
14+
! Right-hand side
15+
b = [1.0_dp, 2.0_dp, 3.0_dp]
16+
17+
! Compute Cholesky factorization: A = L * L^T
18+
call cholesky(A, L, lower=.true., err=state)
19+
if (state%error()) error stop state%print()
20+
21+
! Solve using pre-computed Cholesky factors
22+
call solve_chol(L, b, x, lower=.true., err=state)
23+
if (state%error()) error stop state%print()
24+
25+
print '("Solution: ",*(f8.4,1x))', x
26+
27+
end program example_solve_chol

‎src/lapack/stdlib_linalg_lapack_aux.fypp‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,8 @@ module stdlib_linalg_lapack_aux
4141
public :: stdlib_selctg_${ri}$
4242
#:endfor
4343
public :: handle_potrf_info
44+
public :: handle_potrs_info
45+
public :: handle_posv_info
4446
public :: handle_getri_info
4547
public :: handle_gesdd_info
4648
public :: handle_gesv_info
@@ -1323,6 +1325,65 @@ module stdlib_linalg_lapack_aux
13231325

13241326
end subroutine handle_potrf_info
13251327

1328+
! Cholesky solve (triangular solve with pre-computed factors)
1329+
elemental subroutine handle_potrs_info(this,info,triangle,n,nrhs,lda,ldb,err)
1330+
character(len=*), intent(in) :: this
1331+
character, intent(in) :: triangle
1332+
integer(ilp), intent(in) :: info,n,nrhs,lda,ldb
1333+
type(linalg_state_type), intent(out) :: err
1334+
1335+
! Process output
1336+
select case (info)
1337+
case (0)
1338+
! Success
1339+
case (-1)
1340+
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', &
1341+
triangle,'. should be U/L')
1342+
case (-2)
1343+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
1344+
case (-3)
1345+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size nrhs=',nrhs)
1346+
case (-5)
1347+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid lda=',lda,': should be >=',n)
1348+
case (-7)
1349+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid ldb=',ldb,': should be >=',n)
1350+
case default
1351+
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
1352+
end select
1353+
1354+
end subroutine handle_potrs_info
1355+
1356+
! Cholesky factorization and solve (combined)
1357+
elemental subroutine handle_posv_info(this,info,triangle,n,nrhs,lda,ldb,err)
1358+
character(len=*), intent(in) :: this
1359+
character, intent(in) :: triangle
1360+
integer(ilp), intent(in) :: info,n,nrhs,lda,ldb
1361+
type(linalg_state_type), intent(out) :: err
1362+
1363+
! Process output
1364+
select case (info)
1365+
case (0)
1366+
! Success
1367+
case (-1)
1368+
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', &
1369+
triangle,'. should be U/L')
1370+
case (-2)
1371+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
1372+
case (-3)
1373+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size nrhs=',nrhs)
1374+
case (-5)
1375+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid lda=',lda,': should be >=',n)
1376+
case (-7)
1377+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid ldb=',ldb,': should be >=',n)
1378+
case (1:)
1379+
err = linalg_state_type(this,LINALG_ERROR,'matrix is not positive definite: ', &
1380+
'leading minor of order',info,' is not positive definite')
1381+
case default
1382+
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
1383+
end select
1384+
1385+
end subroutine handle_posv_info
1386+
13261387
elemental subroutine handle_getri_info(this,info,lda,n,err)
13271388
character(len=*), intent(in) :: this
13281389
integer(ilp), intent(in) :: info,lda,n

‎src/linalg/stdlib_linalg.fypp‎

Lines changed: 86 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,9 @@ module stdlib_linalg
4444
public :: mnorm
4545
public :: get_norm
4646
public :: solve
47-
public :: solve_lu
47+
public :: solve_lu
48+
public :: solve_chol
49+
public :: cholesky_solve
4850
public :: solve_lstsq
4951
public :: solve_constrained_lstsq
5052
public :: trace
@@ -457,9 +459,91 @@ module stdlib_linalg
457459
!> [optional] state return flag. On error if not requested, the code will stop
458460
type(linalg_state_type), optional, intent(out) :: err
459461
end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$
462+
#:endfor
463+
#:endfor
464+
end interface solve_lu
465+
466+
! Solve linear system Ax = b using pre-computed Cholesky decomposition (subroutine interface)
467+
interface solve_chol
468+
!! version: experimental
469+
!!
470+
!! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a
471+
!! symmetric positive definite matrix \( A \) that has been pre-factorized using Cholesky.
472+
!! ([Specification](../page/specs/stdlib_linalg.html#solve_chol-solve-a-linear-system-using-cholesky-factors))
473+
!!
474+
!!### Summary
475+
!! Subroutine interface for solving a linear system using pre-computed Cholesky factors.
476+
!!
477+
!!### Description
478+
!!
479+
!! This interface provides methods for computing the solution of a linear matrix system using
480+
!! Cholesky factors. Supported data types include `real` and `complex`. Preallocated space
481+
!! for the solution vector `x` is user-provided. The `lower` argument is REQUIRED and must
482+
!! match the `lower` used during the Cholesky factorization.
483+
!! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`)
484+
!! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
485+
!!
486+
!!@note The solution is based on LAPACK's `*POTRS` routines.
487+
!!
488+
#:for nd,ndsuf,nde in ALL_RHS
489+
#:for rk,rt,ri in RC_KINDS_TYPES
490+
pure module subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$(a,b,x,lower,err)
491+
!> Input matrix a[n,n] containing Cholesky factors from cholesky
492+
${rt}$, intent(in) :: a(:,:)
493+
!> Right hand side vector or array, b[n] or b[n,nrhs]
494+
${rt}$, intent(in) :: b${nd}$
495+
!> Result array/matrix x[n] or x[n,nrhs]
496+
${rt}$, intent(inout), contiguous, target :: x${nd}$
497+
!> Is the lower triangular factor stored? (REQUIRED)
498+
logical(lk), intent(in) :: lower
499+
!> [optional] state return flag. On error if not requested, the code will stop
500+
type(linalg_state_type), optional, intent(out) :: err
501+
end subroutine stdlib_linalg_${ri}$_solve_chol_${ndsuf}$
460502
#:endfor
461503
#:endfor
462-
end interface solve_lu
504+
end interface solve_chol
505+
506+
! One-shot Cholesky factorization and solve (convenience wrapper using POSV)
507+
interface cholesky_solve
508+
!! version: experimental
509+
!!
510+
!! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a
511+
!! symmetric positive definite matrix \( A \). Combines factorization and solve in one call.
512+
!! ([Specification](../page/specs/stdlib_linalg.html#cholesky_solve-one-shot-cholesky-solve))
513+
!!
514+
!!### Summary
515+
!! One-shot factorization and solve for SPD systems (wraps LAPACK POSV).
516+
!!
517+
!!### Description
518+
!!
519+
!! This interface computes both the Cholesky factorization and solves the linear system
520+
!! in a single call. Use this for one-time solves. For repeated solves with the same
521+
!! matrix but different RHS, use `cholesky` + `solve_chol` for better performance.
522+
!! Supported data types include `real` and `complex`.
523+
!! By default, A is not overwritten. Set `overwrite_a=.true.` to allow in-place
524+
!! factorization for better performance.
525+
!!
526+
!!@note The solution is based on LAPACK's `*POSV` routines.
527+
!!
528+
#:for nd,ndsuf,nde in ALL_RHS
529+
#:for rk,rt,ri in RC_KINDS_TYPES
530+
pure module subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$(a,b,x,lower,overwrite_a,err)
531+
!> Input SPD matrix a[n,n]
532+
${rt}$, intent(inout), target :: a(:,:)
533+
!> Right hand side vector or array, b[n] or b[n,nrhs]
534+
${rt}$, intent(in) :: b${nd}$
535+
!> Result array/matrix x[n] or x[n,nrhs]
536+
${rt}$, intent(inout), contiguous, target :: x${nd}$
537+
!> [optional] Use lower triangular factorization? Default = .true.
538+
logical(lk), optional, intent(in) :: lower
539+
!> [optional] Can A data be overwritten and destroyed? Default = .false.
540+
logical(lk), optional, intent(in) :: overwrite_a
541+
!> [optional] state return flag. On error if not requested, the code will stop
542+
type(linalg_state_type), optional, intent(out) :: err
543+
end subroutine stdlib_linalg_${ri}$_cholesky_solve_${ndsuf}$
544+
#:endfor
545+
#:endfor
546+
end interface cholesky_solve
463547

464548
! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized.
465549
interface lstsq

0 commit comments

Comments
 (0)