function [ssI,sI,dI] = inv_band5(ssB,sB,dB)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the elements in the sub-sub diagonal, sub-diagonal and
% diagonal of the inverse of a symmetric positive definite matrix
% B which has bandwidth 5.
%
% Note that the off-band elements in the inverse B are not necessarily
% zero, but we do not compute them here because they are not needed for
% our spline application.
%
% This uses the Hutchinson and de Hoog algorithm described in section 3.2.2 
% Green and Silverman (1994).
%
% INPUTS:
%      ssB    Nx1 vector of sub-subdiagonal elements of B (uses 1st N-2)
%      sB     Nx1 vector of subdiagonal elements of B (uses 1st N-1)
%      dB     Nx1 vector of diagonal elements of B
%
% OUTPUT:
%      ssI    Nx1 vector containing the sub-sub-diagonal elements of inv(B) (uses 1st N-2)
%      sI     Nx1 vector containing the sub-diagonal elements of inv(B) (uses 1st N-1)
%      dI     Nx1 vector containing the diagonals of inv(B)
%
% Written by Karen Chiswell on March 23, 2006
% Modified:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = length(dB);

%       initialize output vectors
dI = dB;
sI = sB;
ssI = ssB;

%       get cholesky factorization of B using chol_band5.m
[ssL,sL,d] = chol_band5(ssB,sB,dB);

%       Compute the central diagonals of inv(B):  see p. 34 of Green & Silverman

%       start in bottom right corner of inv(B)
%       note that although the formulas are for the upper triangular
%          half of inv(B), by symmetry they are also for the lower 
%          triangular part
dI(N) = 1/d(N);              % [invB]_NN
sI(N-1) = -sL(N-1)*dI(N);    % [invB]_N,N-1
dI(N-1) = 1/d(N-1) - sL(N-1)*sI(N-1);      % [invB]_N-1,N-1

%       now for the remaining elements
for i = N-2:-1:1
    ssI(i) = -sL(i)*sI(i+1) - ssL(i)*dI(i+2);       % [invB]_i+2,i
    sI(i)  = -sL(i)*dI(i+1) - ssL(i)*sI(i+1);       % [invB]_i+1,i
    dI(i)  = 1/d(i) - sL(i)*sI(i) - ssL(i)*ssI(i);  % [invB]_i,i
end   % loop on i

% end of inv_band5.m