dom_dec_PCA_fixspace.m

function [B_0,B_G] = dom_dec_PCA_fixspace(X,XFix_0,XFix_G,A,A_G, ...
                                            gamma_dofs,PCA_modes,epsilon)
%function [B_0,B_G] = dom_dec_PCA_fixspace(X,XFix_0,XFix_G,A,A_G, ...
%                                           gamma_dofs,PCA_modes,epsilon)
%
% computing PCA fixspaces according to schema 2 for the RB
% extension
%
% Return values:
%  B_0: extension for `X_{N,2}^0`
%  B_G: extension for `X_{N,2}^G`

% I.Maier 19.07.2011

% compute pca modes
X = PCA_fixspace(X,[XFix_0,XFix_G],A,PCA_modes,[],epsilon);

% orthonormalize XFix_G with respect to L2-Gamma inner product

A_GG = sparse(A_G(gamma_dofs,gamma_dofs));

if isempty(XFix_G)
  XFix_G = [];
else
  for i = 1:(size(XFix_G,2)-1)
    for j = (i+1):size(XFix_G,2)
      if isequal(XFix_G(:,i),XFix_G(:,j))
        XFix_G(:,j) = 0;
      end;
    end;
  end;
  
  opts.shape = 'upper';
  opts.type = 'ict';  
  R_M = ichol(A_GG, opts);
  [~,R,E] = qr(R_M * XFix_G(gamma_dofs,:),0);
  XFix_G = XFix_G(:,E);
  
  ind = find(abs(diag(R))>epsilon);
  if ~isempty(ind)
    XFix_G = XFix_G(:,ind);
    Rind = R(ind,ind);
    XFix_G = XFix_G(:,ind) / Rind;
  end;

end;

% separate RB_0, RB_G extension
if ~isempty(XFix_G)
  Xo = X - XFix_G * (XFix_G' * A_G * X);
else
  Xo = X;
end;
clear('A_G');

K = Xo(gamma_dofs,:)' * A_GG * Xo(gamma_dofs,:);
K = 0.5*(K + K');
clear('A_GG');

[ep,vp] = eigs(K);
evalues = diag(vp);

fi_G = find(abs(evalues)>=epsilon);
fi_0 = find(abs(evalues)<epsilon);

B_0 = Xo * ep(:,fi_0);
B_G = X * ep(:,fi_G);

clear('X');
clear('Xo');

B_0(gamma_dofs,:) = 0;

% compute pod-extension
B_0 = PCA_fixspace(B_0,XFix_0,A,[],[],epsilon);