A=[2 0; 0 10]
b=[0 0]'
x=[4 sqrt(1.8)]'


function x=graddes(A,b,x,iter)
t=linspace(0,2*%pi,60) 

I = diag([1 1])

  for i=1:iter


    x_out(i)=x(1)
    y_out(i)=x(2)
    
    
    r = b-A*x
    if (norm(r)==0) 
      l = 0
    else
      l = (r'*r)/((A*r)'*r)
    end
    x = (I-l*A)*x+l*b
  // fac = norm(I-l*A)*norm(inv(I-l*A))
   fac = max(spec(A))*max(spec(inv(A)))
    f=2*x(1)^2+10*x(2)^2 

linx=linspace(-sqrt(f/2),sqrt(f/2),100)
for i=1:length(linx)
  liny(i)=sqrt((f-2*linx(i)^2)/10)
end

plot2d(linx,liny)
plot2d(linx,-liny)
end
plot2d(x_out,y_out)
//printf("L2-Error: %f\n",norm(x-inv(A)*b))

endfunction


function x=CG(A,b,iter,x)
  timer()
  
  n=length(x)
  P(1:n,1)= b-A*x
  p       = b-A*x
  r       = b-A*x
  X(1:n,1)=x
  for i=1:iter
    l=(r'*p)/(p'*A*p)
    x=x+l*p
    r=r-l*A*p
    helper=0
    for j=1:i
      helper=helper+(r'*A*P(1:n,j))/(P(1:n,j)'*A*P(1:n,j))*P(1:n,j)
  end
    P(1:n,i+1)=r-helper
    p = r-helper
end
printf("L2-Error: %f\n",norm(x-inv(A)*b))
printf("time: %f\n",timer())
endfunction
    
    
function [A,f]=buildPoisson()
n=7
delta_x=1/8^2
  A(1,1)=2*delta_x
  A(2,1)=-1*delta_x
  
  for i=2:n-1
   A(i,i)=2*delta_x
   A(i+1,i)=-1*delta_x
   A(i-1,i)=-1*delta_x
  end
  
  A(n,n)=2*delta_x
  A(n-1,n)=-1*delta_x

  for i=1:n
    f(i)=2*i/8
    end
endfunction


function x = Jacobi(A,b,iter,x)
  timer();
  n = length(b)
  for k=1:n
    d(k)=A(k,k)
    end
  D = diag([d])
  M = inv(D)*(D-A)
  rho=max(spec(M))
  N = inv(D)
  for i=1:iter
    x=M*x+N*b
end
epsilon = abs(x-inv(A)*b)
//printf("Error is %f\n",max(epsilon))
printf("L2-Error: %f\n",norm(x-inv(A)*b))
printf("Time: %f\n",timer());
endfunction

function D = diago(A)
  [n,n]=size(A)
  for i=1:n
    d(i)=A(i,i)
  end
  D=diag([d])
endfunction

function A=rightify(A)
  [n,n]=size(A)
  for i=1:n
    for j=1:i
      A(i,j)=0
  end
end
endfunction

function A=leftify(A)
  [n,n]=size(A)
  for i=1:n
    for j=1:i
      A(j,i)=0
  end
end
endfunction

function x = SOR(A,b,iter,x)
  timer()
  n = length(b)
  
  D = diago(A)
  R = rightify(A)
  L = leftify(A)

  M1 = inv(D)*(D-A)
  rho=max(spec(M1))  
  
  omega = 2/(1+sqrt(1-rho^2))
  printf("Best omega:%f\n",omega)
  M = inv(D+omega*L)*((1-omega)*D-omega*R)
  N = omega*inv(D+omega*L)
//  printf("Maximal eigenvalue of the system matrix is %f\n",max(real(abs(spec(M)))))
  for i=1:iter
    x=M*x+N*b
end
epsilon = abs(x-inv(A)*b)
//printf("Error is %f\n",max(epsilon))
printf("L2-Error: %f\n",norm(x-inv(A)*b))
printf("time: %f\n",timer())
endfunction
