function [R, s, X, f] = sqr(A, t, B)
%
% Sparse QR Factorization
%
% [R, p      ] = sqr(A)
% [R, p, f   ] = sqr(A)
% [R         ] = sqr(A, p)
% [R, f      ] = sqr(A, p)
% [R, p, X   ] = sqr(A, B)
% [R, p, X, f] = sqr(A, B)
% [R, X      ] = sqr(A, p, B)
% [R, X, f   ] = sqr(A, p, B)
%
% A  an m-by-n(m >= n) sparse matrix 
% B  a dense m-by-r matrix 
% R  an m-by-n sparse upper triangular matrix 
% p  a dense row vector of length n on left hand side, or 
%    a dense row/column vector of length n on right hand side.
% X  a dense n-by-r matrix 
% f  a 1-by-1 matrix containing the number of flops executed in sqr.
%
% Let A be an m-by-n(m >= n) sparse matrix. This function 
% computes the sparse QR factorization AP = QR, where P is 
% a column permutation matrix and is represented by a 
% permutation vector p.  If p is not provided as an input
% argument, it is obtained by a multiple minimum degree 
% ordering algorithm. The sparse upper triangular matrix R
% is always returned.
%
% If a dense m-by-r matrix B = [b1, b2, ..., br] is given 
% and the rank of A is n, then the solutions to the sparse 
% linear least squares problems 
%	      ||Axj-bj\||, j = 1, 2, ..., r
% are computed. The dense n-by-r matrix X = [x1, x2, .., xr]
% is returned.
%
% If a dense m-by-r matrix B is given and the rank of A is 
% less than n, then the corresponding dense n-by-r solution 
% matrix X is set to zero.
%

if (~issparse(A)) 
    error('The first input argument must be a sparse matrix.');
end

[m, n] = size(A);

if (m < 1 | n < 1) 
    error('The first input argument must be non-empty.');
end

if (m < n) 
    error('The first input argument does not have enough rows.');
end

solve_ls = 0;

if (nargin == 1)

    if ((nargout == 2 | nargout == 3))
	SYM = A~=0.0;
	H = SYM'*SYM;
   	[flag, num_flops, R, s] = ata_sqr(H, A);
	if (nargout == 3) X = num_flops; end
    else 
        error('Wrong number of input or output arguments.');
    end

elseif (nargin == 2)

    if ((nargout == 1) | (nargout == 2))
	if (issparse(t))
    	    error('The second input argument must be dense.');
	end

	[tm, tn] = size(t);

	if (~((tm == n & tn == 1) | (tm == 1 & tn == n)))
    	    error('The second input argument is invalid.');
	end

	if (tm == n) 
	    if (sort(t) ~= (1:n)')
    	    	error('The second input argument is invalid.');
	    end
	else 
    	    if (sort(t) ~= (1:n))
                error('The second input argument is invalid.');
            end
	end

	SYM = A~=0.0;
	H = SYM'*SYM;
	[flag, num_flops, R] = ata_sqr(H, A, t);
	if (nargout == 2) s = num_flops; end
    elseif ((nargout == 3) | (nargout == 4))
	if (issparse(t))
    	    error('The second input argument must be dense.');
	end

        [tm, tn] = size(t);

        if ((tm ~= m) | (tn < 1))
            error('The second input argument is invalid.');
        end

	solve_ls = 1;

	SYM = A~=0.0;
	H = SYM'*SYM;
	[flag, num_flops, R, s, X] = ata_sqr(H, A, t);
	if (nargout == 4) f = num_flops; end
    else
        error('Wrong number of input or output arguments.');
    end
   
elseif (nargin == 3) 

    if ((nargout == 2) | (nargout == 3))
        if (issparse(t))
            error('The second input argument must be dense.');
        end

        [tm, tn] = size(t);

        if (~((tm == n & tn == 1) | (tm == 1 & tn == n)))
            error('The second input argument is invalid.');
        end

        if (tm == n)
            if (sort(t) ~= (1:n)')
                error('The second input argument is invalid.');
            end
        else
            if (sort(t) ~= (1:n))
                error('The second input argument is invalid.');
            end
        end

        if (issparse(B))
            error('The third input argument must be dense.');
        end

        [bm, bn] = size(B);

        if ((bm ~= m) | ( bn < 1))
            error('The third input argument is invalid.');
        end

	solve_ls = 1;

	SYM = A~=0.0;
	H = SYM'*SYM;
   	[flag, num_flops, R, s] = ata_sqr(H, A, t, B);
	if (nargout == 3) X = num_flops; end
    else
        error('Wrong number of input or output arguments.');
    end

else

    error('Wrong number of input arguments.');

end

if flag
    R = R(1:m,1:n);

    if solve_ls
	disp(' ')
	disp('SQR: The input matrix is rank deficient.')
	disp('SQR: The solution matrix is set to zero.')
	disp(' ')
    end
end
