  clc
  printf("Please use the following functions:\n")
  printf("  ex4  - Compute the steepest descend method\n")
  printf("  ex4a - Compute the conjugate directions method\n")
  printf("  ex6  - Compute the CG method\n")
  printf("$")


  // Method: Steepest Descent
  function [approximate, correct, l2error] = SteepestDescent(A,b,x)
    
    // ************************************
    // * implementation of algorithm 10-1 *
    // ************************************
    
    // compute r (first step)
    r = b - A * x
    
    // compute p value (needed for second step)
    p = 0
    if norm(r) > 0
      p = r / norm(r)
    end

    // compute lambda (second step)
    lambda = (r' * p) / (p' * (A * p))
    
    // compute the approximation (second step)
    approximate = x + lambda * p
    
    // compute the correct solution
    correct = A \ b
    
    // compute the L2 norm
    l2error = norm(approximate - correct,2)
  endfunction


  // Method: Conjugate directions
  function [approximate, correct, l2error, rnew] = CD(A,b,x,r)
    
    // ************************************
    // * implementation of algorithm 10-2 *
    // ************************************
     
    // compute p value (needed for first step)
    p = 0 * r
    if norm(r) > 0
      p = r / norm(r)
    end

    // compute lambda (first step)
    lambda = (r' * p) / (p' * (A * p))
    
    // compute the approximation (second step)
    approximate = x + lambda * p
    
    // update r (third step)
    rnew = r - lambda * A * p
    
    // compute the correct solution
    correct = A \ b
    
    // compute the L2 norm
    l2error = norm(approximate - correct,2)
  endfunction


  function [approximate, correct, l2error, pnew, rnew, anew] = CG(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

    // compute the correct solution
    correct = A \ b
    
    // compute the L2 norm
    l2error = norm(approximate - correct,2) 
  endfunction


  function run(style)
    // given Maelstrom
    A = [2 0; 0 10]
    b = [0 0]'
    x = [4 sqrt(1.8) ]'
    iter=70

    // -------------------- Computational Code -------------------
    
    //plot2d(0,0,-1,"031"," ",[-5,-5,5,5])

    // conjugate directions
    if style == 2
      r = b - A*x
    end

    // conjugate gradient
    if style == 3
      r = b - A*x
      p = r
      a = r' * r
    end
    
    
    for k=1:iter

      // steepest descent
      if style == 1
        [x, cor, err] = SteepestDescent(A,b,x)
        //printf("Iteration %02d leads to an error of %.15f\n",k,err)
      else
        // conjugate directions
        if style == 2
          [x, cor, err,r] = CD(A,b,x,r)
          //printf("Iteration %02d leads to an error of %.15f (r updated to %f)\n",k,err,r)
        else
          // conjugate gradient
          if style == 3
            [x, cor, err, p, r, a] = CG(A,b,x,p,r,a)
            printf("Iteration %02d leads to an error of %f (r=%f, p=%f, a=%.20f)\n",k,err,r,p,a)
            if abs(a) < 0.00000000001
              printf("alpha is zero, returning\n")
              return
            end
          else
            printf("Method not known!\n")
            return
          end
        end
      end
    
      if modulo(k,5) == 0
        printf("Iteration %02d leads to an error of %.15f\n",k,err)
      end
      
      
      px = x(1)
      py = x(2)
      z=1
      
      //lxarc(px-0.5*err, py-0.5*err, err, err,0,360*64)
          
    end
  endfunction

  function ex4
    run(1)
  endfunction

  function ex4a
    run(2)
  endfunction

  function ex6
    run(3)
  endfunction
