choldc Subroutine

public subroutine choldc(a, n, np, p, info)

Note

Given a positive definite symmetric matrix a(1:n,1:n), with physical dimensions np, this routine constructs its Cholesky decomposition, A=L L^T. On input, only the upper triangle of a need to be given; it is not modified. The Cholesky factor L is returned in the lower triangle of a, except for its diagonal elements which are returned in p(1:n).

(c) Numerical recipes

Arguments

Type IntentOptional Attributes Name
real(kind=R8), intent(inout), dimension(np, np) :: a

system matrix
integer(kind=I4), intent(in) :: n

system size
integer(kind=I4), intent(in) :: np

matrix size
real(kind=R8), intent(out), dimension(np) :: p

diagonal elements
integer(kind=I4), intent(out) :: info

information ouput

Called by

proc~~choldc~~CalledByGraph proc~choldc choldc proc~moindres_carres moindres_carres proc~moindres_carres->proc~choldc proc~moindres_carres_lineaire moindres_carres_lineaire proc~moindres_carres_lineaire->proc~choldc program~test_utils test_utils program~test_utils->proc~choldc proc~least_squares_tcheby least_squares_tcheby proc~least_squares_tcheby->proc~moindres_carres_lineaire program~test_least test_least program~test_least->proc~moindres_carres program~test_least->proc~moindres_carres_lineaire program~test_tchebychev test_tchebychev program~test_tchebychev->proc~least_squares_tcheby

Source Code

   subroutine choldc (a, n, np, p, info)
   !================================================================================================
   !<@note Given a positive definite symmetric matrix a(1:n,1:n), with
   !< physical dimensions np, this routine constructs its Cholesky
   !< decomposition, A=L L^T. On input, only the upper triangle of
   !< a need to be given; it is not modified. The Cholesky factor L
   !< is returned in the lower triangle of a, except for its diagonal
   !< elements which are returned in p(1:n).
   !<
   !< (c) Numerical recipes
   !<
   !<@endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in   )                    :: n     !! *system size*
   integer(kind=I4), intent(in   )                    :: np    !! *matrix size*
   integer(kind=I4), intent(  out)                    :: info  !! *information ouput*
   real   (kind=R8), intent(inout), dimension(np, np) :: a     !! *system matrix*
   real   (kind=R8), intent(  out), dimension(np)     :: p     !! *diagonal elements*

      integer(kind=I4) :: i, j, k
      real   (kind=R8) :: ssum

      info = 0
o:    do i = 1, n
         do j = i, n
            ssum = a(i,j)
            do k = i-1, 1, -1
               ssum = ssum -a(i,k)*a(j,k)
            enddo
            if(i==j) then
               if (ssum<=0.) then
                  info = 1
                  p(np) = 0._R8
                  p(1) = -ssum
                  exit o
               endif
               p(i) = sqrt(ssum)
            else
               a(j, i) = ssum/p(i)
            endif
         enddo
      enddo o
   return
   endsubroutine choldc