
#tests of svd decomposition
disp("\n\n test svd - proste macierze 3x2 2x3 3x3 ");
A1=[1 2 3; 1 -1 0]
A2=[1 1;1 2; -3 4]
[U,S,V]=svd(A1)
norm(U*S*V'-A1,1)


[U,S,V]=svd(A2)
norm(U*S*V'-A2,1)


A3=[1 1 2; 1 0 2; 2 1 4]

[U,S,V]=svd(A3)
norm(U*S*V'-A3,1)

svd(A3)
#ostatnia wartosc szczeg. bliska zero...
rank(A3)

#macierz rzedu 1 A=u*v'
disp("\n\n macierz rzedu jeden - 5x4 ");
m=5;n=4;
u=randn(m,1);
v=randn(n,1);
A=u*v'
svd(A)
rank(A)
NA=null(A) #kolumny - baza ortonormalna jadra A
RA=orth(A) #kolumny baza ortonormalna R(A)
[U,S,V]=svd(A)
#porownujemy
NA-V(:,2:n)
RA-u/norm(u) #lub RA+u/norm(u)

#rozwiazujemy niereg LZNK Ax-f
f=rand(m,1);
sol=A\f
g=U'*f; #||Ax-f||_2=|| SV'x-U'f||_2 -> 
y=zeros(n,1);
y(1)=g(1)/S(1,1);
x=V*y
norm(x-sol) #==0?
#oczywiscie za y(k) k>1 mozemy wziac cokolwiek
y(2:n)=randn(n-1,1);
A*V*y-A*x
#ale norma x najmniejsza
norm(x)
norm(V*y)

#zadanie wlasne 
disp("\n\nzadanie wlsne funkcja eig()");
A=[2 1;1 2]
eig(A)
[Q,D]=eig(A)
A-Q*D*Q'
A*Q-Q*D


#metoda potegowa - prosta implementacje (bez war stopu etc)
disp("\n\n testy metody potegowej");
N=2;
M=30;

function [r,x]=pot(A,M=30,x0=randn(length(A),1))
#implementacja metody  potegowej
#ze stala iloscia iteracji (w praktyce powinien byc war stopu)
 y=x=x0/norm(x0);
 for k=1:M,
  ax=A*x; 
  r=ax'*x;
  if(round(k/5)-(k/5)==0)
   printf("r(%d)=%10.7g\n",k,r);
  endif
  x=ax/norm(ax);
 endfor
end 
[r,x]=pot(A)
format long e
r-3

#powtorzmy dla -A
[r,x1]=pot(-A,30)
[r,x2]=pot(-A,31)
r+3
[x1,x2]



#to samo dla losowej A symetrycznej
N=10;
A=randn(N,N);
A=A+A';
lm=max(abs(eig(A))) #najwieksza ww co do modulu
x=randn(N,1);
x=x/norm(x);
[r,x]=pot(A,31)
r-lm  
r+lm  

#znow prosta macierz 
A=randn(3,3);A=A+A';
eig(A)
[r1,x1]=pot(A,71)
#a teraz startujemy z wektora prostopadlego do x
x0=randn(3,1); x0=x0-(x1'*x0)*x1;
x1'*x0
[r2,x2]=pot(A,60,x0)
r1-r2
#mimo startu z wektora prostopadlego znow zbieglismy do wektora wlasnego 
#x1
#a teraz co kilka iteracji reortogonalizujemy
[r2,x2]=pot(A,20,x0);
r2
#reortogonalizujemy
x0=x2;
x0=x0-(x1'*x0)*x1;
[r2,x2]=pot(A,20,x0)

disp("\n\n testy odwrotnej met potegowej");
function [r,x]=odwpot(A,s=0,M=30,x0=randn(length(A),1))
#implementacja metody odwrotnej potegowej
#ze stala iloscia iteracji (w praktyce powinien byc war stopu)
 y=x=x0/norm(x0);
 A=A-s*eye(size(A));
 [l,u,p]=lu(A);
 for k=1:M,
  ax=u\(l\(p*x));
  #ax=A\x; 
  rb=ax'*x;
  x=ax/norm(ax);
  r=1/rb+s;
 endfor
end 

r=odwpot(A,1.001,10)
r-1
r=odwpot(A,1.5,10)
r-1
#widac ze duzo gorzej

r=odwpot(A,3.001,10)
r-3
r=odwpot(A,2.5,10)
r-3
#a co gdy s posrodku?
r=odwpot(A,2,100)
#rozne wartosci...
#w praktyce szansa trafienia dokldnie posrodku minimalna
#ale blisko juz mozemy
r=odwpot(A,2.0001,100)
#za malo iteracji
r=odwpot(A,2.0001,10000)
r=odwpot(A,2.5,100000)
#no ale 10^5 iteracji troche duzo...

#metoda QR
disp("\n\ntesty metody QR dla zadania wlasnego");
function [D,Z]=qre(A,M=30,tol=1e-10)
 Z=eye(size(A));
 for k=1:M,
  [q,r]=qr(A);
  A=r*q;
  Z=Z*q;
 end
 erbd=norm(A-diag(diag(A)),1)
 if(erbd>tol) 
    printf("Warning: A_%d is far from being diagonal ||A-diag(A)||_1=%g\n",M,erbd);
 endif
 if(nargout<2)
  D=diag(A);
 else
  D=A;
 endif
endfunction
d=qre(A)
[D,Z]=qre(A)
norm(A*Z-Z*D,1)
#a teraz gorsza o prawie identycznych ww
A=[10 1;1 10];
eig(A)
d=qre(A)
[D,Z]=qre(A)
norm(A*Z-Z*D,1)
[D,Z]=qre(A,70)


