# chex35.r                         J F Monahan, August 2007
 #                                  revised September 2009
 # Example 3.5 Condition of Hilbert matrix
 #
 # set n = size of matrix
  for (n in(3:8)) {               # not too big, but 2 is a problem 
  print("size of matrix")
  print(n)
 # Hilbert matrix is easy, but not exactly represented
  A <- matrix(0,n,n)              # initialize to zero
  "Hilbert matrix"
  A <- 1/(row(A)+col(A)-1)
 # inverse of matrix has integer values -- exact representation
 # save some constants to make calculation easier
  Is <- 1:n                       # needs n > 2 to work
  bi <- choose(Is+n-1,Is-1)       # save these (use with i)
  bj <- choose(n,Is)              # and these (use with j)
  B <- matrix(0,n,n)              # initialize to zero
    sign <- -1
    for (i in 1:n) {
      for (j in 1:n) {
        sign <- -sign              # (-1)**(i+j)
        B[i,j] <- sign*j*choose(i+j-2,i-1)*bi[i]*choose(j+n-1,n-i)*bj[j]
                      } }   
  print("inverse of Hilbert matrix")
  print(B)
 # is it really the inverse?
  "A %*% B"
  isitI <- A %*% B
  print(c(max(abs(diag(isitI)-1)),max(abs(isitI[diag(n)==0]))))
 # compute p = inf norms (row sum norm)
  normA <- max( apply(abs(A),1,sum) )
  normB <- max( apply(abs(B),1,sum) )
  print("infinity norm of Hilbert matrix and inverse")
  print(c(normA,normB))
  print("true kappa, reciprocal, and estimate")
  tkappa <- normA*normB
  print(c(tkappa,1/tkappa,kappa(B,norm="I",exact=FALSE)))
                                     } # end of loop on n
 # 
  rm(list=ls())