unit2.pas

来自「1、说明： 本书中所有的常用数值算法子过程按书中的章数分别放在以C开头」· PAS 代码 · 共 87 行

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;
procedure QRBKSB(A:matrx2; N:integer; Q:matrx2; B:array of real;
                                              var X:array of real);
procedure QRDCMP(var A:matrx2; M, N:integer;var Q:matrx2);
implementation
procedure QRBKSB(A:matrx2; N:integer; Q:matrx2; B:array of real;
                                              var X:array of real);
var
    I,J:integer;  SUM:real;
begin
    For I:=1 To N do
    begin
        Sum:=0;
        For J:=1 To N do
            Sum:=Sum + Q[I, J] * B[J];
        X[I]:=Sum;
    end;
    For I:=N DownTo 1 do
    begin
        Sum:=X[I];
        For J:=I + 1 To N do
            Sum:=Sum - A[I, J] * X[J];
        If A[I, I] = 0 Then ShowMessage('A is singular matrix.');
        X[I]:=Sum / A[I, I];
    end;
end;

procedure QRDCMP(var A:matrx2; M, N:integer;var Q:matrx2);
var
    I,J,K:integer; S,T,Sgn,H,F:real;
begin
    For I:=1 To M do
    begin
        For J:=1 To M do
            Q[I, J]:=0;
        Q[I, I]:=1;
    end;
    For K:=1 To M - 1  do
    begin
        S:=0;
        For I:=K To M do
            S:=S + Abs(A[I, K]);
        If S <> 0 Then
        begin
            T:=0;
            For I:=K To M do
            begin
                A[I, K]:=A[I, K] / S;
                T:=T + A[I, K] * A[I, K];
            end;
            if A[K,K] >= 0 THEN
                Sgn:=1
            else
                Sgn:=-1;
            T:=-Sqrt(T) * Sgn;
            A[K, K]:=A[K, K] - T;
            H:=-T * A[K, K];
            For J:=K + 1 To N do
            begin
                F:=0;
                For I:=K To M do
                    F:=F + A[I, K] * A[I, J];
                F:=F / H;
                For I:=K To M do
                    A[I, J]:=A[I, J] - A[I, K] * F;
            end;
            For J:=1 To M do
            begin
                F:=0;
                For I:=K To M do
                    F:=F + A[I, K] * Q[I, J];
                F:=F / H;
                For I:=K To M do
                    Q[I, J]:=Q[I, J] - A[I, K] * F;
            end;
            A[K, K]:=T * S;
            For I:=K + 1 To M do
                A[I, K]:=0;
        end;
    end;
end;
end.