QR-Zerlegung mit dem Gram-Schmidt-Verfahren

program GramS(input, output);

    {Gram-Schmidtsches Orthogonalisierungsverfahren}
    
    uses
        stdfile;
        
    const
        MaxArray = 32;
        
    var
        a, r:    array [1..MaxArray, 1..MaxArray] of double;
        n, m:    integer;
        i, j, k: integer;
        sum:     double;
                
    procedure init;
        var
            fName: string;
            f:     text;
            t:     double;
            i:     integer;
        begin {init}
            writeln('Gram-Schmidtsches Orthogonalisierungsverfahren');
            writeln;
            writeln('Welche Datei enthaelt die Matrix? (In grams.dat liegt ein Beispiel.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                if not StdRead(f, t) then begin
                    writeln('Kann erste Dimension nicht lesen.');
                    halt
                end; {if}
                if t > MaxArray then begin
                    writeln('Erste Dimension zu gross (groesser als ', MaxArray, '.');
                    halt
                end; {if}
                n := round(t);
                if not StdRead(f, t) then begin
                    writeln('Kann zweite Dimension nicht lesen.');
                    halt
                end; {if}
                if t > MaxArray then begin
                    writeln('Zweite Dimension zu gross (groesser als ', MaxArray, '.');
                    halt
                end; {if}
                m := round(t);
                for i := 1 to n do begin
                    if not StdReadn(f, m, @a[i, 1]) then begin
                        writeln('Fehler beim Lesen der Matrix');
                        halt
                    end {if}
                end; {for}
                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('Man erhaelt die Matrizen');
            writeln;
            for i := 1 to n do begin
               for j := 1 to m do begin
                   write(a[i, j]:6:3)
               end; {for}
               writeln
            end; {for}
            writeln;
            for i := 1 to n do begin
               for j := 1 to m do begin
                   write(r[i, j]:6:3)
               end; {for}
               writeln
            end; {for}
        end; {result}

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