
function [] = progmain()
mode(0);

ieee(1);

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));// die exakte Loesung
   end  
   //bestimme A und B
   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
   
   
   r=B-A*B;
   u=B;
   rr=r+88;// hier 88 ist egar 
   p=r;
   Iter=0;
   w=sGS(A);//symmetrische Gauss-Seidel Iteration als Vorkondition
   r1=r;
   u1=u;
   p1=w*r1;
   ro=p1'*r1;
   Iter1=0;
   // eine numerische loesung u durch   CG Verfahren
   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
   
   //eine numerische loesung u1 durch   SGS-PCG Verfahren
   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';
disp('k=')
disp(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(A)

A=0;
B=0;
U=0;
u=0;
r=0;
i=0;
end

endfunction 

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

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

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

function w = sGS(A)
   n=size(A,'r');
   for i=1:n
        for j=1:i
         lk(i,j)=A(i,j);
        end
        dia(i,i)=1/A(i,i);
        for j=i:n
         rt(i,j)=A(i,j);
        end
    end
    w=lk*dia*rt;
    w=w^(-1);
endfunction

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);
 endfunction
