function [Q,R,nflops,rk] = compr(A,typ,par)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Compression routine (rank-revealing QR)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opt = 1; % compression routine option
% opt = 1: use a matlab rank-revealing QR routine; may be slow in large matlab test
% opt = 2: use mex routines which call Fortran rank-revealing QR;
%          the Fortran routines should be compiled in Matlab in advance:
%              mex mexmgsclpv.f90 mgsclpv.f90
%              mex mexzmgsclpv.f90 zmgsclpv.f90
% opt = 3: use the routine larank in the paper by T.-L. Lee, et al.:
%          A rank-revealing method with updating, downdating, and applications.
%          Part II, SIAM. J. Matrix Anal. Appl., 31 (2009), pp. 503--525.
%          (Works for Windows, not tested for Unix)
% opt = 4: use Matlab QR and then truncate; for test purpose only!!!
%          If matlab SVD is used, the result is more accurate.

switch opt
    case 1
        [Q,R,nflops,rk] = mgsclpv(A,typ,par);
        rk = min([rk,size(R)]);
        Q = Q(:,1:rk); R = R(1:rk,:);
    case 2
        if length(typ) == 3 && min(typ == 'tol')
            t = int32(0); par = double(par);
        else
            t = int32(1); par = int32(par);
        end
        if isreal(A)
            [Q,R,rk,nflops] = mexmgsclpv(A,t,par);
        else
            [Q,R,rk,nflops] = mexzmgsclpv(A,t,par);
        end
        Q = Q(:,1:rk); R = R(1:rk,:);
        %nflops = double(4*numel(A)*rk-2*rk*rk*(size(A,1)+size(A,2))+4*rk^3/3);
        % this nflops is for unix only, temp
    case 3
        if length(typ) == 3 && min(typ == 'tol')
            if size(A,1)*10 < size(A,2) % to handle a bug in larank
                [rk,U,R,Q] = larank(A',par); R = R'*U';
            else
                [rk,Q,R,V] = larank(A,par); R = R*V';
            end
        else
            [Q,R,V] = svd(A); R = R*V';
            rk = min([par,size(R)]);
            Q = Q(:,1:rk); R = R(1:rk,:);
            %error('This compression option only works with a tolerance')
        end
        % this routine does not seem to provide an explicit count
        nflops = double(4*numel(A)*rk-2*rk*rk*(size(A,1)+size(A,2))+4*rk^3/3);
    case 4
        %%% for test purpose!!!
        if min(size(A)) == 0
            Q = zeros(size(A,1),0);
            R = zeros(0,size(A,2));
            rk = 0; nflops = 0;
            return
        end

        [Q,R,p] = qr(full(A),0);
        for i = 1:length(p); rp(p(i)) = i; end

        if length(typ) == 3 && min(lower(typ) == 'tol')
            if min(size(R)) == 1
                rk = 1;
            else
                rk = find(abs(diag(R)) >= par,1,'last');
                % diag may creat problem for a vector
            end
        else
            rk = par;
        end

        if isempty(rk); rk = 1; end
        nflops = double(4*numel(A)*rk-2*rk*rk*(size(A,1)+size(A,2))+4*rk^3/3);
        rk = min([rk,size(R)]);
        Q = Q(:,1:rk); R = R(1:rk,rp);
end
% fprintf('compression flops: %8.2e\n',nflops);
% maxerr(Q*R,A)