rb_residual_norm.m

function [res,nres] = rb_residual_norm(model, rd, rbsim)
% RB_RESIDUAL_NORM Calculate the residual norm of the rb solution rbsim
%
%  Dominik Wittwar, 2016
%

A_coeff = model.A_coeff();
B_coeff = model.B_coeff();
C_coeff = model.C_coeff();
E_coeff = model.E_coeff();
Q_coeff = model.Q_coeff();
R_coeff = model.R_coeff();

% load reduced solution and estim

PN = rbsim.PN;
estim = rd.estim;

% Directly assemble the small R and Q matrices:
Q = lincomb_sequence(rd.Q_comp, Q_coeff);
R = lincomb_sequence(rd.R_comp, R_coeff);
BN = lincomb_sequence(rd.B_comp,B_coeff);
R_P = R + BN'*PN*BN;

% Assemble matrices 

CCT = 0;
for i=1:length(C_coeff)
    for j = 1:length(C_coeff)
        CCT = CCT + estim.CCT{i,j}*C_coeff(i)*C_coeff(j);
    end
end

WTEATW = 0;
for i = 1:length(E_coeff)
    for j = 1:length(A_coeff)
        WTEATW = WTEATW + estim.WTEATW{i,j}*E_coeff(i)*A_coeff(j);
    end
end

WTEETW = 0;
for i = 1:length(E_coeff)
    for j = 1:length(E_coeff)
        WTEETW = WTEETW + estim.WTEETW{i,j}*E_coeff(i)*E_coeff(j);
    end
end

CATW = 0;
for i = 1:length(C_coeff)
    for j = 1:length(A_coeff)
        CATW = CATW + estim.CATW{i,j}*C_coeff(i)*A_coeff(j);
    end
end

WTAATW = 0;
for i = 1:length(A_coeff)
    for j = 1:length(A_coeff)
        WTAATW = WTAATW + estim.WTAATW{i,j}*A_coeff(i)*A_coeff(j);
    end
end

WTECT = 0;
for i = 1:length(E_coeff)
    for j = 1:length(C_coeff)
        WTECT = WTECT + estim.WTECT{i,j}*E_coeff(i)*C_coeff(j);
    end
end

temp = trace((Q*CCT)^2);
res =  temp + trace((PN*WTEETW)*(PN*WTEETW)) -2*trace(Q*WTECT'*PN*WTECT) + ...
    -2*trace(PN*WTEATW*PN*WTEATW') + 2*trace(BN*(R_P\BN')*PN*WTEATW'*PN*WTEATW*PN) + ...
    +2*trace(Q*CATW*PN*CATW') + trace(PN*WTAATW*PN*WTAATW) + ...
    -2*trace(Q*CATW*PN*BN*(R_P\BN')*PN*CATW') + ...
    -2*trace(PN*WTAATW*PN*BN*(R_P\BN')*PN*WTAATW) + ...
    +trace(BN*(R_P\BN')*PN*WTAATW*PN*BN*(R_P\BN')*PN*WTAATW*PN);

res = sqrt(abs(res));
nres = res/sqrt(temp);

end