riccati_gen_reduced_data.m

Go to the documentation of this file.

function reduced_data = riccati_gen_reduced_data(model, detailed_data)
%reduced_data = riccati_gen_reduced(model, detailed_data)
% 
%  Build the reduced data, i.e. all matrices and components that are
%  required for the calculation of the reduced solution and the residual
%  norm in an online efficient manner.
%
% Andreas Schmidt, 2015
W = detailed_data.RB;
if ~isfield(detailed_data, 'RB_V');
    V = W;
else
    V = detailed_data.RB_V;
end
% 1. Project all the matrices:
model_data = gen_model_data(model);

% For copying all the coefficient functions:
reduced_data = model_data;
reduced_data.RB = V;

reduced_data.VtV = V'*V;
Ak = {};
for k = 1:length(model_data.A_components)
    Ak{k} = W'*model_data.A_components{k}*V;
end
reduced_data.A_components = Ak;
reduced_data.B_components = {};
for k = 1:length(model_data.B_components)
    reduced_data.B_components{k} = W'*model_data.B_components{k};
end
reduced_data.C_comp = {};
for k = 1:length(model_data.C_components)
    reduced_data.C_components{k} = model_data.C_components{k}*V;
end
reduced_data.E_comp = {};
for k = 1:length(model_data.E_components)
    reduced_data.E_components{k} = W'*model_data.E_components{k}*V;
end
reduced_data.x0_components = {};
for k = 1:length(model_data.x0_components)
    reduced_data.x0_components{k} = W'*model_data.x0_components{k};
end

% For the error estimator:
estim.M1 = {};
for i=1:length(model_data.B_components)
    for j=1:length(model_data.B_components)
        estim.M2{i,j} = (W'*model_data.B_components{i})*(model_data.B_components{j}'*W);
    end
end
estim.M3 = {};
for i=1:length(model_data.C_components)
    for j=1:length(model_data.C_components)
        for k=1:length(model_data.E_components)
            for l=1:length(model_data.E_components)
                estim.M4{i,j,k,l} = (W'*model_data.E_components{k}*model_data.C_components{i}')*(model_data.C_components{j}*model_data.E_components{l}'*W);
            end
        end
    end
end

estim.M5 = {};
estim.M6 = {};
for i = 1:length(model_data.A_components)
    for j = 1:length(model_data.C_components)
        for k = 1:length(model_data.C_components)
            for l = 1:length(model_data.E_components)
                estim.M1{i,j,k,l} = W'*model_data.E_components{l}*model_data.C_components{j}'*model_data.C_components{k}*(model_data.A_components{i})'*W;
            end
        end
    end
    for j = 1:length(model_data.E_components)
        estim.M3{i,j} = W'*model_data.E_components{j}*(model_data.A_components{i})'*W;
        estim.M6{i,j} = estim.M3{i,j}';
    end
    for j = 1:length(model_data.A_components)
        estim.M5{i,j} = W'*model_data.A_components{i}*(model_data.A_components{j}')*W;
    end
   
end
for i = 1:length(model_data.C_components)
    for j = 1:length(model_data.C_components)
        for k = 1:length(model_data.C_components)
            for l = 1:length(model_data.C_components)
                estim.M7{i,j,k,l} = trace(model_data.C_components{i}*model_data.C_components{j}'*model_data.C_components{k}*model_data.C_components{l}');
            end
        end
    end
end
for i = 1:length(model_data.E_components)
    for j = 1:length(model_data.E_components)
estim.M8{i,j} = W'*model_data.E_components{i}*(model_data.E_components{j}')*W;
    end
end


% Finally the last components for the error estimator:
estim.normE = normest(model_data.E_components{1});


% Prepare the value for ||BB^T||
estim.normBBT = {};
estim.normBBT_use2norm = false;
if length(model_data.B_components) == 1
    estim.normBBT = normest(model_data.B_components{1});
    estim.normBBT_use2norm = true;
else
    for i = 1:length(model_data.B_components)
        for j = 1:length(model_data.B_components)
            for k = 1:length(model_data.B_components)
                for l = 1:length(model_data.B_components)
                    estim.normBBT{i,j,k,l} = trace(model_data.B_components{k}'*model_data.B_components{l});
                end
            end
        end
    end
end

% Gamma:
if isfield(detailed_data,'gamma')
    estim.gamma = detailed_data.gamma;

    if strcmp(estim.gamma.mode, 'kernel_interpolation')
        % Build the kernel weights
        mus = estim.gamma.interpol_mus;
        gammas = estim.gamma.interpol_gammas;
        k = riccati_kernel_functions(model, estim.gamma.kernel);
        % Calculate the alphas:
        K = zeros(length(gammas));
        for i = 1:length(gammas)
            for j = 1:length(gammas)
                K(i,j) = k(mus(i,:),mus(j,:));
            end
        end

        estim.gamma.alpha = K\gammas;
    elseif strcmp(estim.gamma.mode, 'exponential_bound')
        if length(model_data.E_components) > 1
            error('Only one component in E is currently supported!');
        end
        
        % NOTE: This is a workaround for a real strange and nasty "bug" in
        % Matlabs svds function:
        converged = false;
        while ~converged
            estim.normEInv = 1/svds(model_data.E_components{1},1,0);
            
            if ~isempty(estim.normEInv)
                converged = true;
            end
        end
    end
end

reduced_data.estim = estim;