function x = tri_ge(a,b)
% Function to solve b = a*x by Gaussian elimination where
% the matrix a is a packed tridiagonal matrix
% Inputs
% a Packed tridiagonal matrix, N by N unpacked
% b Column vector of length N
% Output
% x Solution of b = a*x; Column vector of length N
%* Check that dimensions of a and b are compatible
[N,M] = size(a);
[NN,MM] = size(b);
if( N ~= NN | MM ~= 1)
error('Problem in tri_GE, inputs are incompatible');
end
%* Unpack diagonals of triangular matrix into vectors
alpha(1:N-1) = a(2:N,1);
beta(1:N) = a(1:N,2);
gamma(1:N-1) = a(1:N-1,3);
%* Perform forward elimination
for i=2:N
coeff = alpha(i-1)/beta(i-1);
beta(i) = beta(i) - coeff*gamma(i-1);
b(i) = b(i) - coeff*b(i-1);
end
%* Perform back substitution
x(N) = b(N)/beta(N);
for i=N-1:-1:1
x(i) = (b(i) - gamma(i)*x(i+1))/beta(i);
end
x = x.'; % Transpose x to a column vector
return;