function[levl,levc,lerr,nit,pvl,pvc,perr,ok] = power(mat,x0,maxit,tol,levl,lerr)
%%  [levl,levc,lerr,nit,pvl,pvc,perr,ok] = power(mat,x0,maxit,tol,levl,lerr)
%% Perform up to maxit iterations of power method 
%%  Inputs: 
%%   Required: mat - Square matrix. Eigenvalue of maximum modulus is to be
%%		     estimated.
%%   	       x0  - Column vector to use as initial eigenvector guess.
%%   Optional: maxit - maximum number of iterations, maxit=20 if missing
%%		       tol - relative tolerance for error in eigenvalue,
%%		       tol = 1e-8 if missing,
%%		     iterations stop if difference between two successive 
%%			eigenvalue estimates differ by < tol*norm(mat,1)
%%	       levl  - eigenvalue estimate associated with use of x0
%%			as approximate eigenvector. Obtained from
%%		     	previous call to power when x0 = levc.
%%	       lerr  - error estimate associated with levl as eigenvalue
%%		     	estimate.  Obtained from previous call to power.
%%   Ouputs:   levl -  eigenvalue estimate from last iteration performed.
%%			In the current encoding pvc is the normalized 
%%			previous eigenvector estimate and 
%%		   	  levl = mat*pvc(1)/pvc(1)  (ratio of first components)
%%     	       levc -  normalized eigenvector estimate   mat*pvc/norm(mat*pvc,1)
%%     	       lerr -  estimate of error in levl estimated as levl - pvl,   
%%		       	where pvl is eigenvalue estimate of from previous 
%%		     	iteration.  sgn(levl)*norm(mat,1) returned if pvl not
%%		        available. 
%%	     nit  -  Number of iterations performed
%%	     pvl,pvc,per - levl,levc,lerr values from nit-1'st iteration
%		     if available.
%	     ok   -  Logical:  TRUE (1) only if at most maxit iterations 
%		     	performed and abs(lerr)<tol*norm(mat,1).
%
a=mat;
pvc=x0/(sign(x0(1))*norm(x0,2));
if nargin == 2
	maxit=20;
	tol=1e-8;
elseif nargin==3
	tol=1e-8;
end
if nargin<5
pvl=Inf;
perr=Inf;
elseif nargin<6
pvl=levl;
perr=Inf;
elseif nargin<7
pvl=levl;
perr=lerr;
end

mnrm = norm(mat,1);
nit=0;
if mnrm==0
	levl=pvl;
	levc=pvc;
	lerr=0;
	ok = 1;
	return;
else
	ok = 0;
end
while nit<maxit ,
	if nit > 0 ,
		pvc = levc;
		pvl = levl; 
		perr = lerr;
	end
	levc = mat*pvc; 		% Iterate
	levl = eigvlest(levc,pvc);	% Eigenvalue estimate 
	levc = levc./(sign(levc(1))*norm(levc,2));       % Normalize eigenvector estimate
	if nit ~= 0 | nargin>=5 
		lerr = levl-pvl;
	else
		lerr = sign(levl)*mnrm;
	end
	ok = (abs(lerr)< tol*mnrm );
	nit = nit + 1;
	if ok
		 break 
	end
end;	 
