clearglobal
clc
printf("Please use the following functions:\n")
printf("  exercise2b - Condition numbers\n")
printf("  exercise2c - CG-Errors\n")
printf("  exercise2c - preconditioned CG-Errors\n")
printf("$")

function [H,h] = Hilbert(n)
  H = zeros(n,n);
  h = zeros(n,1)
  for i=1:n
    for j=1:n
      H(i,j) = 1 / ( i+j-1);
      h(i) = h(i) + H(i,j);
    end
  end
endfunction


function exercise2b
  a = [2 3 5 10 15 25 50 100]'
  
  for i=1:size(a,1)
    
    [H,h] = Hilbert(a(i))
    
    // sort eigenvalues descending order
    //s = sort(spec(H))

    // condition number
    // largest / smallest
    // c = s(1) / s(a(i))
    c = cond(H)
    
    //printf("%02d has %20.f   %.20f\n",a(i),c,s(a(i)))
    printf("%3d -> %25.f\n",a(i),c)
  end
endfunction


function [approximate, pnew, rnew, anew] = CGstep(A,b,x,p,r,a)

    v = A * p
    lambda = a / (v' * p)
    
    approximate = x + lambda * p
    rnew = r - lambda * v
    anew = rnew' * rnew
    pnew = rnew + anew / a * p

endfunction

// conjugate gradient with initialisation vector
function x = CGx(A,b,x)
  r = b - A*x
  p = r
  a = r' * r
  for k=1:size(A,1)
    v = A * p
    lambda = a / (v' * p)
    x = x + lambda * p
    r = r - lambda * v
    anew = r' * r
    p = r + anew / a * p
    a = anew
  end
endfunction

// preconditioned conjugate gradient
// with initialisation vector
function x = CGxPre(A,b,x,P)
  //P = Pl' * Pl

  r = b - A*x
  p = P * r
  a = r' * p
  for k=1:size(A,1)
    v = A * p
    lambda = a / (v' * p)
    x = x + lambda * p
    r = r - lambda * v
    //if norm(r,'inf') < 0.00000001
    //  return
    //end
    z = P * r
    anew = r' * z
    p = z + anew / a * p
    a = anew
  end
endfunction

// conjugate gradient with initialisation vector 0
function x = CG(A,b)
  x = zeros(size(A,1),1)
  x = CGx(A,b,x)
endfunction

// preconditioned conjugate gradient
// with initialisation vector 0
function x = CGPre(A,b,P)
  x = zeros(size(A,1),1)
  x = CGxPre(A,b,x,P)
endfunction

function exercise2c
  a = [2 3 5 10 15 25 50 100]'

  //
  n=30
  //a = ones(n,1)
  a(1) = 2
  //for i=2:n
  //  a(i) = a(i) + a(i-1)
  //end
  
  
  for i=1:size(a,1)
    
    [H,h] = Hilbert(a(i))
    
    x = CG(H,h)
    cor = ones(size(H,1),1)
    n = norm(cor-x,'inf')       

    printf("%3d -> %1.20f\n",a(i),n)
  end  
endfunction


// Does not assume that A spd
function [P1, P2 ] = Precondition(A)
  D = diag(diag(A))
  L = tril(A) - D
  U = triu(A) - D
  
  du = inv(D+U)
  dl = inv(D+L)
  
  P1 = du * L * dl * U
  P2 = du - du * L * dl

endfunction

function exercise2d
  a = [2 3 5 10 15 25 50 100]'

  for i=1:size(a,1)
    
    [H,h] = Hilbert(a(i))
    
    xo = CG(H,h)
    cor = ones(size(H,1),1)
    no = norm(cor-xo,'inf')

    [ P1, P2 ] = Precondition(H)
    
    // Jacobi-Iteration
    // P = D^{-1}
    //xn = CGPre(H,h,inv(diag(diag(H))))
    xn = CGPre(H,h,P2)
    nn = norm(cor-xn,'inf')

    printf("%3d -> %1.20f\n",a(i),nn)
  end  
endfunction

function exercise2d2
  a = [2 3 5 10 15 25 50 100]'

  //
  n=30
  a = ones(n,1)
  a(1) = 2
  for i=2:n
    a(i) = a(i) + a(i-1)
  end
  
  
  for i=1:size(a,1)
    
    [H,h] = Hilbert(a(i))
    
    xo = CG(H,h)
    cor = ones(size(H,1),1)
    no = norm(cor-xo,'inf')

    [ P1, P2 ] = Precondition(H)
    
    P = 0.001*eye(size(H,1),size(H,1))
    // Jacobi-Iteration
    // P = D^{-1}
    //xn = CGPre(H,h,inv(diag(diag(H))))
    xn = CGPre(H,h,P)
    nn = norm(cor-xn,'inf')


    printf("%3d -> %25.f %1.10f\n",a(i), cond(H),   no)
    printf("%3d -> %25.f %1.10f\n",a(i), cond(P*H), nn)
    printf("-------------------------------------------------------------------\n")
  end  
endfunction
