s_utl_cholesky_solve Function

public pure function s_utl_cholesky_solve(a, b, n) result(x)

Solve a * x = b for x using Cholesky factor returned by stdlib's chol(). a : (n,n) symmetric positive-definite b : (n)

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: a(n,n)

input square matrix, symmetric positive-definite
real(kind=wp), intent(in) :: b(n)

input right-hand side vector (n)
integer(kind=i4), intent(in) :: n

length of vectors

Return Value real(kind=wp), (n)

output solution vector (n)


Calls

proc~~s_utl_cholesky_solve~~CallsGraph proc~s_utl_cholesky_solve s_utl_cholesky_solve chol chol proc~s_utl_cholesky_solve->chol

Called by

proc~~s_utl_cholesky_solve~~CalledByGraph proc~s_utl_cholesky_solve s_utl_cholesky_solve proc~f_lin_mahalanobis_core f_lin_mahalanobis_core proc~f_lin_mahalanobis_core->proc~s_utl_cholesky_solve proc~f_lin_mahalanobis f_lin_mahalanobis proc~f_lin_mahalanobis->proc~f_lin_mahalanobis_core proc~s_nlp_hclust_core s_nlp_hclust_core proc~s_nlp_hclust_core->proc~f_lin_mahalanobis_core proc~s_nlp_kmeans_core s_nlp_kmeans_core proc~s_nlp_kmeans_core->proc~f_lin_mahalanobis_core interface~fsml_mahalanobis fsml_mahalanobis interface~fsml_mahalanobis->proc~f_lin_mahalanobis proc~s_nlp_hclust s_nlp_hclust proc~s_nlp_hclust->proc~s_nlp_hclust_core proc~s_nlp_hkmeans_core s_nlp_hkmeans_core proc~s_nlp_hkmeans_core->proc~s_nlp_hclust_core proc~s_nlp_hkmeans_core->proc~s_nlp_kmeans_core proc~s_nlp_kmeans s_nlp_kmeans proc~s_nlp_kmeans->proc~s_nlp_kmeans_core interface~fsml_hclust fsml_hclust interface~fsml_hclust->proc~s_nlp_hclust interface~fsml_kmeans fsml_kmeans interface~fsml_kmeans->proc~s_nlp_kmeans proc~s_nlp_hkmeans s_nlp_hkmeans proc~s_nlp_hkmeans->proc~s_nlp_hkmeans_core interface~fsml_hkmeans fsml_hkmeans interface~fsml_hkmeans->proc~s_nlp_hkmeans

Source Code

pure function s_utl_cholesky_solve(a, b, n) result(x)

! ==== Description
!! Solve a * x = b for x using Cholesky factor returned by stdlib's chol().
!! a : (n,n) symmetric positive-definite
!! b : (n)

! ==== Declarations
  integer(i4), intent(in)  :: n      !! length of vectors
  real(wp)   , intent(in)  :: a(n,n) !! input square matrix, symmetric positive-definite
  real(wp)   , intent(in)  :: b(n)   !! input right-hand side vector (n)
  real(wp)                 :: x(n)   !! output solution vector (n)
  real(wp)                 :: c(n,n) !! lower-triangular Cholesky factor of a
  real(wp)                 :: s      !! temporary accumulator for partial sums in substitutions
  integer(i4)              :: i, j   !! loop indices for rows and columns

! ==== Instructions

  ! compute lower-triangular Cholesky factor c such that a = c * transpose(c)
  c = chol(a, lower = .true., other_zeroed = .true.)

  ! forward substitution: c * y = b  ->  y
  i = 1
  s = b(i)
  x(i) = s / c(i, i)
  do i = 2, n
     s = b(i)
     do j = 1, i - 1
        s = s - c(i, j) * x(j)
     enddo
  x(i) = s / c(i, i)
  enddo

  ! backward substitution: transpose(c) * x = y  ->  x
  i = n
  s = x(i)
  x(i) = s / c(i, i)
  do i = n-1, 1, -1
     s = x(i)
     do j = i + 1, n
        s = s - c(j, i) * x(j)
     enddo
     x(i) = s / c(i, i)
  enddo

end function s_utl_cholesky_solve