1==1;
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=%g\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

function [r,x,it,ritv,itv]=invpow(A, sh=0,verb=0,tol=1e-12, x0=rand(length(A(1,:)),1),maxit=200)
#function [r,x,it]=powerr(A, tol, x0,maxit)
#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=%g\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

#zad 1 
A1=[3,1;1,3]
d=eig(A1) #pary wl (4,[1;1])  (2,[1,-1])
rat=max(d)/min(d);rat=rat*rat
[r,x,it1,ritv]=powerr(A1);
r
r-4
err1=abs(ritv-4);
x
norm(x-[x(1);x(1)])#znamy pare wlasna!
#czy wyniki z ponizszego wektora zbizaja sie do rat?
err1(1:it1-1)./err1(2:it1)

A2=[-4,-3;-3,-4]
d=eig(A2) #wektory wlasne [1;1] [1;-1]
rat=max(abs(d))/min(abs(d));rat=rat*rat
[r,x,it2,ritv]=powerr(A2);
r
r+7
err2=abs(ritv+7);
x
norm(x-[x(1);x(1)])#znamy pare wlasna!
err2(1:it2-1)./err2(2:it2)

A3=[100,1;1,100]
d=eig(A3) #pary wl (101,[1;1])  (-99,[1,-1])
rat=max(d)/min(d);rat=rat*rat
[r,x,it3,ritv]=powerr(A3);
r
r-101
err3=abs(ritv-101);
err3(180:it3-1)./err3(181:it3)


semilogy(1:it1,err1,1:it2,err2,1:it3,err3);
#prosze zgadnac ktory wykres odpowiada ktoremu przypadkowi....

#semilogy(1:it1,err1,";1;",1:it2,err2,";2;",1:it3,err3,";3;")

#pause(1);
#porownaj ilosci iteracji


#zad 2 testy odwrotnej potegowej
[r,x,it1,ritv]=invpow(A3, sh=100.5);
r 
er1=abs(ritv-101);
[r,x,it2,ritv]=invpow(A3, sh=99.3); 
r
er2=abs(ritv-99);
[r,x,it3,ritv]=invpow(A3, sh=99.1);
r
er3=abs(ritv-99);
semilogy(1:it1,er1,";1;",1:it2,er2,";2;",1:it3,er3,";3;")




