function p3
  values = [5,6,7,11,12,13];

  for n=1:length(values)
    run_all_tests(values(n));
  end
end
% --------------------------------------------------------------
function run_all_tests(n)
  disp('')

  % The matrix in the problem is called a 'Hilbert matrix'.
  % Depending on the size of n, it has an exponentially big
  % condition number, which means we should expect increasingly
  % bad solutions.

  fprintf('n=%d, cond(hilb(%d))=%g\n', n, n, cond(hilb(n)))
  disp('')
  gauss_test(n, @gauss_vanilla, 'no piv.');
  gauss_test(n, @gauss_partial_direct, 'partial piv. (swap)');
  gauss_test(n, @gauss_partial_indirect, 'partial piv. (indirect)');
  gauss_test(n, @gauss_scaled_partial, 'scaled partial piv.');
  %gauss_test(n, @backslash, 'matlab-backslash');
end
% --------------------------------------------------------------
function sol = gauss_test(n, linsolve, description)

  mat = hilb(n);

  %mat = [5,7,9;1,2,4;9,3,19];
  %n = 3;

  real_sol = ones(n,1);
  rhs = mat*real_sol;
  sol = linsolve(mat, rhs);

  residual = mat*sol - rhs;
  err = sol-real_sol;
  fprintf('%-25s - |A*x-b| = %15g - rel.error [%%] = %15g\n', ...
    description,...
    sqrt(residual'*residual), ...
    100 * sqrt(err'*err)/sqrt(real_sol'*real_sol))
end
% --------------------------------------------------------------
function sol = backslash(A, b)
  sol = A \ b;
end
% --------------------------------------------------------------
function sol = gauss_vanilla(A, b)
  [n, dummy] = size(A);

  for step = 1:n-1
    pivot = A(step, step);
    for i = step+1:n
      factor = A(i,step)/pivot;
      A(i,step:n) = A(i,step:n) - factor*A(step,step:n);
      b(i) = b(i) - factor*b(step);
    end
  end

  sol = b;
  for step = 0:n-1
    row = n-step;
    sol(row) = sol(row) - A(row,(row+1):n)*sol(row+1:n);
    sol(row) = sol(row) / A(row,row);
  end
end
% --------------------------------------------------------------
function sol = gauss_partial_direct(A, b)
  % A version of Gauss elimination with partial pivoting
  % using direct row swapping.

  [n, dummy] = size(A);

  for step = 1:n-1
    % find argmax_{step <= piv_row <= n} |a(piv_row, step)|
    max_row = 0;
    max_value = 0;

    for piv_row = step:n
      if (abs(A(piv_row, step))>max_value) 
        max_value = abs(A(piv_row, step));
        max_row = piv_row;
      end
    end

    % swap rows
    if (step ~= max_row)
      temp = A(max_row, step:n);
      A(max_row, step:n) = A(step, step:n);
      A(step, step:n) = temp;

      temp = b(max_row);
      b(max_row) = b(step);
      b(step) = temp;
    end

    pivot = A(step,step);
    for i = step+1:n
      factor = A(i,step)/pivot;
      A(i,step:n) = A(i,step:n) - factor*A(step,step:n);
      b(i) = b(i) - factor*b(step);
    end
  end

  sol = b;
  for step = 0:n-1
    row = n-step;
    sol(row) = sol(row) - A(row,(row+1):n)*sol(row+1:n);
    sol(row) = sol(row) / A(row,row);
  end
end
% --------------------------------------------------------------
function sol = gauss_partial_indirect(A, b)
  % A version of Gauss elimination with partial pivoting
  % using indirect access.

  [n, dummy] = size(A);

  permut = 1:n;

  for step = 1:n-1
    % find argmax_{step <= piv_row <= n} |a(piv_row, step)|
    max_row = 0;
    max_value = 0;

    for piv_row = step:n
      if (abs(A(permut(piv_row), step))>max_value) 
        max_value = abs(A(permut(piv_row), step));
        max_row = piv_row;
      end
    end

    % save values for postcond. check
    prev_idx = permut(step);
    target_idx = permut(max_row);

    % swap out rows in permutation
    if (max_row ~= step)
      temp = permut(max_row);
      permut(max_row) = permut(step);
      permut(step) = temp;
    end

    pms = permut(step);
    pivot = A(pms, step);

    % verify postconditions
    assert(pms == target_idx)
    assert(permut(max_row) == prev_idx)

    % subtract out pivot row
    for i = step+1:n
      pmi = permut(i);
      factor = A(pmi,step)/pivot;
      A(pmi,step:n) = A(pmi, step:n) - factor*A(pms, step:n);
      b(pmi) = b(pmi) - factor*b(pms);
    end
  end

  % this prints a triangular matrix, if enabled:
  if (0)
    for i=1:n
      A(permut(i),1:n)
    end
  end

  % perform backsubstitution
  sol = zeros(n,1);
  for step = 0:n-1
    row = n-step;
    prow = permut(row);
    
    sol(row) = b(prow);

    for j = row+1:n
      sol(row) = sol(row) - A(prow,j)*sol(j);
    end
    sol(row) = sol(row) / A(prow,row);
  end
end
% --------------------------------------------------------------
function sol = gauss_scaled_partial(A, b)
  % A version of Gauss elimination with scaled partial pivoting.

  [n, dummy] = size(A);

  permut = 1:n;

  % build row maxima table
  rowmaxes = zeros(n,1);
  for i=1:n
    rowmaxes(i) = max(abs(A(i,1:n)));
  end

  for step = 1:n-1
    % find argmax_{step <= piv_row <= n} |a(piv_row, step)|/rowmaxes(piv_row)
    max_row = 0;
    max_value = 0;

    for piv_row = step:n
      value = abs(A(permut(piv_row),step))/rowmaxes(permut(piv_row));
      if (value>max_value) 
        max_value = value;
        max_row = piv_row;
      end
    end

    % save values for postcond. check
    prev_idx = permut(step);
    target_idx = permut(max_row);

    % swap out rows in permutation
    if (max_row ~= step)
      temp = permut(max_row);
      permut(max_row) = permut(step);
      permut(step) = temp;
    end

    pms = permut(step);
    pivot = A(pms, step);

    % verify postconditions
    assert(pms == target_idx)
    assert(permut(max_row) == prev_idx)

    % subtract out pivot row
    for i = step+1:n
      pmi = permut(i);
      factor = A(pmi,step)/pivot;
      A(pmi,step:n) = A(pmi, step:n) - factor*A(pms, step:n);
      b(pmi) = b(pmi) - factor*b(pms);
    end
  end

  % this prints a triangular matrix, if enabled:
  if (0)
    for i=1:n
      A(permut(i),1:n)
    end
  end

  % perform backsubstitution
  sol = zeros(n,1);
  for step = 0:n-1
    row = n-step;
    prow = permut(row);
    
    sol(row) = b(prow);

    for j = row+1:n
      sol(row) = sol(row) - A(prow,j)*sol(j);
    end
    sol(row) = sol(row) / A(prow,row);
  end
end
% --------------------------------------------------------------
function assert(m)
  if (m == 0)
    error('assertion failed')
  end
end
