function A=hilbert1(n)
  for i=1:n
    for j=1:n
      A(i,j)=1/(i+j-1)
    end
  end
endfunction

function b=hilb(n)
  a=hilbert1(n)
  b(1:n)=0
  for i=1:n
    for j=1:n
      b(i)=b(i)+a(i,j)
    end
  end
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 = GaussSeidel(A,b,iter,x)
  n = length(b)
  D = diago(A)
  R = rightify(A)
  L = leftify(A)

  
  M = -inv(D+L)*R
  N = inv(D+L)
  for i=1:iter
    x=M*x+N*b
end
epsilon = abs(x-inv(A)*b)
epsilon1 = abs(x-[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]')
printf("Error is %f %f\n",max(epsilon), max(epsilon1))
endfunction


function x=CG1(A,b,iter,x)
    timer()
 n=length(x)
 p=b-A*x
 r=b-A*x
 a_new = r'*r
 a_old = r'*r
 for i=1:iter
  a_old = a_new
    if (a_new == 0)
      break
  end
    v=A*p
    l = a_old/(v'*p)
    
    x=x+l*p
    r=r-l*v
    a_new=(r'*r)
    p=r+a_new/a_old*p

end
sol(1:n)=1
printf("Numerical L2-Error: %f\n",max(abs(x-inv(A)*b)))
printf("Real L2-Error: %f\n",max(abs(x-sol)))
printf("time: %f\n",timer())
endfunction


function x=PCG(A,b,iter,x)
    timer()
 n=length(x)
  D = diago(A)
  R = rightify(A)
  L = leftify(A)
  P= inv(D+R)*D*inv(D+L)
 // P=inv(D)
//  P=chol(P)'*chol(P)
  //PL=chol(P)*A*chol(P)'
    q=cond(P)
 // vgl = max(spec(PL))/min(spec(PL))
  printf("Konditionszahl P %f \n",q)
    q=cond(A)
 // vgl = max(spec(A))/min(spec(A))
  printf("Konditionszahl A %f \n",q)
//  P=inv(D)
  r=b-A*x
  p=P*r
  a_old =r'*p
  a_new =r'*p
 
 
 for i=1:iter
  a_old = a_new
    if (a_new ==0)
      break
  end
    v=A*p
    l = a_old/(v'*p)
    
    x=x+l*p
    r=r-l*v
    z=P*r
    a_new=(r'*z)
    
    p=z+a_new/a_old*p

end
sol(1:n)=1
printf("Numerical L2-Error: %f\n",max(abs(x-inv(A)*b)))
printf("Real L2-Error: %f\n",max(abs(x-sol)))
printf("time: %f\n",timer())
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



function u=exercise2(n)
  A=hilbert1(n)
  b=hilb(n)
  x=cond(A)
  vgl = max(spec(A))/min(spec(A))
  printf("Konditionszahl %f %f\n",x,vgl)
  u=A\b
endfunction

function u=exercise3(n,iter)
  A=hilbert1(n)
  b=hilb(n)
  z(1:n)=0
  u=CG1(A,b,iter,z)
endfunction

function u=exercise4(n,iter)
  A=hilbert1(n)
  b=hilb(n)
  z(1:n)=0
  u=PCG(A,b,iter,z)
endfunction
