# orthit4.s                         J F Monahan, June 2001
 # orthogonal iteration for a few eigenvalues and vectors
# fourth problem
 n <- 5   # dimension of matrix
 p <- 3   # of cols or eigenvals to get
 A <- c( .6561, .2916, .0486, .0036, .0001,
         .2401, .4116, .2646, .0756, .0081, 
         .0625, .2500, .3750, .2500, .0625,
         .0081, .0756, .2646, .4116, .2401, 
         .0001, .0036, .0486, .2916, .6561)
 A <- matrix(A,n,n)
 "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())