%
%   reduction compares the rank of AT with the given rankAT
%
%   The vector ips describes the linear dependent columns of [AT, BT].
%   If you delete this columns 
%
%   W=[AT, BT];
%   W(ips)=[];
%
%   then
%
%   Q*R == W
%
%   ok - 1 if the Matrix AT has the expected rank = rankAT
%
function [Q,R,ips,ok]=reduction(AT,rankAT,BT,threshold)
ok=1;
[Q,RAT,E]=qr(full(AT));
rankRAT=sum(abs(diag(RAT))>threshold);
if rankRAT ~= rankAT
    disp('reduction: matrix hasn''t expected rank - z(x) not acceptable or');
    disp('           DAE doesn''t regular');
    ok=0;
    %
    %   only to assign something
    %
    R=[];
    ips=[];
    return
end;
%
%   ips = Permutation of E^T
%
[ips,jj]=find(E');
R=RAT(:,ips);
%
%   ips = Permutation of E
%
[ips,jj]=find(E);
%
%   ips contains the linear dependent columns of AT
%
ips=ips(rankRAT+1:end);
%
% delete of linear dependent columns of R
%
R(:,ips)=[];
%
%   if BT is not empty
%
if size(BT,1) ~=0
    QATT=Q'*BT;
    QATT_2=QATT(rankRAT+1:end,:);
    [Q_2,R_2,E_2]=qr(full(QATT_2));
    rankQATT_2=sum(abs(diag(R_2))>threshold);
    if rankAT+rankQATT_2 < size(AT,1)
        disp('reduction: ???');
        ok=0;
        return
    end;
    [ips2,jj]=find(E_2');
    Q=block_mult(Q,Q_2);
    Qhelp=[QATT(1:rankRAT,:);R_2(:,ips2)];
    [ips2,jj]=find(E_2);
    diff=size(AT,1)-rankRAT;
    Qhelp(:,ips2(diff+1:end))=[];
    R=[R,Qhelp];
    %
    %   summerizing of the unused rows in the Jacobian
    %
    ips=[ips;ips2(diff+1:end)+size(AT,2)];
end;