Thomas-Algorithmus

program Thomas(input, output);

    {Loesung von Ax = c mit tridiagonaler Matrix}
    
    uses
        stdfile;
    
    const
        MaxArray = 32;
        
    var
        n:             integer;
        a:             array [1..MaxArray, 1..MaxArray] of double;
        c, e, f, g, x: array [1..MaxArray] of double;
        i, j, k:       integer;
                        
    procedure init;
        var
            fName: string;
            f:     text;
            t:     double;
            i:     integer;
        begin {init}
            writeln('Tridiagonaler Thomasalgorithmus');
            writeln;
            writeln('In welcher Datei steht die Matrix mit Loesungsvektor?');
            writeln('(In thomas.dat liegt ein Beispiel.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                if not StdRead(f, t) then begin
                    writeln('Kann Dimension nicht lesen.');
                    halt
                end; {if}
                if t > MaxArray then begin
                    writeln('Dimension zu gross (groesser als ', MaxArray,
                        ').');
                    halt
                end; {if}
                n := round(t);
                for i := 1 to n do begin
                    if not StdReadn(f, n, @a[i, 1]) then begin
                        writeln('Fehler beim Lesen der ', i, '-ten Zeile.');
                        halt
                    end {if}
                end; {for}
                if not StdReadn(f, n, @c[1]) then begin
                    writeln('Kann Loesungsvektor nicht lesen.');
                    halt
                end; {if}
                writeln;
                close(f)
            end else begin {if}
                writeln('Kann Datei "', fName, '" nicht oeffnen.');
                halt
            end {if}
        end; {init}

    procedure result;
        var
            i: integer;
        begin {result}
            writeln('Das Ergebnis lautet');
            writeln;
            for i := 1 to n do begin
                writeln(x[i]:7:3)
            end {for}
        end; {result}

    begin {Thomas}
        init;
        for i := 2 to n do begin
            e[i] := a[i, i - 1]
        end; {for}
        for i := 1 to n do begin
            f[i] := a[i, i]
        end; {for}
        for i := 1 to n - 1 do begin
            g[i] := a[i, i + 1]
        end; {for}
        {Dekomposition}
        for k := 2 to n do begin
            e[k] := e[k] / f[k - 1];
            f[k] := f[k] - e[k] * g[k - 1]
        end; {for}
        {Vorwaertselimination}
        for k := 2 to n do begin
            c[k] := c[k] - e[k] * c[k - 1]
        end; {for}
        {Ruecksubstitution}
        x[n] := c[n] / f[n];
        for k := n - 1 downto 1 do begin
            x[k] := (c[k] - g[k] * x[k + 1]) / f[k]
        end; {for}
        result     
    end. {Thomas}