matrix.~pas

以面向对象方法实现的数值算法类库

{-------------------------------------------------------------------------------
2004.07.02 定义矩阵类
  定义矩阵的基本运算
--------------------------------------------------------------------------------
2004.07.03
  定义矩阵之间的乘法：矩阵相乘
  矩阵相乘以后的结果：目前相乘后的结果是在相乘函数中动态创建的对象，但是存在内存
释放的问题。如何解决？--可以人工释放内存，调用Free方法--以前所作多相流数据处理程
序存在同样的问题。即：两个向量相乘或求互相关、自相关的时候。
-------------------------------------------------------------------------------}

unit Matrix;

interface

uses classes, MathTypes, DB, Series, Grids;

Type

  TMatrix = class
  protected
    {--------------------------------------------------------------------------}
    FData:          TFloatMatrix;              // 数据
    FRowCount:      integer;                   // 矩阵的大小 行数
    FColCount:      integer;                   // 矩阵的大小 列数
    {--------------------------------------------------------------------------}
    FFormat:        string;                    // 显示格式
    FDisplayGrid:   TStringGrid;               // 显示器
    FStartRow:      Integer;
    FStartCol:      Integer;
    {--------------------------------------------------------------------------}
    FSeries:        TCustomSeries;             // 曲线
    FXCol:          Integer;
    FYCol:          Integer;
    {--------------------------------------------------------------------------}
    FScope:         integer;                   // 取值范围，用于填充随机数
    FPause:         Boolean;                   // 消元过程暂停
    {--------------------------------------------------------------------------}
    FDataSet:       TDataSet;                  // 数据库
    FStartField:    Integer;                   // 首字段ID
    FFieldCount:    Integer;                   // 字段数，如果为0,则根据列数从首字段依次读入  
    {--------------------------------------------------------------------------}
  public
    { 构造与析构方法-----------------------------------------------------------}
    constructor Create;
    destructor  Destroy; override;
    { 数据赋值 ----------------------------------------------------------------}
    procedure   CopyFrom(Source: TMatrix);
    procedure   CopyFromArray(x: TFloatMatrix);
    { 数据维护 ----------------------------------------------------------------}
    procedure   Resize(RowCount, ColCount: integer);
    {--------------------------------------------------------------------------}
    procedure   LoadFromFile(FileName: string);
    procedure   SaveToFile(FileName: string);
    procedure   LoadFromDatabase(ADataSet: TDataSet; FieldID: integer = 0;
      nFieldCount: integer = 1; Start: integer = 0; nCount: integer = -1); overload;
    procedure   LoadFromDatabase; overload;
    {--------------------------------------------------------------------------}
    procedure   ShowInGrid;
    procedure   ShowAsSeries;
    {--------------------------------------------------------------------------}
    procedure   FillSample;
    procedure   FillAsUnit;
    {--------------------------------------------------------------------------}
    function    CombineCol(bMat: TMatrix): TMatrix;          // 合并两个矩阵
    procedure   CombineMatrix(bMatrix, AResult: TMatrix);
    { 常见的数值方法 ----------------------------------------------------------}
    procedure   GassJordan_CompletePivot; overload;
    procedure   GassJordan_CompletePivot(bMatrix: TMatrix); overload;
    procedure   GaussJordan(bMatrix: TMatrix);
    procedure   MultiPlyMatrix(BMat, AResult: TMatrix);
    function    MultiPly(BMat: TMatrix): TMatrix;

  published
    property ColCount: integer read FColCount;
    property RowCount: integer read FRowCount;
    property Values: TFloatMatrix read FData write FData;
    property DisplayFormat: string read FFormat write FFormat;
    property SampleScope: Integer read FScope write FScope;

    property DisplayGrid: TStringGrid read FDisplayGrid write FDisplayGrid;
    property StartRow: Integer read FStartRow write FStartRow;
    property StartCol: Integer read FStartCol write FStartCol;

    property ASeries: TCustomSeries read FSeries write FSeries;
    property XCol: integer read FXCol write FXCol;
    property YCol: integer read FYCol write FYCol;

    property Pause: Boolean read FPause write FPause;

  end;


implementation

Uses SysUtils,Dialogs;

{ TMatrix }

procedure TMatrix.CopyFrom(Source: TMatrix);
Var
    i,j: Integer;
begin
    inherited;
    if Source is TMatrix then begin
      Resize((Source as TMatrix).RowCount,(Source as TMatrix).ColCount);
      for i:=0 to FRowCount-1 do
        for j:=0 to FColCount-1 do
          FData[i][j] := (Source as TMatrix).FData[i][j];
    end;
end;

procedure TMatrix.CopyFromArray(x: TFloatMatrix);
Var
    i,j,M,N: Integer;
begin
    M := Length(x);
    N := Length(x[0]);
    Resize(M, N);
    for i:=0 to FRowCount-1 do
      for j:=0 to FColCount-1 do
        FData[i][j] := x[i][j];
end;

function TMatrix.CombineCol(bMat: TMatrix): TMatrix;
Var
    M,N,i,j: integer;
begin
    M := FRowCount;
    if M<>bMat.RowCount then begin
      Raise
        EMathError.CreateFmt('行数不相等(%d,%d)=0',[M, bMat.FRowCount]);
    end;
    N := FColCount + bMat.ColCount;
    Result := TMatrix.Create;
    Result.Resize(M,N);
    for i:=0 to M-1 do
      for j:=0 to N-1 do
        if j<FColCount then Result.FData[i][j] := FData[i][j]
        else                Result.FData[i][j] := bMat.FData[i][j-FColCount];
