#  orthit.r
 #
 # orthogonal iteration for largest eigenvalues & vectors
 #
 # uses modified Gram-Schmidt
   mgs <- function(X) {
     R <- matrix(0,ncol(X),ncol(X))     # matrices of zeros to start
     Q <- matrix(0,nrow(X),ncol(X))     #
     for (k in 1:ncol(X)) {
       s <- sqrt( sum(X[,k]^2) )        # norm
       R[k,k] <- s
       Q[,k] <- X[,k]/s                  # normalize
       if(k == ncol(X)) next
       l <-((k+1):ncol(X))
       R[k,l] <- t(Q[,k])%*%X[,l]           # get coefficients
       X[,l] <- X[,l]-outer(Q[,k],R[k,l])}    # residuals
     mgs <- list(Q,R)}                     # return a list
  #
     X <- cbind(rep(1,6),(1:6),(1:6)^2)
     X
     mgsX <- mgs(X)
     mgsX[[1]]
     mgsX[[2]]
  #  assign matrix
    print("first matrix")               # sparse transition matrix
    A <- matrix( c(.9,0,.1,0,0,0,0, .3,.6,.0,.1,0,0,0,
      0.,.2,.7,0,.1,0,0, 0,0,.3,.6,0,.1,0, 0,0,0,.2,.7,0,.1,
      0,0,0,0,.3,.7,0, 0,0,0,0,0,.2,.8)  ,7,7)   
    A                                   # print it out
    Z <- (diag(1,7))[,1:3]              # starting orthonormal vectors
    for (i in (1:40)) {
    V <- A %*% Z                        # multiply
    qrnew <- mgs(V)                     # normalize via MGS
    Z <- qrnew[[1]]
    if( (i %% 5) == 0 ) {               # print a few
    print(i)
    print(Z)
    print(qrnew[[2]])} }
    print("second matrix")
    A <- matrix(c(1,-3,-2,1,-3,10,-3,6,-2,-3,3,-2,1,6,-2,1),4,4)
    A                                   # print it out
    Z <- (diag(1,4))[,1:3]              # starting orthonormal vectors
    for (i in (1:40)) {
    V <- A %*% Z                        # multiply
    qrnew <- mgs(V)                     # normalize via MGS
    Z <- qrnew[[1]]
    if( (i %% 5) == 0 ) {               # print a few
    print(i)
    print(Z)
    print(qrnew[[2]])} }
    print("third matrix")
    A <- matrix(c(8,-3,-2,0,6,5, -3,6,1,-6,0,-2, -2,1,5,-5,2,0,
          0,-6,-5,8,-3,-2, 6,0,2,-3,6,1, 5,-2,0,-2,1,5),6,6)
    A                                   # print it out
    Z <- (diag(1,6))[,1:3]              # starting orthonormal vectors
    for (i in (1:40)) {
    V <- A %*% Z                        # multiply
    qrnew <- mgs(V)                     # normalize via MGS
    Z <- qrnew[[1]]
    if( (i %% 5) == 0 ) {               # print a few
    print(i)
    print(Z)
    print(qrnew[[2]])} }
    print("fourth matrix")
    A <- matrix(0,5,5)
    A <- dbinom(row(A)-1,4,prob=(2*col(A)-1)/10)
    A                                   # print it out
    Z <- (diag(1,5))[,1:3]              # starting orthonormal vectors
    for (i in (1:40)) {
    V <- A %*% Z                        # multiply
    qrnew <- mgs(V)                     # normalize via MGS
    Z <- qrnew[[1]]
    if( (i %% 5) == 0 ) {               # print a few
    print(i)
    print(Z)
    print(qrnew[[2]])} }
  rm(list=ls())     # done
  q()