function [X, R, p] = slls(A, B, T)
%
% Sparse Linear Least Squares Problems Solver
%
% X         = slls(A, B)
% [X, R]    = slls(A, B)
% [X, R, p] = slls(A, B)
% R         = slls(A)
% [R, p]    = slls(A)
% X         = slls(A, B, R)
%
% A  an m-by-n real matrix 
% B  an m-by-r real matrix 
% R  an m-by-n real upper triangular matrix 
% X  an n-by-r real matrix 
% p  a 1-by-n real matrix 
%
% X = slls(A, B) returns X = A\B, the least squares solution to A*X = B.
% The rank of A is determined without numeric column pivoting. If the 
% rank of A is k, then the solution X has at most k nonzeros per column.
% In solving A*X = B, the sparse QR factorization A(:,p) = Q*R is computed,
% where Q is an m-by-m orthogonal matrix, R is an m-by-n upper triangular 
% matrix and p contains a column minimum degree ordering of A used for 
% preserving the sparsity of R.  The matrix Q is not explicitly computed.
% Matrices R and p are not returned.
%
% [X, R] = slls(A, B) is the same as the function X = slls(A, B) except 
% that (1) sparse QR factorization A = Q*R is computed, i.e., no column 
% permutation on A is performed, and (2) the upper triangular matrix R, 
% which satisfies A = Q*R, is returned.
%
% [X, R, p] = slls(A, B) is exactly the same as the function X = slls(A, B) 
% except that matrices R and p are returned.
%
% R = slls(A) returns the sparse upper triangular matrix R that satisfies 
% A = Q*R, where Q is an m-by-m orthogonal matrix which is not explicitly 
% computed.
%
% [R, p] = slls(A) returns the sparse upper triangular matrix R that 
% satisfies A(:,p) = Q*R, where p contains a column minimum degree ordering 
% of A used for preserving the sparsity of R and Q is an m-by-m orthogonal 
% matrix which is not explicitly computed.
%
% X = slls(A, B, R) returns the least squares solution to A*X = B.
% The matrix A has been factored as Q*R, where Q is an m-by-m orthogonal 
% matrix and R is an m-by-n upper triangular matrix. If Q is not available 
% and R is available, this function obtains the least squares solution to 
% A*X = B by solving the semi-normal equations R'*R*X = A'*B, and then 
% performing two steps of iterative refinement. The matrix A must have 
% full column rank.
%
% Copyright (C) 1995, Chunguang Sun
% Advanced Computing Research Institute
% Cornell Theory Center, Cornell University.
%

if (nargin < 1) 
    error('SLLS requires one or two input arguments.'); 
end

[m, n] = size(A);
if (m < 1 | n < 1), error('Input argument one is empty.'); end

if (nargin == 1)
    % R = slls(A) or
    % [R, p] = slls(A)

    if (nargout <= 1)      
	% R = slls(A)

        S = spones(A);
        H = S'*S+speye(n);

        if (issparse(A))
            [f, U] = mexslls(H, A);
        else
            [f, U] = mexslls(H, sparse(A));
        end

        X = U(1:m,1:n);
	if (~issparse(A)) X = full(X); end
    else 
	% [R, p] = slls(A)

        p = colmmd(A);
        AP = A(:,p);
        S = spones(AP);
    	H = S'*S+speye(n);

        if (issparse(AP))
            [f, U] = mexslls(H, AP);
        else
            [f, U] = mexslls(H, sparse(AP));
        end

        X = U(1:m,1:n);
	if (~issparse(A)) X = full(X); end

        R = p;
    end
elseif (nargin == 2)
    [mb,h] = size(B);

    if (mb ~= m)
    	error('Input argument one and input argument two do not match.');
    end

    % X = slls(A, B)
    % [X, R] = slls(A, B)
    % [X, R, p] = slls(A, B)

    if (nargout == 2)      
	% [X, R] = slls(A, B)

        AB = [A B];
        S = spones(AB);
        H = S'*S+speye(n+h);

        if (issparse(AB))
            [f, U, cols] = mexslls(H, AB);
        else
            [f, U, cols] = mexslls(H, sparse(AB));
        end

    	rank = sum(cols<=n);

    	if (rank == n)
            X = U(1:n,1:n)\U(1:n,(n+1):(n+h));
            if (~issparse(A) | ~issparse(B)), X = full(X); end
    	else
            if (rank < min(m,n))
            	disp(['Warning: Rank deficient, rank = ' int2str(rank)]);
            end

            rcols = cols(1:rank);
            Y = U(1:rank,rcols)\U(1:rank,(n+1):(n+h));

            if (issparse(A) & issparse(B))
            	X = sparse([], [], [], n, h, nnz(Y));
            else
            	X = zeros(n,h);
            end

            X(rcols,:) = Y;
    	end

        R = U(1:m,1:n);
	if (~issparse(A)) R = full(R); end
    else
	% X = slls(A, B)
	% [X, R, p] = slls(A, B)

    	p = colmmd(A);
    	APB = [A(:,p) B];
    	S = spones(APB);
    	H = S'*S+speye(n+h);

    	if (issparse(APB))
    	    [f, U, cols] = mexslls(H, APB);
    	else
    	    [f, U, cols] = mexslls(H, sparse(APB));
    	end

        rank = sum(cols<=n);

        if (rank == n)
            X = U(1:n,1:n)\U(1:n,(n+1):(n+h));
            if (~issparse(A) | ~issparse(B)), X = full(X); end
            X(p,:) = X;
        else
            if (rank < min(m,n))
                disp(['Warning: Rank deficient, rank = ' int2str(rank)]);
            end

            rcols = cols(1:rank);
            Y = U(1:rank,rcols)\U(1:rank,(n+1):(n+h));

            if (issparse(A) & issparse(B))
                X = sparse([], [], [], n, h, nnz(Y));
            else
                X = zeros(n,h);
            end

            X(rcols,:) = Y;
            X(p,:) = X;
        end

    	if (nargout >= 3)
	    R = U(1:m,1:n);
            if (~issparse(A)) R = full(R); end
    	end
    end
else
    % X = slls(A, B, T)

    [mb,h] = size(B);

    if (mb ~= m)
    	error('Input argument one and input argument two do not match.');
    end

    [mt,nt] = size(T);
    if (mt < n | nt ~= n) 
        error('The third input argument is invalid.');
    elseif (mt > n) 
        T = T(1:n,1:n);
    end

    TT = T';
    AT = A';
    
    X = T\(TT\(AT*B));           % solve the semi-normal equations

    X = X + T\(TT\(AT*(B-A*X))); % the first step of iterative refinement
    X = X + T\(TT\(AT*(B-A*X))); % the second step of iterative refinement
end

f1 = flops + f;
flops(f1);
