Skip to content

Commit d685294

Browse files
perazzperazz
authored
linalg: Cholesky factorization (#840)
2 parents 91dcc50 + 140fc7f commit d685294

File tree

9 files changed

+496
-1
lines changed

9 files changed

+496
-1
lines changed

‎doc/specs/stdlib_linalg.md

Lines changed: 95 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1167,6 +1167,101 @@ Exceptions trigger an `error stop`, unless argument `err` is present.
11671167
{!example/linalg/example_svdvals.f90!}
11681168
```
11691169

1170+
1171+
## `cholesky` - Compute the Cholesky factorization of a rank-2 square array (matrix)
1172+
1173+
### Status
1174+
1175+
Experimental
1176+
1177+
### Description
1178+
1179+
This subroutine computes the Cholesky factorization of a `real` or `complex` rank-2 square array (matrix),
1180+
\( A = L \cdot L^T \), or \( A = U^T \cdot U \). \( A \) is symmetric or complex Hermitian, and \( L \),
1181+
\( U \) are lower- or upper-triangular, respectively.
1182+
The solver is based on LAPACK's `*POTRF` backends.
1183+
1184+
### Syntax
1185+
1186+
`call ` [[stdlib_linalg(module):cholesky(interface)]] `(a, c, lower [, other_zeroed] [, err])`
1187+
1188+
### Class
1189+
Subroutine
1190+
1191+
### Arguments
1192+
1193+
`a`: Shall be a rank-2 square `real` or `complex` array containing the coefficient matrix of size `[n,n]`. It is an `intent(inout)` argument, but returns unchanged if the argument `c` is present.
1194+
1195+
`c` (optional): Shall be a rank-2 square `real` or `complex` of the same size and kind as `a`. It is an `intent(out)` argument, that returns the triangular Cholesky matrix `L` or `U`.
1196+
1197+
`lower`: Shall be an input `logical` flag. If `.true.`, the lower triangular decomposition \( A = L \cdot L^T \) will be performed. If `.false.`, the upper decomposition \( A = U^T \cdot U \) will be performed.
1198+
1199+
`other_zeroed` (optional): Shall be an input `logical` flag. If `.true.`, the unused part of the output matrix will contain zeroes. Otherwise, it will not be accessed. This saves cpu time. By default, `other_zeroed=.true.`. It is an `intent(in)` argument.
1200+
1201+
`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
1202+
1203+
### Return values
1204+
1205+
The factorized matrix is returned in-place overwriting `a` if no other arguments are provided.
1206+
Otherwise, it can be provided as a second argument `c`. In this case, `a` is not overwritten.
1207+
The `logical` flag `lower` determines whether the lower- or the upper-triangular factorization should be performed.
1208+
1209+
Results are returned on the applicable triangular region of the output matrix, while the unused triangular region
1210+
is filled by zeroes by default. Optional argument `other_zeroed`, if `.false.` allows the expert user to avoid zeroing the unused part;
1211+
however, in this case, the unused region of the matrix is not accessed and will usually contain invalid values.
1212+
1213+
Raises `LINALG_ERROR` if the underlying process did not converge.
1214+
Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes.
1215+
Exceptions trigger an `error stop`, unless argument `err` is present.
1216+
1217+
### Example
1218+
1219+
```fortran
1220+
{!example/linalg/example_cholesky.f90!}
1221+
```
1222+
1223+
## `chol` - Compute the Cholesky factorization of a rank-2 square array (matrix)
1224+
1225+
### Status
1226+
1227+
Experimental
1228+
1229+
### Description
1230+
1231+
This is a `pure` functional interface for the Cholesky factorization of a `real` or
1232+
`complex` rank-2 square array (matrix) computed as \( A = L \cdot L^T \), or \( A = U^T \cdot U \).
1233+
\( A \) is symmetric or complex Hermitian, and \( L \), \( U \) are lower- or upper-triangular, respectively.
1234+
The solver is based on LAPACK's `*POTRF` backends.
1235+
1236+
Result matrix `c` has the same size and kind as `a`, and returns the lower or upper triangular factor.
1237+
1238+
### Syntax
1239+
1240+
`c = ` [[stdlib_linalg(module):chol(interface)]] `(a, lower [, other_zeroed])`
1241+
1242+
### Arguments
1243+
1244+
`a`: Shall be a rank-2 square `real` or `complex` array containing the coefficient matrix of size `[n,n]`. It is an `intent(inout)` argument, but returns unchanged if argument `c` is present.
1245+
1246+
`lower`: Shall be an input `logical` flag. If `.true.`, the lower triangular decomposition \( A = L \cdot L^T \) will be performed. If `.false.`, the upper decomposition \( A = U^T \cdot U \) will be performed.
1247+
1248+
`other_zeroed` (optional): Shall be an input `logical` flag. If `.true.`, the unused part of the output matrix will contain zeroes. Otherwise, it will not be accessed. This saves cpu time. By default, `other_zeroed=.true.`. It is an `intent(in)` argument.
1249+
1250+
### Return values
1251+
1252+
Returns a rank-2 array `c` of the same size and kind as `a`, that contains the triangular Cholesky matrix `L` or `U`.
1253+
1254+
Raises `LINALG_ERROR` if the underlying process did not converge.
1255+
Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes.
1256+
Exceptions trigger an `error stop`, unless argument `err` is present.
1257+
1258+
### Example
1259+
1260+
```fortran
1261+
{!example/linalg/example_chol.f90!}
1262+
```
1263+
1264+
11701265
## `.inv.` - Inverse operator of a square matrix
11711266

11721267
### Status
@@ -1196,7 +1291,6 @@ If an exception occurred on input errors, or singular matrix, `NaN`s will be ret
11961291
For fine-grained error control in case of singular matrices prefer the `subroutine` and the `function`
11971292
interfaces.
11981293

1199-
12001294
### Example
12011295

12021296
```fortran

‎example/linalg/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,3 +35,5 @@ ADD_EXAMPLE(svd)
3535
ADD_EXAMPLE(svdvals)
3636
ADD_EXAMPLE(determinant)
3737
ADD_EXAMPLE(determinant2)
38+
ADD_EXAMPLE(cholesky)
39+
ADD_EXAMPLE(chol)

‎example/linalg/example_chol.f90

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
! Cholesky factorization: function interface
2+
program example_chol
3+
use stdlib_linalg, only: chol
4+
implicit none
5+
6+
real, allocatable, dimension(:,:) :: A,L,U
7+
8+
! Set real matrix
9+
A = reshape( [ [6, 15, 55], &
10+
[15, 55, 225], &
11+
[55, 225, 979] ], [3,3] )
12+
13+
! Decompose (lower)
14+
L = chol(A, lower=.true.)
15+
16+
! Compare decomposition
17+
print *, maxval(abs(A-matmul(L,transpose(L))))
18+
19+
! Decompose (upper)
20+
U = chol(A, lower=.false.)
21+
22+
! Compare decomposition
23+
print *, maxval(abs(A-matmul(transpose(U),U)))
24+
25+
end program example_chol

‎example/linalg/example_cholesky.f90

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
! Cholesky factorization: subroutine interface
2+
program example_cholesky
3+
use stdlib_linalg, only: cholesky
4+
implicit none
5+
6+
real, dimension(3,3) :: A,L,U
7+
8+
! Set real matrix
9+
A = reshape( [ [6, 15, 55], &
10+
[15, 55, 225], &
11+
[55, 225, 979] ], [3,3] )
12+
13+
! Decompose (lower)
14+
call cholesky(A, L, lower=.true.)
15+
16+
! Compare decomposition
17+
print *, maxval(abs(A-matmul(L,transpose(L))))
18+
19+
! Decompose (upper)
20+
call cholesky(A, U, lower=.false.)
21+
22+
! Compare decomposition
23+
print *, maxval(abs(A-matmul(transpose(U),U)))
24+
25+
end program example_cholesky

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ set(fppFiles
3333
stdlib_linalg_inverse.fypp
3434
stdlib_linalg_state.fypp
3535
stdlib_linalg_svd.fypp
36+
stdlib_linalg_cholesky.fypp
3637
stdlib_optval.fypp
3738
stdlib_selection.fypp
3839
stdlib_sorting.fypp

‎src/stdlib_linalg.fypp

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,8 @@ module stdlib_linalg
1616
implicit none
1717
private
1818

19+
public :: chol
20+
public :: cholesky
1921
public :: det
2022
public :: operator(.det.)
2123
public :: diag
@@ -49,6 +51,85 @@ module stdlib_linalg
4951
! Export linalg error handling
5052
public :: linalg_state_type, linalg_error_handling
5153

54+
interface chol
55+
!! version: experimental
56+
!!
57+
!! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \).
58+
!! ([Specification](../page/specs/stdlib_linalg.html#chol-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
59+
!!
60+
!!### Summary
61+
!! Pure function interface for computing the Cholesky triangular factors.
62+
!!
63+
!!### Description
64+
!!
65+
!! This interface provides methods for computing the lower- or upper- triangular matrix from the
66+
!! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
67+
!! Supported data types include `real` and `complex`.
68+
!!
69+
!!@note The solution is based on LAPACK's `*POTRF` methods.
70+
!!
71+
#:for rk,rt,ri in RC_KINDS_TYPES
72+
pure module function stdlib_linalg_${ri}$_cholesky_fun(a,lower,other_zeroed) result(c)
73+
!> Input matrix a[m,n]
74+
${rt}$, intent(in) :: a(:,:)
75+
!> [optional] is the lower or upper triangular factor required? Default = lower
76+
logical(lk), optional, intent(in) :: lower
77+
!> [optional] should the unused half of the return matrix be zeroed out? Default: yes
78+
logical(lk), optional, intent(in) :: other_zeroed
79+
!> Output matrix with Cholesky factors c[n,n]
80+
${rt}$ :: c(size(a,1),size(a,2))
81+
end function stdlib_linalg_${ri}$_cholesky_fun
82+
#:endfor
83+
end interface chol
84+
85+
interface cholesky
86+
!! version: experimental
87+
!!
88+
!! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \).
89+
!! ([Specification](../page/specs/stdlib_linalg.html#cholesky-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
90+
!!
91+
!!### Summary
92+
!! Pure subroutine interface for computing the Cholesky triangular factors.
93+
!!
94+
!!### Description
95+
!!
96+
!! This interface provides methods for computing the lower- or upper- triangular matrix from the
97+
!! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
98+
!! Supported data types include `real` and `complex`.
99+
!! The factorization is computed in-place if only one matrix argument is present; or returned into
100+
!! a second matrix argument, if present. The `lower` `logical` flag allows to select between upper or
101+
!! lower factorization; the `other_zeroed` optional `logical` flag allows to choose whether the unused
102+
!! part of the triangular matrix should be filled with zeroes.
103+
!!
104+
!!@note The solution is based on LAPACK's `*POTRF` methods.
105+
!!
106+
#:for rk,rt,ri in RC_KINDS_TYPES
107+
pure module subroutine stdlib_linalg_${ri}$_cholesky_inplace(a,lower,other_zeroed,err)
108+
!> Input matrix a[m,n]
109+
${rt}$, intent(inout), target :: a(:,:)
110+
!> [optional] is the lower or upper triangular factor required? Default = lower
111+
logical(lk), optional, intent(in) :: lower
112+
!> [optional] should the unused half of the return matrix be zeroed out? Default: yes
113+
logical(lk), optional, intent(in) :: other_zeroed
114+
!> [optional] state return flag. On error if not requested, the code will stop
115+
type(linalg_state_type), optional, intent(out) :: err
116+
end subroutine stdlib_linalg_${ri}$_cholesky_inplace
117+
118+
pure module subroutine stdlib_linalg_${ri}$_cholesky(a,c,lower,other_zeroed,err)
119+
!> Input matrix a[n,n]
120+
${rt}$, intent(in) :: a(:,:)
121+
!> Output matrix with Cholesky factors c[n,n]
122+
${rt}$, intent(out) :: c(:,:)
123+
!> [optional] is the lower or upper triangular factor required? Default = lower
124+
logical(lk), optional, intent(in) :: lower
125+
!> [optional] should the unused half of the return matrix be zeroed out? Default: yes
126+
logical(lk), optional, intent(in) :: other_zeroed
127+
!> [optional] state return flag. On error if not requested, the code will stop
128+
type(linalg_state_type), optional, intent(out) :: err
129+
end subroutine stdlib_linalg_${ri}$_cholesky
130+
#:endfor
131+
end interface cholesky
132+
52133
interface diag
53134
!! version: experimental
54135
!!

0 commit comments

Comments
 (0)