linalg: Cholesky factorization

doc/specs/stdlib_linalg.md

@@ -1167,6 +1167,101 @@ Exceptions trigger an `error stop`, unless argument `err` is present.
 {!example/linalg/example_svdvals.f90!}
 ```
 
+
+## `cholesky` - Compute the Cholesky factorization of a rank-2 square array (matrix)
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the Cholesky factorization of a `real` or `complex` rank-2 square array (matrix), 
+\( A = L \cdot L^T \), or \( A = U^T \cdot U \). \( A \) is symmetric or complex Hermitian, and \( L \), 
+\( U \) are lower- or upper-triangular, respectively. 
+The solver is based on LAPACK's `*POTRF` backends.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):cholesky(interface)]] `(a, c, lower [, other_zeroed] [, err])`
+
+### Class
+Subroutine
+
+### Arguments
+
+`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.
+
+`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`.
+
+`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.
+
+`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.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
+
+### Return values
+
+The factorized matrix is returned in-place overwriting `a` if no other arguments are provided.
+Otherwise, it can be provided as a second argument `c`. In this case, `a` is not overwritten.
+The `logical` flag `lower` determines whether the lower- or the upper-triangular factorization should be performed. 
+
+Results are returned on the applicable triangular region of the output matrix, while the unused triangular region 
+is filled by zeroes by default. Optional argument `other_zeroed`, if `.false.` allows the expert user to avoid zeroing the unused part; 
+however, in this case, the unused region of the matrix is not accessed and will usually contain invalid values. 
+
+Raises `LINALG_ERROR` if the underlying process did not converge.
+Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes.
+Exceptions trigger an `error stop`, unless argument `err` is present.
+
+### Example
+
+```fortran
+{!example/linalg/example_cholesky.f90!}
+```
+
+## `chol` - Compute the Cholesky factorization of a rank-2 square array (matrix)
+
+### Status
+
+Experimental
+
+### Description
+
+This is a `pure` functional interface for the Cholesky factorization of a `real` or
+`complex` rank-2 square array (matrix) computed as \( A = L \cdot L^T \), or \( A = U^T \cdot U \). 
+\( A \) is symmetric or complex Hermitian, and \( L \), \( U \) are lower- or upper-triangular, respectively. 
+The solver is based on LAPACK's `*POTRF` backends.
+
+Result matrix `c` has the same size and kind as `a`, and returns the lower or upper triangular factor. 
+
+### Syntax
+
+`c = ` [[stdlib_linalg(module):chol(interface)]] `(a, lower [, other_zeroed])`
+
+### Arguments
+
+`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.
+
+`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.
+
+`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.
+
+### Return values
+
+Returns a rank-2 array `c` of the same size and kind as `a`, that contains the triangular Cholesky matrix `L` or `U`.
+
+Raises `LINALG_ERROR` if the underlying process did not converge.
+Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes.
+Exceptions trigger an `error stop`, unless argument `err` is present.
+
+### Example
+
+```fortran
+{!example/linalg/example_chol.f90!}
+```
+
+
 ## `.inv.` - Inverse operator of a square matrix
 
 ### Status
@@ -1196,7 +1291,6 @@ If an exception occurred on input errors, or singular matrix, `NaN`s will be ret
 For fine-grained error control in case of singular matrices prefer the `subroutine` and the `function`
 interfaces.
 
-
 ### Example
 
 ```fortran

example/linalg/CMakeLists.txt

@@ -35,3 +35,5 @@ ADD_EXAMPLE(svd)
 ADD_EXAMPLE(svdvals)
 ADD_EXAMPLE(determinant)
 ADD_EXAMPLE(determinant2)
+ADD_EXAMPLE(cholesky)
+ADD_EXAMPLE(chol)

example/linalg/example_chol.f90

@@ -0,0 +1,25 @@
+! Cholesky factorization: function interface 
+program example_chol
+  use stdlib_linalg, only: chol
+  implicit none
+
+  real, allocatable, dimension(:,:) :: A,L,U 
+
+  ! Set real matrix
+  A = reshape( [ [6, 15, 55], &
+                 [15, 55, 225], &
+                 [55, 225, 979] ], [3,3] )
+
+  ! Decompose (lower)
+  L = chol(A, lower=.true.)
+
+  ! Compare decomposition 
+  print *, maxval(abs(A-matmul(L,transpose(L))))
+
+  ! Decompose (upper)
+  U = chol(A, lower=.false.)
+
+  ! Compare decomposition 
+  print *, maxval(abs(A-matmul(transpose(U),U)))
+
+end program example_chol

example/linalg/example_cholesky.f90

