function [] = prog1main()

kk=[3,5,9,17];

 for it=1:4 
   k=kk(it);
   n=2*(k^2)-k;
   h=1/(k-1);
   hh=h^2;
   
   X=XX(k,n,h);
   for I=1:2*k^2-k
        U(I)=sin(pi*X(I,1))+cos(pi*X(I,2));
   end  
   for I=1:k+1
       A(I,I)=1;
       A(n+1-I,n+1-I)=1;
       B(I)=g(X(I,1),X(I,2));
       B(n+1-I)=g(X(n+1-I,1),X(n+1-I,2));
   end
   for I=k+2:2*k-1
       A(I,I)=4/hh;
       A(I,I-1)=-1/hh;
       A(I,I+1)=-1/hh;
       A(I,I-k)=-1/hh;
       A(I,I+k)=-1/hh;
       B(I)=f(X(I,1),X(I,2));
   end   
   i=2*k-1;
   for I=1:2*(k-2)
       i=i+1;
       A(i,i)=1;
       B(i)=g(X(i,1),X(i,2));
       i=i+1;
       A(i,i)=1;
       B(i)=g(X(i,1),X(i,2));
       for J=1:k-2
           A(i+J,i+J)=4/hh;
           A(i+J,i+J-1)=-1/hh;
           A(i+J,i+J+1)=-1/hh;
           A(i+J,i+J-k)=-1/hh;
           A(i+J,i+J+k)=-1/hh;
           B(i+J)=f(X(i+J,1),X(i+J,2));
       end
       i=i+(k-2);
   end
   
      
   B=B';
   r=B-A*B;
   u=B;
   rr=r+88;
   p=r;
   Iter=0;
   w=sGS(A);
   r1=r;
   u1=u;
   p1=w*r1;
   ro=p1'*r1;
   Iter1=0;
   while Fehlernorm(r,rr)>10^(-7)
      Ap=A*p;
      a=(r'*p)/(p'*Ap);
      u=u+a*p;
      rr=r;
      r=r-a*Ap;
      be=(Ap'*r)/(p'*Ap);
      p=r-be*p;
      Iter=Iter+1;
    end
   
   while r1'*(w*r1)>10^(-7)
     aa=A*p1;
     alfa=ro/(aa'*p1);
     u1=u1+alfa*p1;
     r1=r1-alfa*aa;
     q=w*r1;
     roo=q'*r1;
     p1=q+roo/ro*p1; 
     ro=roo;
     Iter1=Iter1+1;
   end  
u=u';
B=B';
k
disp('CG , Anzahl der Iteration')
disp(Iter)
disp('Feler')
disp(Fehlernorm(U,u))
disp('SGS-PCG, Anzahl der Iteration')
disp(Iter1)
disp('Feler zur exakten Loesung')
disp(Fehlernorm(U,u1))
% disp(U)
% disp(u)
%disp(B)  
%disp(A)
A=0;
B=0;
U=0;
u=0;
r=0;
i=0;
end


function X=XX(k,n,h)
  for I=1:n
       i= fix((I-1)/k);
       j= mod((I-1),k) ;
       X(I,1)=h*i;
       X(I,2)=h*j;
  end

function f=f(x,y)
      f=(pi^2)*(sin(pi*x)+cos(pi*y));

function g=g(x,y)    
      g=sin(pi*x)+cos(pi*y);

function w = sGS(A)
   n=length(A);
   for i=1:n
        for j=1:i
         link(i,j)=A(i,j);
        end
        dia(i,i)=1/A(i,i);
        for j=i:n
         recht(i,j)=A(i,j);
        end
    end
    w=link*dia*recht;
    w=w^(-1);
    
function w = Jaco(A)
   n=length(A);
   for i=1:n
       w(i,i)=1/A(i,i);
   end
    
function normm=Fehlernorm(x,y)
     n=length(x);
     normm=0;
     for i=1:n
         normm=normm+(x(i)-y(i))^2;
     end
     normm=sqrt(normm);
     


