

A=[1,1;1,0.99]
b=[2;1]
x1=A\b
A(2,2)-=0.01
x2=A\b

norm(x1-x2)
cond(A)


%macierz osobliwa
A=[1 1; 2 2]'
x=A\b
A*x-b

b=[3;3]
x=A\b
A*x-b
A*[1;1]-b

%inv()
A=[1, 2, 3;-1, 1, 1;1, 1, 1]
iA=inv(A)
norm(iA*A - eye(size(A)),1)
norm(A*iA-eye(size(A)),1)
s=randn(3,1);f=A*s;
x=iA*f;x-s

%sprawdzamy rozklady
A=rand(10,10);
[l,u,p]=lu(A);
norm(p*A-l*u,1)
[q,r]=qr(A);
[norm(q*r-A,1),norm(q'*q-eye(size(A)),1)]
B=A'*A;
R=chol(B);
norm(B-R'*R,1)

%czas
N=5;
M=12;
tim=zeros(M,6);
for k=1:M,
  A=rand(N,N);
  A=A'*A;
  b=rand(N,1);
tim(k,1)=N;
tic;inv(A)*b;tim(k,2)=toc;
tic;A\b;tim(k,3)=toc;
tic;[l,u,p]=lu(A);u\(l\(p*b));tim(k,4)=toc;
tic;[q,r]=qr(A);r\(q'*(b));tim(k,5)=toc;
tic;R=chol(A);R\(R'\b);tim(k,6)=toc;
N*=2;
end
tim
tim(2:end,:)./tim(1:M-1,:)

%rozwiazywanie ukladow z macierzami g-trojkatnymi
%U1 z samymi jedynkami a U2=-U1+2*I
N=10;
M=11;
er=zeros(M,5);
for k=1:M,
  U1=single(triu(ones(N,N)));
  U1=triu(ones(N,N));
  U2=2*eye(N)-U1;
  s=rand(N,1);
  f1=U1*s;f2=U2*s;
  x1=U1\f1;
  x2=U2\f2;
  er(k,:)=[N,norm(U1*x1-f1),norm(x1-s),norm(U2*x2-f2),norm(x2-s)];
  N+=5;
endfor
er
N=11;
U1=(triu(ones(N,N)));
inv(U1)
U2=2*eye(N)-U1;
inv(U2)



%rozwiazywanie ukladow z macierza Hilberta
N=1;
M=10;
er=zeros(M,5);
for k=1:M,
  H=single(hilb(N));
%  H=hilb(N);
  s=rand(N,1);
  f=H*s;
  x1=H\f;
  [Q,R]=qr(H);x2=R\(Q'*f);
  x3=invhilb(N)*f;
  er(k,1:3)=[N,norm(H*x1-f),norm(x1-s)];
  er(k,4:5)=[norm(H*x2-f),norm(x2-s)];
  er(k,6:7)=[norm(H*x3-f),norm(x3-s)];
  N++;
endfor
er




function  x1 = solvepowA(Q,R,b,c,a,f)

%macierz Schura
%A nieosobliwa znamy jej rozklad np QR
%A1 macierz powstala z A poprzez dodanie wiersza i kolumny
%A1=[A,b;c',a]
%rozwiaz uklad z A1 korzystajac z rozkladu A
%x1=[x;y] f=[f1;f2]
%Ax+b*y=f1 --> x+(A\b)y=x+w*y=(A\f1)=z czyli
%x=z-w*y
%c'*x+ay=g czyli podstawiajac x
%c'*(z-w*y)+ay=f2 czyli (a -c'*w)y=f2-c'*z 
%wyliczamy y i potem x=z-w*y
%
N=length(c);
f1=f(1:N);
f2=f(end);
%A1=[A,b;c',a]
w=R\(Q'*b); % Aw=b
z=R\(Q'*f1); % Az=f1
al=a - c'*w;
if(abs(al)>1e-14)
  y=(f2-c'*z)/al;
  x=z-w*y;
else
 disp("macierz numerycznie osobliwa?"); 
end
x1=[x;y];
end

N=2;
A=randn(N,N);
%A=eye(2,2);
[Q,R]=qr(A);
b=randn(N,1);c=randn(N,1);
%b=c=[1;2];
a=randn;
f=randn(N+1,1);
A1=[A,b;c',a];
s=A1*randn(3,1);f=A1*s;
x1 = solvepowA(Q,R,b,c,a,f)
norm(x1-s)