@@ -0,0 +1,25 @@
+! Cholesky factorization: subroutine interface 
+program example_cholesky
+  use stdlib_linalg, only: cholesky
+  implicit none
+
+  real, dimension(3,3) :: A,L,U 
+
+  ! Set real matrix
+  A = reshape( [ [6, 15, 55], &
+                 [15, 55, 225], &
+                 [55, 225, 979] ], [3,3] )
+
+  ! Decompose (lower)
+  call cholesky(A, L, lower=.true.)
+
+  ! Compare decomposition 
+  print *, maxval(abs(A-matmul(L,transpose(L))))
+
+  ! Decompose (upper)
+  call cholesky(A, U, lower=.false.)
+
+  ! Compare decomposition 
+  print *, maxval(abs(A-matmul(transpose(U),U)))
+
+end program example_cholesky

src/CMakeLists.txt

@@ -33,6 +33,7 @@ set(fppFiles
     stdlib_linalg_inverse.fypp
     stdlib_linalg_state.fypp
     stdlib_linalg_svd.fypp 
+    stdlib_linalg_cholesky.fypp
     stdlib_optval.fypp
     stdlib_selection.fypp
     stdlib_sorting.fypp

src/stdlib_linalg.fypp

@@ -16,6 +16,8 @@ module stdlib_linalg
   implicit none
   private
 
+  public :: chol
+  public :: cholesky
   public :: det
   public :: operator(.det.)
   public :: diag
@@ -49,6 +51,85 @@ module stdlib_linalg
   ! Export linalg error handling
   public :: linalg_state_type, linalg_error_handling
 
+  interface chol
+    !! version: experimental 
+    !!
+    !! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \). 
+    !! ([Specification](../page/specs/stdlib_linalg.html#chol-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
+    !! 
+    !!### Summary 
+    !! Pure function interface for computing the Cholesky triangular factors. 
+    !!
+    !!### Description
+    !! 
+    !! This interface provides methods for computing the lower- or upper- triangular matrix from the 
+    !! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
+    !! Supported data types include `real` and `complex`.    
+    !! 
+    !!@note The solution is based on LAPACK's `*POTRF` methods.
+    !!     
+     #:for rk,rt,ri in RC_KINDS_TYPES
+     pure module function stdlib_linalg_${ri}$_cholesky_fun(a,lower,other_zeroed) result(c)
+         !> Input matrix a[m,n]
+         ${rt}$, intent(in) :: a(:,:)
+         !> [optional] is the lower or upper triangular factor required? Default = lower
+         logical(lk), optional, intent(in) :: lower
+         !> [optional] should the unused half of the return matrix be zeroed out? Default: yes
+         logical(lk), optional, intent(in) :: other_zeroed
+         !> Output matrix with Cholesky factors c[n,n]
+         ${rt}$ :: c(size(a,1),size(a,2))                  
+     end function stdlib_linalg_${ri}$_cholesky_fun
+     #:endfor    
+  end interface chol  
+
+  interface cholesky
+    !! version: experimental 
+    !!
+    !! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \). 
+    !! ([Specification](../page/specs/stdlib_linalg.html#cholesky-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
+    !! 
+    !!### Summary 
+    !! Pure subroutine interface for computing the Cholesky triangular factors. 
+    !!
+    !!### Description
+    !! 
+    !! This interface provides methods for computing the lower- or upper- triangular matrix from the 
+    !! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
+    !! Supported data types include `real` and `complex`.    
+    !! The factorization is computed in-place if only one matrix argument is present; or returned into 
+    !! a second matrix argument, if present. The `lower` `logical` flag allows to select between upper or 
+    !! lower factorization; the `other_zeroed` optional `logical` flag allows to choose whether the unused
+    !! part of the triangular matrix should be filled with zeroes.
+    !! 
+    !!@note The solution is based on LAPACK's `*POTRF` methods.
+    !!         
+     #:for rk,rt,ri in RC_KINDS_TYPES
+     pure module subroutine stdlib_linalg_${ri}$_cholesky_inplace(a,lower,other_zeroed,err) 
+         !> Input matrix a[m,n]
+         ${rt}$, intent(inout), target :: a(:,:)
+         !> [optional] is the lower or upper triangular factor required? Default = lower
+         logical(lk), optional, intent(in) :: lower
+         !> [optional] should the unused half of the return matrix be zeroed out? Default: yes
+         logical(lk), optional, intent(in) :: other_zeroed
+         !> [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_inplace
+
+     pure module subroutine stdlib_linalg_${ri}$_cholesky(a,c,lower,other_zeroed,err) 
+         !> Input matrix a[n,n]
+         ${rt}$, intent(in) :: a(:,:)
+         !> Output matrix with Cholesky factors c[n,n]
+         ${rt}$, intent(out) :: c(:,:)         
+         !> [optional] is the lower or upper triangular factor required? Default = lower
+         logical(lk), optional, intent(in) :: lower
+         !> [optional] should the unused half of the return matrix be zeroed out? Default: yes
+         logical(lk), optional, intent(in) :: other_zeroed
+         !> [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
+     #:endfor
+  end interface cholesky
+
   interface diag
     !! version: experimental
     !!