# chlsky.s                         J F Monahan, June 2001
 # Cholesky factorization of a positive definite matrix
A<- c(4, 2, 2, 4, 2, 5, 7, 0, 2, 7, 19, 11, 4, 0, 11, 25 )
 A <- matrix(A,4,4)
 "matrix"
 A
 Lt <- chol(A)
 "Cholesky factor -- transposed to be upper triangular"
 Lt
 det <-prod(diag(Lt))
 "determinant"
 det*det
  adjust <- function(x) {
     if( x <= 0 ) c( 0, -2**31 ) else {
          i <- 0
          while (x > 16 ) { 
             x <- x/16
             i <- i + 4
                          }
          while (x < 1 ) { 
             x <- x*16
             i <- i - 4
                          }
          c(x,i)
                                      }
                                      }
 adjust(det*det)
                 
 b <-c( -1, 1, 2.5, .25)
 "Right hand side"
 b
 y <- solve(t(Lt),b)
 "intermediate"
 y
 x <- backsolve(Lt,y)
 "Solution"
 x
 rm(list=ls())