LR-Zerlegung nach Cholesky

program CholeskyLR(input, output);

        {LR-Zerlegung einer symmetrischen positiv definiten Matrix
        nach Cholesky}

    uses
        stdfile;

    const
        MaxArray = 32; {<= 32767}

    var
        i, j, k, n: integer;
        sum:        double;
        a:          array [0..MaxArray,0..MaxArray] of double;

    procedure init;
        var
            fName: string;
            f:     text;
            t:     double;
        begin {init}
            writeln('Darstellung der LR-Zerlegung einer symmetrischen positiv');
            writeln('definiten Matrix nach Cholesky.');
            writeln;
            writeln('In welcher Datei liegen die Matritzen? (In cholesky.dat liegt ein Beispiel.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                if not StdRead(f, t) then begin
                    writeln('Konnte die Dimension nicht lesen.');
                    halt
                end; {if}
                n := round(t);
                if n > MaxArray then begin
                    writeln('Die Dimension ueberschreitet die vorgegebene Arraygroesse.');
                    halt
                end; {if}
                for i:=0 to n-1 do begin
                    if not StdReadn(f, n, @a[i, 0]) then begin
                        writeln('Konnte ', i+1, '-te Zeile nicht lesen.');
                    end; {if}
                end; {for}
            end else begin {if}
                writeln('Die Datei "', fName, '" kann nicht geoeffnet werden.');
            end; {if}
        end; {init}

    procedure result; begin
        for i:=0 to n-1 do begin
            for j:=0 to i do begin
                write(a[i, j]:8:3, ' ')
            end; {for}
            for j:=i+1 to n-1 do begin
                write(0.0:8:3, ' ')
            end; {for}
            writeln;
        end; {for}
    end; {result}

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