function [TA, TT, TR, NULL, BASIS] = restrictmatrix(varargin)
% [ TA, T, TR, NULL, BASIS ] = restrictmatrix( T, A, [options] )
% Restricts matrices T to a subspace A.
% Also tests, if the matrices are invariant on A.
%
% Input: 
%   T               square matrix or cell array of square matrices
%   A               rectangular matrix giving the columns of the subspace A, to which T will be restricted
%
% Options:
%   'verbose',val   verbose level
%   'epsilon',val   controls the offset to test whether A is invariant or not
%   'smallsize'     Removes unnecessary columns of A (starting from the last column) prior to computation. Returns rank(A) x rank(A) matrix. 
% 
% Output:
%   TA              dimA x dimA matrices. Restriction of T to A
%   TT              dimT x dimT matrices. T in the basis of A complemented to a basis of RR^dimT
%   TR              dim(T-A) x dim(T-A) matrices. Lower right corner of matrix T
%   NULL            Lower left corner of T. If T is A invariant, then NULL consists only of zeros.
%   BASIS           the complemented basis of A. Thus BASIS\T*BASIS == T
%   
% E.g.: vdisp(restrictmatrix([1 1 0; 0 1 1; 1 0 1],[1 -1 0; 0 1 -1]'))
%
% Written by: tommsch, 2018

% XX returns wrong stuff, if size(A,2)>rank(A,2)
% XX Was ist Tout wenn dimA>rankA>dimT ? 



if(~iscell(varargin{1})); TT{1}=varargin{1}; else; TT= varargin{1}; end; %make T a cell array
BASIS=varargin{2}; %subspace
J=size(TT,2); %number of matrices T
dim=size(TT{1},1);
epsilon=parsem('epsilon',varargin,10e-12);
verbose=parsem('verbose',varargin,1);
smallsize=parsem('smallsize',varargin);

RANK=rank(BASIS);
CORANK=dim-RANK;

if(smallsize) %remove linear dependent vectors from A
    for j=size(BASIS,2):-1:1;
        if(rank(BASIS(:,[1:j-1 j+1:end]))==RANK);  BASIS(:,j)=0; end;
    end;
    BASIS=removezero(BASIS,2); %remove the zero columns
end

%make A full rank
SZE=size(BASIS,2); %number of columns of A
BASIS=[BASIS eye(dim)]; %add a base of \RR^dim 
for j=size(BASIS,2):-1:SZE+1; %size(A,2)-dim+1:size(A,2);
    if(rank(BASIS(:,[1:j-1 j+1:end]))==dim);  %remove linear dependent vectors.
        BASIS(:,j)=0; 
    end;
end
ADD=removezero(BASIS(:,SZE+1:end),2); %remove the zero columns from the unit base
BASIS=[BASIS(:,1:SZE) ADD];

if(rank(BASIS(:,1:SZE))~=RANK || rank(BASIS)~=dim);
    error('restrictmatrix: Wrong computation. Err No: 2349852 \n');
end;

TA=cell(size(TT));
TR=cell(size(TT));
NULL=cell(size(TT));
MAXERR=0;
for j=1:J %go through all matrices in T
    TTval=BASIS\TT{j}*BASIS;
    NULL{j}=TTval(end-CORANK+1:end,1:SZE); %the residuum in the lower left corner. Equals zero if the matrix is invariant
    TR{j}=TTval(end-CORANK+1:end,SZE+1:end); 
    TT{j}=TTval; %#ok<AGROW>
    TA{j}=TTval(1:SZE,1:SZE);
end 

for j=1:J
    try
        val=sum(abs(NULL{j}(:)));
        MAXERR=max(double(val), MAXERR);
        if(val>epsilon*dim);
            vprintf('Not invariant for matrix %i. Error = %g >> 0\n',j,sum(abs(NULL{j}(:))),'cpr','err','imp',[1 verbose]);
        end
    catch
        vprintf('Could not test for invariance of matrices. It is likely that the matrices are not invariant. Be careful and check output NULL.\n','cpr','err','imp',[1 verbose]);
        break;
    end
end
vprintf('Maximal error during computation = %g \n',MAXERR,'imp',[2 verbose]);

    
end

function dummy; end %#ok<DEFNU>  %Generates an error, if the 'end' of a function is missing.