# c3nwdemo.r         New demo for Chapter 3
 #                           September 2007, revised 2009
 #
 # simpler example for condition number for solving systems
 #  of linear equations
 #
 # here's a matrix
  A <- matrix(c(4,-2,2,4,-2,2,-2,-2,2,-2,11,5,4,-2,5,6),4,4)
  A                               # write it out
 #
 # Cholesky factorization
 #
  Lt <- chol(A)
  Lt                              # write it out
 #
 # rhs
 #
  b <- c(2,0,0,2)
 #
 # solve these equations
 #
  y <- forwardsolve( t(Lt), b)
  x <- backsolve( Lt, y )
  x                               # write it out
 #
 # check the solution
 #
  A %*% x
 #
 # get infinity norm (from homework exercise)
 #
  normA <- max( apply( abs(A), 1, sum) )
 #
 # get infinity norm of inverse for condition number
 #
  kappa <- normA*max( apply( abs(chol2inv(Lt)),1,sum) )
  kappa
 #
 # small change in matrix A
 #
  E <- matrix( 1.e-4, 4, 4)   # recycling
 #
 # Cholesky factor of changed matrix
 #
  Lte <- chol( A+E )
  Lte
 #
 # solve equations
 #
  y <- forwardsolve(t(Lte),b)
  xe <- backsolve(Lte,y)
  xe
 #
 # norm of difference
 # 
  max( abs(x-xe) ) 
 #
 #
 # small change in rhs
 #
  be <- b + 1.e-4             # recycling
 #
 # solve equations
 #
  y <- forwardsolve(t(Lt),be)
  xe <- backsolve(Lt,y)
  xe
 #
 # norm of difference
 # 
  max( abs(x-xe) ) 
 #
 # done
  rm(list=ls())