matrix rref( matrix A)
/*
   reduced row echelon form, with partial pivoting
*/
{
  int i, j, k, m = rows(A), n = cols(A);
 
  double e, tol = max(m,n)*eps*max(sum(abs(A))); /* max abs row sum */

  matrix R = A;  vector t;

  for( j = i = 1; j <= n && i <= m; ++j) /* for each column */
  {
    t = abs(R(i:m,j));

    e = max(t(:));

    if( e < tol) /* set remainder of this column to all zeros */
    {
      R(i:m,j) = zeros( m-i+1, 1);
      continue;
    }

    k = find(t==e) + i - 1; /* index of max element */

    R([i k],j:n) = R([k i],j:n); /* swap rows i and k */

    R(i,j:n) = R(i,j:n)/R(i,j); /* normalize row i */

    for( k = 1; k <= m; ++k) /* put zeros above and below (i,j) element) */
    {
      if( k == i) continue;

      R(k,j:n) = R(k,j:n) - R(k,j)*R(i,j:n);
    }

    ++i; /* next row */
  }

  return R;
}
