// To compute the error measurements,
// use the Ex[1-4] functions
// * Ex1 - trivial iteration
// * Ex2 - Jacobi iteration
// * Ex3 - Gauss-Seidel iteration
// * Ex4 - SOR method (converges to wrong value)


// linear system of equations
function [A,b,x] = LSE()
  A = [ 0.7 -0.4; -0.2 0.5]
  b = [ 0.3 0.3]'
  x = [ 21 -19]'
endfunction

// maximum norm of a vector
function n = mynorm (a)
  n = max( abs(a))
endfunction

// trivial iteration scheme
function ret = trivial(A,b,x,n)
  s = size(A,1)
  Anew = eye(s,s) - A
  for i=0:n-1
    x = Anew * x + b;
  end
  ret = x
endfunction


// Jacobi Iteration scheme
function ret = jacobi(A,b,x,n)
  // diagonal elements
  D = diag(diag(A))
  
  // inverse of diagonal elements
  invD = inv(D)
  
  // iteration matrix M
  M = invD * (D - A)
  
  for i=0:n-1
    x = M * x + invD * b;
  end
  ret = x
endfunction


// Gauss-Seidel Iteration scheme
function ret = gaussseidel(A,b,x,n)
  // diagonal elements
  D = diag(diag(A))
  
  // lower and upper triangle without diagonal
  L = tril(A) - D
  U = triu(A) - D
  

  // iteration matrices N and M
  N = inv(D+L)
  M = -N * U
  
  for i=0:n-1
    x = M * x + N * b;
  end
  ret = x
endfunction



// SOR Iteration scheme
function ret = SOR(A,b,x,n)
  // diagonal elements
  D = diag(diag(A))
  
  // lower and upper triangle without diagonal
  L = tril(A) - D
  U = triu(A) - D

  
  // matrix for the computation of the
  // optimal relaxation parameter
  M = inv(D) * (D-A)
  
  // compute maximum absolute eigenvalue
  p = mynorm(spec(M))
  
  // relaxation weight
  w = 2 / ( 1 + sqrt( 1-p*p))
  
  M = inv( D + w*L ) * [(1-w) * D - w*U]
  N = w * inv( D + w*L )
  
  //printf("wopt=%f\n",wopt)
  
  for i=0:n-1
    x = M * x + N * b;
  end
  ret = x
endfunction















// Trivial Iteration
function Ex1()
  n=7
  [A,b,x] = LSE()
  xreal = A \ b
  disp(xreal)
  for i=0:n
    iter = i * 10;
    ret = trivial(A,b,x,iter)
    nor = mynorm(xreal - ret)
    printf("\n%d iterations lead to an absolute error of %f\n",iter, nor);
  end
endfunction

// Jacobi Iteration
function Ex2()
  n=7
  [A,b,x] = LSE()
  xreal = A \ b
  disp(xreal)
  for i=0:n
    iter = i * 5;
    ret = jacobi(A,b,x,iter)
    nor = mynorm(xreal - ret)
    printf("\n%d iterations lead to an absolute error of %f\n",iter, nor);
    disp(ret)
  end
endfunction

// Gauss-Seidel Iteration
function Ex3()
  n=5
  [A,b,x] = LSE()
  xreal = A \ b
  disp(xreal)
  for i=0:n
    iter = i * 5;
    ret = gaussseidel(A,b,x,iter)
    nor = mynorm(xreal - ret)
    printf("\n%d iterations lead to an absolute error of %f\n",iter, nor);
    disp(ret)
  end
endfunction


// SOR Iteration
function Ex4()
  n=3
  [A,b,x] = LSE()
  xreal = A \ b
  disp(xreal)
  for i=0:n
    iter = i*5;
    ret = SOR(A,b,x,iter)
    nor = mynorm(xreal - ret)
    printf("\n%d iterations lead to an absolute error of %f\n",iter, nor);
    disp(ret)
  end
endfunction

