/*  -*- MaTX -*-
 *
 *  NAME
 *      orth() - Orthogonal basis
 *
 *  SYNOPSIS
 *      Q = orth(A)
 *         Matrix Q;
 *         Matrix A;
 *
 *      Q = orth(A,tol)
 *         Matrix Q;
 *         Matrix A;
 *         Real tol;
 *
 *  DESCRIPTION
 *      orth(A) returns an orthonormal basis for the range of A.
 *      The columns of Q span the range of A and the number of 
 *      columns of Q is the rank(A).
 *
 *      tol is used to determine the rank of A.
 *
 *  SEE ALSO
 *      null, kernel, and rank
 */

Func Matrix orth(A, tol, ...)
	Matrix A;
	Real tol;
{
	Integer rk;
	Matrix Q,P,R;

	error(nargchk(1, 2, nargs, "orth"));
	
	if (nargs == 1) {
		tol = EPS*frobnorm(A);
	}

	{Q,R,P} = qr_p(A);

	rk = Integer(sum(abs(diag2vec(R)) .> tol));

	if (rk > 0) {
		Q = Q(:,1:rk);

		Q = -Q;
		Q(:,rk) = -Q(:,rk);
	} else {
		Q = [];
	}
	return Q;
}