function [r,x,it,ritv,itv]=powerr(A, verb=0,tol=1e-12, x0=rand(length(A(1,:)),1),maxit=200)
#function [r,x,it]=powerr(A, tol, x0,maxit)
#implementation of power method
#Input 
#A matrix 
#verb = 0 without displaying temparary output ; 1 displaying  Rayleigh quotients; 2 displying R. quotient and iterations
#x0 -starting vector
#tol - stopping criteria |r_k-r_{k-1}|<tol
#maxit - maximal no of iterations
#output 
#r - Raileigh quotient - approx of eigenvalue
#x- approximation of eigenvector for the largest eigenvalue
#it - no of iterations
#ritv- vector with succesive Reyleigh quotients
#itv - matrix with columns x_k - x_k power iterations 

  ritv=zeros(maxit,1);
  itv=zeros(length(A(:,1)),maxit); 
  itv(:,1)=xp=x0/norm(x0,2);
  x=A*xp;
  it=1;
  n=norm(x,2);
  ritv(1)=r=x'*xp; 
  rp=r+1; 
  xpp=zeros(size(x0));
  er=tol+1;
  while((er>tol)&&(it<=maxit))
     xpp=xp;
     itv(:,it+1)=xp=x/n;
     x=A*xp;
     it++;
     n=norm(x,2);
     rp=r;
     ritv(it)=r=x'*xp;
     er=abs(r-rp);
     if(verb)
       printf("[%d] r=%5.16f\n",it,r);
       if(verb>1)
        disp(x/n);
       endif
     endif 
  endwhile
  if(it>maxit)
    printf("\n\nWarning. Maximal no of iterations (%d)  is exceeded\n\n",maxit);
  endif
  x=x/n;
  ritv=ritv(1:it);
  itv=itv(:,1:it);
endfunction
