INVERSE_TRIDIAG.m

function [T]=INVERSE_TRIDIAG(A)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%          Berechnung der Inversen der Tridiagonalmatrix A            %%%
%%% mit HMGTI-Algorithmus -- nach dem paper von El-Mikkawy & Karawia    %%%
%%% (Inversion of general tridiagonal matrices, 2006)                   %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function [T]=INVERSE_TRIDIAG(A)
%
% Input:   A     Tridiagonalmatrix in IR^{nxn}
% Output:  T     Matrix in IR^{nxn} - Inverse der Tridiagonalmatrix A, s.d.
%                                      AT=Id


n = size(A,1);    % Dimension
d = diag(A);      % Diagonalen auslesen
a = diag(A,1);
b = [0 ;diag(A,-1)];

alpha = zeros(1,n+1);
beta = zeros(1,n+1);

alpha(1) = 1;     % Initialisieren
alpha(2) = d(1);
beta(n+1) = 1;
beta(n) = d(n);

T = zeros(n,n);

% Berechnung alpha
for i=2:n
    alpha(i+1) = d(i)*alpha(i)-b(i)*a(i-1)*alpha(i-1);
end

if alpha(n+1)==0
    error('no inverse exists')
end

% Berechnung beta
beta(1) = alpha(n+1);

for i=n-1:-1:2
    beta(i) = d(i)*beta(i+1)-b(i+1)*a(i)*beta(i+2);
end

if (min(abs(alpha))==0 | min(abs(beta))==0)
    error('Failure')
end

% Bestimmung der Werte auf der Diagonalen
t11 = d(1)-(b(2)*a(1)*beta(3))/beta(2);
T(1,1) = 1/t11;

tnn = d(n)-(b(n)*a(n-1)*alpha(n-1))/alpha(n);
T(n,n) = 1/tnn;

for i=2:n-1
    tii = d(i)-(b(i)*a(i-1)*alpha(i-1))/alpha(i)-(b(i+1)*a(i)*beta(i+2))/beta(i+1);
    T(i,i) = 1/tii;
end


% Berechnung der restlichen Eintraege
% fuer j>i (obere Dreiecksmatrix)
for j=2:n
    for i=1:j-1
        c1 = 1;
        for k=1:j-i
            c1=c1*a(j-k);
        end
        T(i,j) = (-1)^(j-i)*c1*alpha(i)/alpha(j)*T(j,j);
    end
end

% fuer i>j (untere Dreiecksmatrix)
for i=2:n
    for j=1:i-1
        c2=1;
        for k=1:i-j
            c2=c2*b(j+k);
        end
        T(i,j) = (-1)^(i-j)*c2*beta(i+1)/beta(j+1)*T(j,j);
    end
end