function x=mnhf_thomas(A,b)
%MNHF_THOMAS implements the Thomas algorithm.
%
% A - n-by-n tridiagonal matrix 
%     (do not verify tridiagonal structure)
% b - n-by-1 RHS column vector
%
% Usage  mnhf_thomas(A,b)

[r,c] = size(A);   % Compute the size of A.
[rr,cc] = size(b); % Compute the size of b.
if cc==1
 if isequal(r,c)
  if isequal(c,rr)
   
   % Extract elements along main/upper/lower diagonals.
   d = diag(A,0);
   a = diag(A,-1);
   c = diag(A,1);
   
   % Eliminate elements below the main diagonal.
   for j = 2:r
    d(j) = d(j)-a(j-1)/d(j-1)*c(j-1);
    b(j) = b(j)-a(j-1)/d(j-1)*b(j-1);
   end
   
   % Perform back substitution.
   b(r) = b(r)/d(r);
   for j = r-1:-1:1
    b(j) = (b(j)-c(j)*b(j+1))/d(j);
   end
   
   % Return solution.
   x = b;
   
  else
   fprintf('Size mismatch between matrix and vector.') 
  end
 else
  fprint('Matrix is not square.')
 end
else
 fprint('b is not a column vector.');
end