rb_residual_norm.m

function [n,nn,nnn] = rb_residual_norm(model, rd, rbsim)
% RB_RESIDUAL_NORM Calculate the residual norm of the rb solution rbsim
%
%  Andreas Schmidt, 2015
%
descr = model.descr;

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();

% Assemble all matrices:
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);
invR = inv(R);

R_ECQCA = zeros(size(PN,1));
for i = 1:length(A_coeff)
    for j = 1:length(C_coeff)
        for k = 1:length(C_coeff)
            for l = 1:length(E_coeff)
                R_ECQCA = R_ECQCA + A_coeff(i)*C_coeff(j)*C_coeff(k)*...
                    E_coeff(l)*estim.M1b{l,k}*Q*estim.M1a{j,i};
            end
        end
    end
end

R_EE = zeros(size(PN,1));
for i = 1:length(E_coeff)
    for j = 1:length(E_coeff)
        R_EE = R_EE + E_coeff(i)*E_coeff(j)*estim.M8{i,j};
    end
end

R_AA = zeros(size(PN,1));
for i = 1:length(A_coeff)
    for j = 1:length(A_coeff)
        R_AA = R_AA + A_coeff(i)*A_coeff(j)*estim.M5{i,j};
    end
end

R_BB = zeros(size(PN,1));
for i = 1:length(B_coeff)
    for j = 1:length(B_coeff)
        R_BB = R_BB + B_coeff(i)*B_coeff(j)*rd.B_comp{i}*invR*rd.B_comp{j}';
    end
end

R_EA = zeros(size(PN,1));
for i = 1:length(A_coeff)
    for j = 1:length(E_coeff)
        R_EA = R_EA + A_coeff(i)*E_coeff(j)*estim.M3{i,j};
    end
end

R_ECCE = zeros(size(PN,1));
for i = 1:length(C_coeff)
    for j = 1:length(C_coeff)
        for k = 1:length(E_coeff)
            for l = 1:length(E_coeff)
                R_ECCE = R_ECCE + C_coeff(i)*C_coeff(j)*E_coeff(k)*...
                    E_coeff(l)*estim.M1b{l,j}*Q*estim.M1b{k,i}';
            end
        end
    end
end

R_CCTQCCTQ = zeros(size(estim.M7{1,1}, 1));
for i = 1:length(C_coeff)
    for j = 1:length(C_coeff)
        for k = 1:length(C_coeff)
            for l = 1:length(C_coeff)
                 R_CCTQCCTQ = R_CCTQCCTQ + estim.M7{i,j}*Q*estim.M7{k,l}*Q*C_coeff(i)*C_coeff(j)*...
                    C_coeff(k)*C_coeff(l);
            end
        end
    end
end 

d = trace(R_CCTQCCTQ);
res = d + 4*trace(R_ECQCA*PN) - 4*trace(PN*R_EE*PN*R_BB*PN*R_EA) + ...
    -2*trace(R_ECCE*PN*R_BB*PN) + trace(R_EE*PN*R_BB*PN*R_EE*PN*R_BB*PN) ...
    + 2*trace(R_EE*PN*R_AA*PN) + 2*trace(R_EA'*PN*R_EA'*PN);

n = sqrt(abs(res));
nn = n/sqrt(d);

% Set up the norm of A and the norm of B:
if nargout > 2
    normA = 0;
    for k = 1:length(A_coeff)
        for l = 1:length(A_coeff)
            normA = normA + A_coeff(k)*A_coeff(l)*estim.normA{k,l};
        end
    end
    normA = sqrt(normA);
    normB = 0;
    for k = 1:length(B_coeff)
        for l = 1:length(B_coeff)
            normB = normB + B_coeff(k)*B_coeff(l)*estim.normB{k,l};
        end
    end
    normB = sqrt(normB);
    


    nPhat = sqrt(trace(estim.WtW*PN*estim.WtW*PN));
    normE = abs(E_coeff(1)*estim.normE);
    nnn = n / (2*nPhat*normE*normA + sqrt(d) + normE^2*nPhat^2*norm(invR)*normB^2);
end
end