function [r,x,it,ritv,itv]=invpow(A, sh=0,verb=0,tol=1e-12, x0=rand(length(A(1,:)),1),maxit=200)
#implementation of inverse power method i.e power method applied to (A-shI)^{-1}
#Input 
#A matrix 
#sh = shift; default =0
#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);   B=A-sh*eye(size(A));
  [Q,R]=qr(B);
  itv(:,1)=xp=x0/norm(x0,2);
  x=R\(Q'*xp);
  it=1;
  n=norm(x,2);
  r=x'*xp; 
  ritv(1)=r=1/r+sh;#r approx of eigenvalue of A
  rp=r+1; 
  xpp=zeros(size(x0));
  er=tol+1;
  while((er>tol)&&(it<=maxit))
     xpp=xp;
     itv(:,it+1)=xp=x/n;
     x=R\(Q'*xp);
     it++;
     n=norm(x,2);
     rp=r;
     r=x'*xp;
     #converting to approx of the eigenvalues of A
     ritv(it)=r=1/r+sh;
     er=abs(r-rp);

     if(verb)
       printf("[%d] r=%17.14f\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

