LR-Zerlegung nach Crout

program Crout(input, output);

    uses
        stdfile;
    
    const
        MaxArray = 32;
        
    var
        n:       integer;
        a:       array [1..MaxArray, 1..MaxArray] of double;
        i, j, k: integer;
        sum:      double;
                        
    procedure init;
        var
            fName: string;
            f:     text;
            t:     double;
            i:     integer;
        begin {init}
            writeln('LR-Dekomposition nach Crout');
            writeln;
            writeln('In welcher Datei steht die Matrix? (In crout.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(a[i, j]:7:3, ' ')
                end; {for}
                writeln
            end {for}
        end; {result}

    begin {Crout}
        init;
        for j := 2 to n do begin
            a[1, j] := a[1, j] / a[1, 1];
        end; {for}
        for j := 2 to n - 1 do begin
            for i := j to n do begin
                sum := 0;
                for k := 1 to j - 1 do begin
                   sum := sum + a[i, k] * a[k, j]
                end; {for}
                a[i, j] := a[i, j] - sum
            end; {for}
            for k := j + 1 to n do begin
                sum := 0;
                for i := 1 to j - 1 do begin
                    sum := sum + a[j, i] * a[i, k]
                end; {for}
                a[j, k] := (a[j, k] - sum) / a[j, i]
            end; {for}
        end; {for}
        sum := 0;
        for k := 1 to n - 1 do begin
            sum := sum + a[n, k] * a[k, n]
        end; {for}
        a[n, n] := a[n, n] - sum;
        result     
    end. {Crout}