ldl Subroutine

    subroutine ldl(n,a,z,sigma,w)

    implicit none

    integer,intent(in)                  :: n     !! order of the coefficient matrix `a`
    real(wp),intent(in)                 :: sigma !! scalar factor by which the modifying dyade \(z z^T\) is multiplied.
    real(wp),dimension(*),intent(inout) :: a     !! ***In:*** positive definite matrix of dimension `n`;
                                                 !! only the lower triangle is used and is stored column by
                                                 !! column as one dimensional array of dimension `n*(n+1)/2`.
                                                 !!
                                                 !! ***Out:*** updated \(LDL^T\) factors
    real(wp),dimension(*),intent(inout) :: w     !! working array of dimension `n` (used only if \( \sigma \lt 0 \) ).
    real(wp),dimension(*),intent(inout) :: z     !! vector of dimension `n` of updating elements.

    integer :: i , ij , j
    real(wp) :: t , v , u , tp , beta , alpha , delta , gamma

    if ( abs(sigma)>zero ) then
        ij = 1
        t = one/sigma
        if ( sigma<=zero ) then
            ! prepare negative update
            do i = 1 , n
                w(i) = z(i)
            end do
            do i = 1 , n
                v = w(i)
                t = t + v*v/a(ij)
                do j = i + 1 , n
                    ij = ij + 1
                    w(j) = w(j) - v*a(ij)
                end do
                ij = ij + 1
            end do
            if ( t>=zero ) t = epmach/sigma
            do i = 1 , n
                j = n + 1 - i
                ij = ij - i
                u = w(j)
                w(j) = t
                t = t - u*u/a(ij)
            end do
        end if
        ! here updating begins
        do i = 1 , n
            v = z(i)
            delta = v/a(ij)
            if ( sigma<zero ) tp = w(i)
            if ( sigma>zero ) tp = t + delta*v
            alpha = tp/t
            a(ij) = alpha*a(ij)
            if ( i==n ) return
            beta = delta/tp
            if ( alpha>four ) then
                gamma = t/tp
                do j = i + 1 , n
                    ij = ij + 1
                    u = a(ij)
                    a(ij) = gamma*u + beta*z(j)
                    z(j) = z(j) - v*u
                end do
            else
                do j = i + 1 , n
                    ij = ij + 1
                    z(j) = z(j) - v*a(ij)
                    a(ij) = a(ij) + beta*z(j)
                end do
            end if
            ij = ij + 1
            t = tp
        end do
    end if

    end subroutine ldl