function [] = prog1main()


x=[100 100.01 100.02 100.03 100.04];
y=[1 1.0001 1.0004 1.0009 1.0016]';
disp('Loesung aus Cholesky-Zerlegung')
disp(Aufgabe19(x,y))
disp('Loesung aus Householder-Spiegelungen')
disp(Aufgabe20a(x,y))
disp('Loesung aus Givens-Rotationen')
disp(Aufgabe20b(x,y))


 
 
function [a]=Aufgabe19(x,y)
  A=gegeben(x);
  b=A'*y;
  A=A'*A;
  a=Cholesky(A,b);
 
 function [a]=Aufgabe20a(x,y)
  A=gegeben(x);
  b=y;
  a=HS(A,b);
 
 function [a]=Aufgabe20b(x,y) 
  A=gegeben(x);
  b=y;
  a=GivensR(A,b);
  
 function [A] = gegeben(x)
  for i=1:length(x)
    A(i,1)=x(i)^(2);
    A(i,2)=x(i);
    A(i,3)=1;
  end
 
 function [x] = Cholesky(A,b)
  n=length(A);
  %Cholesky Zerlegung
  for k=1:n 
     sum1=0;
     for j=1:k-1
         sum1=sum1+L(k,j)^2;
     end    
     L(k,k)=(A(k,k)-sum1)^(1/2);
     for i=k+1:n
         sum2=0;
         for j=1:k-1
         sum2=sum2+L(i,j)*L(k,j);
         end  
           L(i,k)=(A(i,k)-sum2)/L(k,k);
     end  
  end
  %vorwaert
  for i=1:n 
      sum=0;
      for j=1:i-1
          sum=sum+L(i,j)*z(j);
      end
      z(i)=(b(i)-sum)/L(i,i);
  end    
  L=L';
  %zurueckwaert
  for i=n:-1:1 
      sum=0;
      for j=i+1:n
          sum=sum+L(i,j)*x(j);
      end
      x(i)=(z(i)-sum)/L(i,i);
  end
  
 function [x] = GivensR(A,b)
 m=5;
 n=3;
  %Givens-Rotationen
  for i=1:(n-1)
   for j=i+1:m   
    r=sqrt(A(i,i)^2+A(j,i)^2);
    c=A(i,i)/r;
    s=A(j,i)/r;
    
    t=b(i);
    b(i)=c*t+s*b(j);
    b(j)=-s*t+c*b(j);
    
    for k=i:n
        p=A(i,k);
        q=A(j,k);
        A(i,k)=c*p+s*q;
        A(j,k)=-s*p+c*q;
    end
   end 
  end 
  %zurueckwaert
  x(n)=b(n)/A(n,n);
   for i=(n-1):-1:1
      sum=A(i,i+1:n)*x(i+1:n)';
      x(i)=(b(i)-sum)/A(i,i);
   end 

 function [x] = HS(A,b)
   m=5;
   n=3;
   %Householder-Spiegelungen
   for i=1:n
       Q(1:m,1:m)=0;
     for k=1:m
       Q(k,k)=1;
     end
     I=Q(1:m-i+1,1:m-i+1);
     a=0;
     for j=i:m
         a=A(j,i)^2+a;
     end
     a=sqrt(a);
     v=A(i:m,i);
     v(1)=v(1)-a;
     p=v*v';
     q=v'*v;
     QQ=I-2*p/q;
     Q(i:m,i:m)=QQ;
     A=Q*A;
     b=Q*b;
   end
   %zurueckwaert
   x(n)=b(n)/A(n,n);
   for i=(n-1):-1:1
      sum=A(i,i+1:n)*x(i+1:n)';
      x(i)=(b(i)-sum)/A(i,i);
   end 
   
   