Matrixinversion nach Gauß-Jordan

program GJordan(input, output);

    {Berechnung der inversen Matrix nach Gauss-Jordan}
    
    uses
        stdfile;
    
    const
        MaxArray = 32;
        
    var
        n:       integer;
        a, b:    array [1..MaxArray, 1..MaxArray] of double;
        i, j, k: integer;
        t:       double;
                        
    procedure init;
        var
            fName: string;
            f:     text;
            t:     double;
            i:     integer;
        begin {init}
            writeln('Berechnung der inversen Matrix nach Gauss-Jordan');
            writeln;
            writeln('In welcher Datei steht die Matrix? (In gjordan.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}
                writeln;
                close(f)
            end else begin {if}
                writeln('Kann Datei "', fName, '" nicht oeffnen.');
                halt
            end {if}
        end; {init}

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

    begin {GJordan}
        init;
        for k := 1 to n do begin
            for j := 1 to n do begin
                b[k, j] := 0
            end; {for}
            b[k, k] := 1
        end; {for}
        for k := 1 to n do begin
            t := a[k, k];
            for j := 1 to n do begin
                a[k, j] := a[k, j] / t;
                b[k, j] := b[k, j] / t
            end; {for}
            for i := 1 to n do begin
                if i <> k then begin
                    t := a[i, k];
                    for j := 1 to n + 1 do begin
                        a[i, j] := a[i, j] - t * a[k, j];
                        b[i, j] := b[i, j] - t * b[k, j]
                    end {for}
                end {if}
            end {for}
        end; {for} 
        result     
    end. {GJordan}