end;

constructor TMatrix.Create;
begin
    inherited Create;
    {--------------------------------------------------------------------------}
    FData := Nil;
    FColCount := 0;
    FRowCount := 0;
    {--------------------------------------------------------------------------}
    FFormat := '%10.4f';
    FDisplayGrid := Nil;
    FStartRow := 1;
    FStartCol := 1;
    {--------------------------------------------------------------------------}
    FSeries := Nil;
    FXCol := 0;
    FYCol := 1;
    {--------------------------------------------------------------------------}
    FPause := False;
    FScope := 100;
    {--------------------------------------------------------------------------}
    FDataSet := Nil;
    FStartField := 0;
    FFieldCount := 0;
end;

destructor TMatrix.Destroy;
Var
    i: Integer;
begin
    for i:=0 to FRowCount-1 do FData[i] := Nil;
    FData := Nil;
    inherited;
end;

procedure TMatrix.FillAsUnit;
Var
    i,j: integer;
begin
    for i:=0 to FRowCount-1 do
      for j:=0 to FColCount-1 do
        if i=j then FData[i][j] := 1
        else        FData[i][j] := 0;
end;

procedure TMatrix.FillSample;
Var
    i,j: integer;
begin
    Randomize;
    for i:=0 to FRowCount-1 do
      for j:=0 to FColCount-1 do
        FData[i][j] := Random(FScope);
end;

{  高斯-若当全主元消去法,针对单一矩阵------------------------------------------}
procedure TMatrix.GassJordan_CompletePivot(bMatrix: TMatrix);
Var
    i,j,ii,jj,iP,jP: Integer;
    Big, aPiv, invp, tmp: double;
    iDone: array of Boolean;
    jDone: array of Boolean;
begin
    SetLength(iDone, FColCount);
    SetLength(jDone, FColCount);
    for i:=0 to FColCount-1 do iDone[i] := False;
    for i:=0 to FColCount-1 do jDone[i] := False;
    for i:=0 to FRowCount-1 do begin
      // 选主元--从(0,0)==>(m,m),跳过已经消元的列
      Big := 0;
      iP := 0;
      jP := 0;
      for ii:=0 to FRowCount-1 do begin
        if not iDone[ii] then begin
          for jj:=0 to FRowCount-1 do begin
            if not jDone[jj]then begin
              if Big < abs(FData[ii][jj]) then begin
                // 标记主元
                iP := ii; jP := jj; Big := Abs(FData[ii][jj]);
              end;
            end;
          end;
        end;
      end;
      // 奇异性判断
      if Big=0 then begin
        Raise
          EMathError.CreateFmt('奇异矩阵(%d,%d)=0',[iP, jP]);
      end;
      iDone[iP] := True;
      jDone[jP] := True;
      aPiv := FData[iP][jP];    // 主元
      // 将主元所在行除以主元，将主元归一
      for jj:=0 to FColCount-1 do begin
        if jj<>jP then FData[iP][jj] := FData[iP][jj] / aPiv
        else           FData[iP][jj] := 1;
      end;
      // 消元--将主元所在列的其他元素消成0;消元的过程共M行
      for ii:=0 to FRowCount-1 do begin
        if ii<>iP then begin  // 非主元所在行才需要进行消元
          if FData[ii][jP]<>0 then begin  // 与主元同列的元素不为0，需要消元
            invp := FData[ii][jP];
            for jj:=0 to FColCount-1 do begin
              if jj<>jP then FData[ii][jj] := FData[ii][jj] - FData[iP][jj] * invP
              else           FData[ii][jj] := 0;
            end;
          end;
        end;
      end;// 消元过程结束
      if FPause then begin
        ShowInGrid;
        ShowMessage(Format('%d,%d',[iP,jP]));
      end;
    end;
    // 整理，将结果按顺序排列好
    iP := 0;
    for j:=0 to FRowCount-1 do begin
      // 查找该列的主元
      for i:=0 to FRowCount do begin
        if FData[i][j]=1 then begin
          iP := i;
          break;
        end;
      end;
      // 交换
      if iP<>j then begin
        for jj:=0 to FColCount-1 do begin
          tmp := FData[j][jj];
          FData[j][jj] := FData[iP][jj];
          FData[iP][jj] := tmp;
        end;
      end;
    end;
end;

{  高斯-若当全主元消去法,针对双矩阵-------------------------------------------}
procedure TMatrix.GassJordan_CompletePivot(bMatrix: TMatrix);
Var
    i,j,ii,jj,iP,jP: Integer;
    Big, aPiv, invp, tmp: double;
    iDone: array of Boolean;
    jDone: array of Boolean;
begin
    SetLength(iDone, FColCount);
    SetLength(jDone, FColCount);
    for i:=0 to FColCount-1 do iDone[i] := False;
    for i:=0 to FColCount-1 do jDone[i] := False;
    for i:=0 to FRowCount-1 do begin
      // 选主元--从(0,0)==>(m,m),跳过已经消元的列
      Big := 0;
      iP := 0;
      jP := 0;
      for ii:=0 to FRowCount-1 do begin
        if not iDone[ii] then begin
          for jj:=0 to FRowCount-1 do begin
            if not jDone[jj]then begin
              if Big < abs(FData[ii][jj]) then begin
                // 标记主元
                iP := ii; jP := jj; Big := Abs(FData[ii][jj]);
              end;
            end;
          end;
        end;
      end;