

#zadania na rozklay uwarunkowanie etc

#zadanie porownaj czas rozwiazywania Ax=b poprzez \ lu qr chol
M=4;
N=10;
czas=zeros(M,5);
for k=1:M,
 czas(k,1)=N;
 A=rand(N,N);b=randn(N,1);
 C=A'*A; #C=C'>0
 tic;x1=A\b;czas(k,2)=toc;
 tic;R=chol(C);y=R\(R'\b);czas(k,3)=toc;
 tic;x=A\b;[l,u,p]=lu(A);x2=u\(l\(p*b));czas(k,4)=toc;
 tic;x=A\b;[q,r]=qr(A);x3=r\(q'*b);czas(k,5)=toc;
 N*=5;
endfor
disp("Czas rozwiazywania Ax=b poprzez rozne rozklady dla N losowej");
disp("Kol: N backslash lu cholesky qr (Householder)");
czas
#powtarzamy w pojedynczej precyzji
M=3;
N=10;
czasp=zeros(M,5);
for k=1:M,
 czasp(k,1)=N;
 A=single(rand(N,N));b=single(randn(N,1));
 C=A'*A; #C=C'>0
 tic;x1=A\b;czasp(k,2)=toc;
 tic;R=chol(C);y=R\(R'\b);czasp(k,3)=toc;
 tic;x=A\b;[l,u,p]=lu(A);x2=u\(l\(p*b));czasp(k,4)=toc;
 tic;x=A\b;[q,r]=qr(A);x3=r\(q'*b);czasp(k,5)=toc;
 N*=5;
endfor
disp("Czas rozwiazywania Ax=b poprzez rozne rozklady dla N losowej w pojed precyzji");
disp("Kol: N backslash lu cholesky qr (Householder)");
czasp 
disp("nie ma roznicy - pewnie backslash konwertuje do double....");



#zadanie porownaj czas rozwiazywania Ux=b poprzez \ i mnozenie przez m odwrotna
#i z \ ale z U w formacie rzadkim (sparse)
#U- macierz g-trojkatna z 2 diagonalna z 1 na diagonali i -1 na 
#
M=5;
N=10;
czas=zeros(M,4);
for k=1:M,
 czas(k,1)=N;
 U=eye(N)-diag(ones(N-1,1),1);b=randn(N,1);iU=inv(U);
 sU=sparse(U); #ta sama macierz format rzadki - tylko pameitamy o niezerowych 
#elementach macierzy
 tic;x1=U\b;czas(k,2)=toc;
 tic;x1=sU\b;czas(k,4)=toc;
 tic;x1=iU*b;czas(k,3)=toc;
 N*=5;
endfor
disp("Czas rozwiazywania Ux=b U g-tr 2diag 1 na diagonali i -1 na superdiag");
disp("kols: N backslash mnozenie przez odwrotna backslash dla U sparse");
czas
 



#zadanie porownaj rozwiazania Ukx=b poprzez \ 
#dla U1 macierzy g-tr z samym 1 - U2=U1-2*Id (1 nad diagonala, -1 na diagonali
#dla s znanego (losowego) policz bk=Uk*s xk=Uk\bk i normy xk-s i Axk-bk
#
M=14;N=15;
er=rer=zeros(M,5);
for k=1:M,
 er(k,1)=rer(k,1)=N;
 U1=ones(N,N);U1=triu(U1);U2=U1-2*eye(size(U1));
 s=randn(N,1);
 b1=U1*s;b2=U2*s;
 x1=U1\b1;er(k,2)=norm(x1-s);er(k,3)=norm(U1*x1-b1);
 x2=U2\b2;er(k,4)=norm(x2-s);er(k,5)=norm(U2*x2-b2);
 rer(k,[2,4])=er(k,[2,4])/norm(s); 
 rer(k,3)=er(k,3)/norm(b1);rer(k,5)=er(k,5)/norm(b2);
 N+=4;
endfor
format long e
disp("U1 macierz g-tr z 1 na diag i ponad; u2=U1-2*I proownanie bledow");
disp("x1=U1 bsl b1 x2=U2 bsl b2  bk=Uk*s s losowy");
disp("kols: N ||x1-s||_2 ||U1*x1-b1||_2 ||x2-s||_2 ||U2*x1-b2||_2");
eru=er
disp("Blad wzgledny tzn N ||x1-s||_2/||s||_2 etc");
reru=rer
#dla U1 wszystko OK ale dla U2 blad residualny (ostatnia kolumna ok ale 
#prawdziwy blad tragiczny dla N=57 blad rzedu 6.... 
format
semilogy(rer(:,1),rer(:,2),";blad U1;",rer(:,1),rer(:,4),";U2;");pause(1);

#zadanie porownaj rozwiazania Hx=b poprzez \ i mnozenie przez odwrotna
#obliczona funkcja invhilb H macierz Hilberta - hilb(N) ja generuje
#dla s znanego (losowego) policz bk=Uk*s xk=Uk\bk i normy xk-s i Axk-bk
#
M=9;N=5;
er=rer=zeros(M,5);
for k=1:M,
 er(k,1)=rer(k,1)=N;
 H=hilb(N);iH=invhilb(N);
 s=randn(N,1);
 b=H*s;
 x1=H\b;er(k,2)=norm(x1-s);er(k,3)=norm(H*x1-b);
 x2=iH*b;er(k,4)=norm(x2-s);er(k,5)=norm(H*x2-b);
 rer(k,[2,4])=er(k,[2,4])/norm(s); 
 rer(k,[3 5])=er(k,[3 5])/norm(b);
 N+=1;
endfor
format long e
disp("H macierz Hilberta iH jej odwrotna -  porownanie bledow algorytmu");
disp("rozwiazywania Hx=b : H bsl b vs iH*b");
disp("x1=H bsl b x2=iH*b b=H*s s losowy");
disp("kols: N ||x1-s||_2 ||H*x1-b||_2 ||x2-s||_2 ||H*x1-b||_2");
er
disp("Blad wzgledny tzn N ||x1-s||_2/||s||_2 etc");
rer
format
#a w pojedynczej precyzji?
M=5;N=3;
er=rer=zeros(M,5);
for k=1:M,
 er(k,1)=rer(k,1)=N;
 H=single(hilb(N));iH=single(invhilb(N));
 s=randn(N,1);s=single(s);
 b=H*s;
 x1=H\b;er(k,2)=norm(x1-s);er(k,3)=norm(H*x1-b);
 x2=iH*b;er(k,4)=norm(x2-s);er(k,5)=norm(H*x2-b);
 rer(k,[2,4])=er(k,[2,4])/norm(s); 
 rer(k,[3 5])=er(k,[3 5])/norm(b);
 N+=1;
endfor
format long e
disp("to samo w pojedynczej precyzji");
er
disp("blad wzgledny..");
rer
#juz dla N=6 warning
format


