{-------------------------------------------------------------------------------
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;