function p4
  gauss_test(@gauss_vanilla, 'no piv.');
  gauss_test(@gauss_partial_direct, 'partial piv. (swap)');
  gauss_test(@gauss_partial_indirect, 'partial piv. (indirect)');
  gauss_test(@gauss_scaled_partial, 'scaled partial piv.');
  %gauss_test(@backslash, 'matlab-backslash');
end
% --------------------------------------------------------------
function sol = gauss_test(linsolve, description)
  mat = [
  -9,11,-21,63,-252;
  70,-69,141,-421,1684;
  -575,575,-1149,3451,-13801;
  3891,-3891,7782,-23345,93365;
  1024,-1024,2048,-6144,24572];

  rhs = [-356;2385;-19551;132274;34812];
  true_sol = [1;-1;1;-1;1];
  assert(mat*true_sol == rhs)

  sol = linsolve(mat, rhs);

  residual = mat*sol - rhs;
  err = sol-true_sol;
  disp(sprintf('%-25s - |A*x-b| = %15g - |x-x_true| = %15g', ...
    description,...
    sqrt(residual'*residual), ...
    sqrt(err'*err))...
    )
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
