# orthit1.s                         J F Monahan, June 2001
 # orthogonal iteration for a few eigenvalues and vectors
# first problem
 n <- 7   # dimension of matrix
 p <- 3   # of cols or eigenvals to get
 A <- 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)
 A <- matrix(A,n,n)/10
 "matrix"
 A
 I <- diag(1,n)[,1:p]      # just p columns
 Z <- I
 for( i in 1:50) {
 V <- A %*% Z
 Z <- qr.qy(qr(V),I)
 R <- t(Z) %*% V 
 if( (i %% 5) == 0 ) {     # every fifth iteration
                      print("iteration"); print(i)
                      print("Z matrix"); print(Z)
                      print("R matrix"); print(R)
                     }
                }
 rm(list=ls())