#!/usr/bin/octave
%Mn 16/17 lab na lu
%zad 1-3 trywialne

%zad 4
%bledy dla 2 m grnotr U1=same 1 na diag i ponadiag
% U2=1 na diagonali ; -1 nad diagonala
%format long e
nn=[6 16,23,30,52,60,100];
err=zeros(length(nn),7);
for k=1:length(nn),
 n=nn(k);
 U1=triu(ones(n,n));
 U2=-U1+2*eye(n,n);
 if(n<=10)
  iU1=inv(U1)#2 diagonalna
  iU2=inv(U2)#pelna g-tr w prawym rogu spora liczba....
 endif
 s=rand(n,1);ns=norm(s);
 f1=U1*s;nf1=norm(f1);
 f2=U2*s;nf2=norm(f2);
 x1=U1\f1;
 x2=U2\f2;
 e1=norm(x1-s)/ns;
 r1=norm(U1*x1-f1)/nf1;
 c1=norm(U1,1)*norm(inv(U1),1);
 e2=norm(x2-s)/ns;
 r2=norm(U2*x2-f2)/nf2;
 c2=norm(U2,1)*norm(inv(U2),1);
 err(k,:)=[n r1 e1 c1  r2 e2 c2] ;
endfor


#zad 5 jak 4 ale dla m Hilberta Hn=(1/(i+j-1))_{i,j=1}^n
printf("testy dla m Hilberta\n");
format long e
nn=[4,16,23,30,52,60];

for k=1:length(nn),
 n=nn(k)
 H=hilb(n);
 s=rand(n,1);ns=norm(s);
 f1=H*s;nf1=norm(f1);
 x1=H\f1;
 e1=norm(x1-s)/ns
 r1=norm(H*x1-f1)/nf1
endfor

%to samo ale rozwiazujac uklad mnozac przez macierz odwrotna
%uzyskana przez inhilb()
printf("testy dla m Hilberta - uzywajc inhilb\n");
format long e
nn=[4,10,11,12,13,14,15,16,23,30,52];
printf("[n] ||H*invhilb(n)I|| |H*iH-I|| err res \n"); 
for k=1:length(nn),
 n=nn(k);
 H=hilb(n);
 s=rand(n,1);ns=norm(s);
 f2=H*s;nf2=norm(f2);
 iH=inv(H);inH=invhilb(n);
 erm=norm(H*inH-eye(n),1);
 erm1=norm(H*iH-eye(n),1); 
 x2=invhilb(n)*f2;
 e2=norm(x2-s)/ns;
 r2=norm(H*x2-f2)/nf2;
 printf("[%2d] %7.5g %7.5g \t \t %7.5g \t%7.5g\n",n,erm,erm1,e2,r2); 
endfor



#zad 8
#LU dla m 3diagonalnej bez permutacji wierszy

function x=solve3diag(a,b,c,f)
#Funkcja rozwiazuje ukald AX=f (o ile rozklad LU A istnieje)
#Wejscie:
#macierz A 3-diagonalna z wektorem a jako gl diagonala
#    b- poddiagonala a c naddiagonala
#f- wektor prawej strony
#Wyjscie:
#x- rozwiazanie Ax=f

n=length(a);
x=zeros(n,1);
#rozklad elementy l(k) w b a U w a (c sie nie zmienia)
for k=1:n-1,
 b(k)=b(k)/a(k);
 a(k+1)=a(k+1)-b(k)*c(k);
 f(k+1)=f(k+1)-b(k)*f(k);
endfor

x(n)=f(n)/a(n);
for k=(n-1):-1:1,
 x(k)=(f(k)-c(k)*x(k+1))/a(k);
endfor

endfunction

printf("\n testy dla A 3-diagonalnej \n");

nn=[10,100,1000,1000];

for k=1:length(nn),
 n=nn(k)
#a=2*ones(n,1);b=c=-ones(n-1,1);
a=rand(n,1);b=c=rand(n-1,1);
A=diag(a)+diag(b,-1)+diag(c,1);
s=rand(n,1);f=A*s;
tic;y=A\f;t1=toc
e1=norm(y-s)
r1=norm(A*y-f)
tic;x=solve3diag(a,b,c,f);t2=toc
e2=norm(x-s)
r2=norm(A*x-f)
endfor

function [x1,res]=solveAp1(q,r,b,c,al,f1)
%rozwiazujemy uklad rownan
%dla macierzy A1=[A b;c' al] tzn
%A1**x1=f1 znajac rozklad QR macierzy A 
%jak kosztem O(n^2) rozwiazac go?
%x1=[x;xnp] to
%x=inv(A)*f-inv(A)*b*xnp
%(al - c'*inv(A)*b)xnp=f-c'*inv(A)*f --->xnp
%podstawiamy do x 
%INPUT
%q,r - QR factors of A
%b,c - vectors nx1 
%al - skalr na pozycji (n+1)x(n+1)
%f1 - vector (n+1)x1 rhs of system of equations A1*x1=f1
%OUTPUT
%x1 solution
%res = ||A1*x1-f1||_2

f=f1(1:end-1);
iAb=r\(q'*b);
iAf=r\(q'*f);
beta=al-c'*iAb;
if(abs(beta)<1e-12)
 printf("A1 may be  singular or badly conditioned\n");
else
 xnp=(f1(end)-c'*iAf)/beta;
end
x=iAf-iAb*xnp;
x1=[x;xnp];
if(nargout>1)
  res=norm([q*r*x+b*xnp-f;c'*x+al*xnp-f1(end)]);
end%if
end%function
disp("test met rozw. ukladu z [A,b;c',al] znajac rozklad QR macierzy A\n");
n=10;
A=rand(n,n);
b=randn(n,1);c=randn(n,1);al=1;
[q,r]=qr(A);
A1=[A,b;c',al];
f1=randn(n+1,1);
[x1,res]=solveAp1(q,r,b,c,al,f1);
norm(A1*x1-f1)
%porownajmy czas...
n=10;M=5;t1=t2=zeros(M,1);
for k=1:M,
A=rand(n);
b=randn(n,1);c=randn(n,1);al=1;
[q,r]=qr(A);
A1=[A,b;c',al];
f1=randn(n+1,1);
tic;solveAp1(q,r,b,c,al,f1);t1(k)=toc;
tic;A1\f1;t2(k)=toc;
printf("[%d] naszalg=%g backslash=%g \n",n,t1(k),t2(k));
n*=2;
end%for
%