function S_red=SelectionMatrixSymmetricReduction(n_size,model_order)

%% Reduction of saved variables for symmetric matrices by selection matrix S_red

% S_red is block diagonal of S_red_order0,S_red_order1,S_red_order2,S_red_order3

% S_red_order0=1;
% S_red_order1=eye(n_size);
% S_red_order2 and S_red_order3 are more complicated 

% Inputs:
% n_size: size of kernel, e.g. M+1
% model order: kernel orders, e.g. [0 1 2 3]

% Outputs:
% S_red: selection matrix with zeros and ones

%% Model order 0

if any(model_order==0)
    S_red_order0=sparse(1);
else
    S_red_order0=[];
end


%% Model order 1

if any(model_order>0)
    S_red_order1=speye(n_size);
else
    S_red_order1=[];
end

%% Model order 2

if any(model_order>1)

    row_order2=zeros(n_size.^2,2);
    no=0;
    for ind1=1:n_size
        for ind2=1:n_size
                no=no+1;
                row_order2(no,:)=[ind1 ind2];
            end
    end
    
    [row_order2_unique,ia,ib]=unique(sort(row_order2,2),'rows');
    S_red_order2=sparse(1:size(row_order2,1),ib,1);
    
    
    
else
    S_red_order2=[];
end

%% Model order 3

if any(model_order>2)
    
    row_order3=zeros(n_size.^3,3);
    no=0;
    for ind1=1:n_size
        for ind2=1:n_size
            for ind3=1:n_size
                no=no+1;
                row_order3(no,:)=[ind1 ind2 ind3];
            end
        end
    end
    
    [row_order3_unique,ia,ib]=unique(sort(row_order3,2),'rows');
    S_red_order3=sparse(1:size(row_order3,1),ib,1);


else
    S_red_order3=[];
end

if any(model_order>3)
	error('Kernel order larger than 3 not supported right now');
end

    S_red=blkdiag(S_red_order0,S_red_order1,S_red_order2,S_red_order3);

end