function [ eta , V_ygivenx ] = blip0(V_xx,V_xy,V_yy)
% BLIP0(V_xx,V_yx,V_yy) sets up basic arrays for computing BLiP and
% prediction covariance.
% INPUTS:
% V_xx = n by n covariance matrix of given x
% V_xy = m by n covariance matrix of x with y
% V_yy = m by m covariance matrix of y to be predicted
% OUTPUTS:
% eta: m by n matrix; compute prediction with yhat = eta*x .
% V_ygivenx = m by m matrix of covariance of y-yhat.
% NOTE: Uses cholesky of V_xx; if V_xx is numerically
%   singular, then cholesky stops early.  This is probably
%   OK for must time series applications.
[n m] = size(V_xx);
[T,p] = chol(V_xx);
if p>0
   q = p-1;
else
   q = n;
end
zeta = backsolve(T,V_xy(1:q,:),1,1);
eta = backsolve(T,zeta,1,0);
if q<n
   eta = [ eta ; zeros(n-q,1) ];
end
V_ygivenx = V_yy - zeta'*zeta;