#  Mathematics 172
#  Numerical Analysis 2
#  April 20, 1999
#
#  The basic power method for eigenvalues -- returns
#  an approximation to the (real) dominant eigenvalue of a 
#  matrix and a corresponding eigenvector
#
#  usage:  
# 
#  PowerM(A,tol,maxiter);
#
#    where:  A = the matrix (should be a square matrix of course!)
#            tol = stopping tolerance, we continue until successive
#                  iterates produce approximate eigenvectors differing
#                  by less than tol in the infinity-norm.
#            maxiter = maximum number of iterations
#
with(linalg): with(stats):
PowerM := proc(A,tol,maxiter)

   local i, j, oldvec, newvec, newvecn, n, notdone;
   
   n := rowdim(A);
#   
#  the next lines generate a random starting vector for the power method
#  iteration, and normalize it (infinity norm = 1)
#   
   randomize();
   oldvec := convert([stats[random,normald](n)],array);  
   oldvec := scalarmul(oldvec,1/norm(oldvec,infinity));
#
#  power method iteration;
#  at each step, newvecn is a vector of infinity norm = 1
#  in the direction of newvec.  We DO NOT simply divide
#  by the infinity norm, since that produces alternating
#  signs in the case of a negative dominant eigenvalue(!)
#
   notdone := true;
   i := 1;
   while (i <= maxiter) and notdone do
     newvec := multiply(A,oldvec);
     for j to n do
        if abs(newvec[j]) = norm(newvec,infinity) then
           break
        fi
     od;
     newvecn:=scalarmul(newvec,1/newvec[j]);
     if norm(evalm(newvecn-oldvec),infinity) < tol then
        notdone := false
     else
        i := i + 1;
        oldvec:=eval(newvecn);
     fi;
   od;
   if i <= maxiter then
     RETURN(newvec[1]/oldvec[1],eval(newvecn))
   else
     ERROR("power method did not converge")
   fi
   end: