function soln = geSolve(A, b)
  numRows = size(A)(1);
  numCols = size(A)(2);
  if (numRows == numCols)
    ndims = numRows;

    scales = max(abs(A'))';
    aug = [A b];
    rowSeq = 1:ndims;
    numColsProcessed = 0;
    while (numColsProcessed < ndims - 1)
      jj = numColsProcessed + 1;
      ratioVec = abs(aug(rowSeq(jj:ndims), jj)) ./ scales(rowSeq(jj:ndims));
      offset = find(ratioVec == max(ratioVec))(1) - 1;
      if (offset > 0)
        saveVal = rowSeq(jj);
        rowSeq(jj) = rowSeq(jj + offset);
        rowSeq(jj + offset) = saveVal;
      end

      pivot = aug(rowSeq(jj), jj);
      for ii = rowSeq(jj+1:ndims)
        aug(ii, jj+1:ndims+1) = aug(ii, jj+1:ndims+1) - ...
               aug(ii,jj) / pivot * aug(rowSeq(jj), jj+1:ndims+1);
        aug(ii, jj) = 0;
      end
      numColsProcessed = numColsProcessed + 1;
    end

    A = aug(1:ndims, 1:ndims);
    b = aug(1:ndims, ndims + 1);
    lastxEqn = rowSeq(ndims);
    soln = [zeros(ndims-1, 1); b(lastxEqn)/A(lastxEqn, ndims)];
    for jj = (ndims - 1):(-1):1
      idx = rowSeq(jj);
      soln(jj) = (b(idx) - sum(A(idx, (jj+1):ndims)' .* soln(jj+1:ndims))) ...
                    / A(idx,jj);
    end

  else
    disp('Your matrix must be square.')
  end
end