int qrp( matrix A, Q, R, vector p)
/*
 * QR with column interchanges, Q * R = A(:,p), return value is rank(A)
 */
{
    int i, j, k, q, m = rows(A), n = cols(A);
    double alpha, beta, s, t, tol = 10 * eps * max(+/abs(A')); /* 1-norm */
    vector v;

    if( narg != 4) {
	print("qrp( A, Q, R, p) takes four arguments");
	return 0;
    }

    R = A;
    Q = eye( m, m);
    p = [1:n];
    q = min( [m, n]);

    for( j = 1; j <= q; ++j) {		/* go across the columns */

	k = j;				/* find largest remaining column */
	t = +/abs(R(j:m,j)');		/*  using 1-norm */
	for( i = j+1; i <= n; ++i)
	    if( (s = +/abs(R(j:m,i)')) > t) {
		k = i;
		t = s;
	    }

	if( t < tol)			/* remaining matrix is negligable */
	    return j-1;

	if( k != j) {			/* swap indices and cols */
	    p([j k]) = p([k j]);
	    R(:,[j k]) = R(:,[k j]);
	}

	if( j == m)			/* done, full rank */
	    break;

	v = R(j:m,j)/t;			/* setup Householder transformation */
	alpha = norm(v);
	beta = alpha * ( alpha + abs(v(1)));
	v(1) = v(1) + (v(1) > 0 ? alpha : -alpha);

	R(j,j) = v(1) > 0 ? -t * alpha : t * alpha;	/* transform R */
	R(j+1:m,j) = 0;
	for( i = j+1; i <= n; ++i)
	    R(j:m,i) = R(j:m,i) - v * ((v'*R(j:m,i))/beta);

	for( i = 1; i <= m; ++i)			/* transform Q */
	    Q(i,j:m) = Q(i,j:m) - v' * ((Q(i,j:m)*v)/beta);
    }
    return q;	/* full rank */
}
