# chlsky.r                         J F Monahan, June 2001
 #                                  revised August 2007
 # Cholesky factorization of a positive definite matrix
 A <- matrix(c(4,2,2,4, 2,5,7,0, 2,7,19,11, 4,0,11,25), 4, 4)
 "matrix"
 A
 Lt <- chol(A)
 "Cholesky factor -- transposed to be upper triangular"
 Lt
 det <-prod(diag(Lt))
 "determinant"
 det*det
 # use function 'adjust' previously constructed
  adjust <- function(x) {
     if( x <= 0 ) c( 0, -2**31 ) else {  # flag if not positive
          i <- 0
          while (x > 16 ) { 
             x <- x/16
             i <- i + 4
                          }
          while (x < 1 ) { 
             x <- x*16
             i <- i - 4
                          }
          c(x,i)
                                      }
                                      }
 # now compute -- it may have already overflowed
 adjust(det*det)
                 
 b <-c( -1, 1, 2.5, .25)
 "Right hand side"
 b
 y <- forwardsolve(t(Lt),b)            # lower triangular system
 "intermediate"
 y
 x <- backsolve(Lt,y)                  # upper triangular system
 "Solution"
 x
 rm(list=ls())