#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
! Cholesky factorization of a matrix, based on LAPACK *POTRF functions
submodule (stdlib_linalg) stdlib_linalg_cholesky
use stdlib_linalg_constants
use stdlib_linalg_lapack, only: potrf
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
implicit none(type,external)
character(*), parameter :: this = 'cholesky'
contains
elemental subroutine handle_potrf_info(info,triangle,lda,n,err)
character, intent(in) :: triangle
integer(ilp), intent(in) :: info,lda,n
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 (-4)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid lda=',lda,': is < n = ',n)
case (1:)
err = linalg_state_type(this,LINALG_ERROR,'cannot complete factorization:',info, &
'-th order leading minor is not positive definite')
case default
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
end select
end subroutine handle_potrf_info
#:for rk,rt,ri in RC_KINDS_TYPES
! Compute the Cholesky factorization of a symmetric / Hermitian matrix, A = L*L^T = U^T*U.
! The factorization is returned in-place, overwriting matrix A
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
!> Local variables
type(linalg_state_type) :: err0
integer(ilp) :: lda,n,info,j
logical(lk) :: lower_,other_zeroed_
character :: triangle
${rt}$, parameter :: zero = 0.0_${rk}$
!> Check if the lower or upper factor is required.
!> Default: use lower factor
lower_ = .true.
if (present(lower)) lower_ = lower
triangle = merge('L','U',lower_)
!> Check if the unused half of the return matrix should be zeroed out (default).
!> Otherwise it is unused and will contain garbage.
other_zeroed_ = .true.
if (present(other_zeroed)) other_zeroed_ = other_zeroed
!> Problem size
lda = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
! Check sizes
if (n<1 .or. lda<1 .or. lda<n) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR, &
'invalid matrix size: a(m,n)=',[lda,n])
else
! Compute factorization
call potrf(triangle,n,a,lda,info)
call handle_potrf_info(info,triangle,lda,n,err0)
! Zero-out the unused part of matrix A
clean_unused: if (other_zeroed_ .and. err0%ok()) then
if (lower_) then
forall (j=2:n) a(:j-1,j) = zero
else
forall (j=1:n-1) a(j+1:,j) = zero
endif
endif clean_unused
endif
! Process output and return
call linalg_error_handling(err0,err)
end subroutine stdlib_linalg_${ri}$_cholesky_inplace
! Compute the Cholesky factorization of a symmetric / Hermitian matrix, A = L*L^T = U^T*U.
! The factorization is returned as a separate matrix
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
type(linalg_state_type) :: err0
integer(ilp) :: lda,n,ldc,nc
! Check C sizes
lda = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
ldc = size(c,1,kind=ilp)
nc = size(c,2,kind=ilp)
if (lda<1 .or. n<1 .or. lda<n .or. ldc<n .or. nc<n) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR, &
'invalid matrix sizes: a=',[lda,n], &
' c=',[ldc,nc])
else
! Copy data in
c(:n,:n) = a(:n,:n)
! Get cholesky factors
call stdlib_linalg_${ri}$_cholesky_inplace(c,lower,other_zeroed,err0)
end if
! Process output and return
call linalg_error_handling(err0,err)
end subroutine stdlib_linalg_${ri}$_cholesky
! Compute the Cholesky factorization of a symmetric / Hermitian matrix, A = L*L^T = U^T*U.
! Function interface
pure module function stdlib_linalg_${ri}$_cholesky_fun(a,lower,other_zeroed) result(c)
!> Input matrix a[n,n]
${rt}$, intent(in) :: a(:,:)
!> Output matrix with Cholesky factors c[n,n]
${rt}$ :: c(size(a, 1),size(a, 2))
!> [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
c = a
call stdlib_linalg_${ri}$_cholesky_inplace(c,lower,other_zeroed)
end function stdlib_linalg_${ri}$_cholesky_fun
#:endfor
end submodule stdlib_linalg_cholesky
