matrix.~pas

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

      // 奇异性判断
      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;
      // 处理第二矩阵
      for jj:=0 to bMatrix.FColCount-1 do
        bMatrix.FData[iP][jj] := bMatrix.FData[iP][jj]  / aPiv;
      // 消元--将主元所在列的其他元素消成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;
            // 处理第二矩阵
            for jj:=0 to bMatrix.FColCount-1 do
              bMatrix.FData[ii][jj] :=
                bMatrix.FData[ii][jj] - bMatrix.FData[iP][jj] * invP;
          end;
        end;
      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;
        // 处理第二矩阵
        for jj:=0 to bMatrix.FColCount-1 do begin
          tmp := bMatrix.FData[j][jj];
          bMatrix.FData[j][jj] := bMatrix.FData[iP][jj];
          bMatrix.FData[iP][jj] := tmp;
        end;
      end;
    end;
end;

procedure TMatrix.LoadFromFile(FileName: string);
Var
    sf: TextFile;
    i,j,M,N: integer;
begin
    AssignFile(sf,FileName);
    Reset(sf);
    Readln(sf, M, N);
    Resize(M, N);
    for i:=0 to M-1 do begin
      for j:=0 to N-1 do read(sf, FData[i][j]);
      readln(sf);
    end;
    CloseFile(sf);
end;

function TMatrix.MultiPly(BMat: TMatrix): TMatrix;
Var
    M, L, N: integer;
    i,j,k: integer;
    v: double;
begin
    M := RowCount;
    L := ColCount;
    if L<>BMat.FRowCount then
      Raise
        EMathError.CreateFmt('矩阵的维数不匹配：A(%d,%d),B(%d,%d)',
          [M, L, BMat.RowCount, BMat.ColCount]);
    N := BMat.ColCount;
    Result := TMatrix.Create;
    Result.Resize(M,N);
    for i:=0 to M-1 do begin
      for j:=0 to N-1 do begin
        v := 0;
        for k:=0 to L-1 do v := v + FData[i][k] * BMat.Values[k][j];
        Result.FData[i][j] := v;
      end;
    end;
end;

procedure TMatrix.Resize(RowCount, ColCount: integer);
Var
    i: Integer;
begin
    if (RowCount<>FRowCount) or (FColCount<>ColCount) then begin
      for i:=0 to FRowCount-1 do FData[i] := Nil;
      FData := Nil;
      FRowCount := RowCount;
      FColCount := ColCount;
      SetLength(FData, FRowCount);
      for i:=0 to FRowCount-1 do SetLength(FData[i], FColCount);
    end;
end;

procedure TMatrix.SaveToFile(FileName: string);
Var
    sf: TextFile;
    i,j: integer;
begin
    AssignFile(sf,FileName);
    Rewrite(sf);
    Writeln(sf, FRowCount:10, FColCount:10);
    for i:=0 to FRowCount-1 do begin
      for j:=0 to FColCount-1 do write(sf, format(FFormat,[FData[i][j]]));
      writeln(sf);
    end;
    CloseFile(sf);
end;

procedure TMatrix.ShowInGrid;
Var
    i,j: Integer;
begin
    if FDisplayGrid<>Nil then begin
      for i:=0 to FRowCount-1 do
        for j:=0 to FColCount-1 do
          FDisplayGrid.Cells[j+StartCol,i+StartRow] := Format(FFormat,[FData[i][j]]);
    end;
end;

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

procedure TMatrix.GaussJordan(bMatrix: TMatrix);
Var
    i,j,ii: integer;
    jP: integer;
    big: double;
    temp: double;
    av: double;
    ak: double;
begin
    for i:=0 to FRowCount-1 do begin
      // 选主元
      big := 0;
      jP := 0;
      for ii:=i to FRowCount-1 do begin
        if big < abs(FData[ii][i]) then begin
          big := abs(FData[ii][i]);
          jP := ii;
        end;
      end;
      // 调换
      if jP<>i then begin
        // 第一
        for j:=0 to FColCount-1 do begin
          temp := FData[jP][j];
          FData[jP][j] := FData[i][j];
          FData[i][j] := temp;
        end;
        // 第二
        for j:=0 to bMatrix.FColCount-1 do begin
          temp := bMatrix.FData[jP][j];
          bMatrix.FData[jP][j] := bMatrix.FData[i][j];
          bMatrix.FData[i][j] := temp;
        end;
      end;
      // 归一
      av := FData[i][i];
      // 第一矩阵
      FData[i][i] := 1;
      for j:=i+1 to FColCount-1 do begin
        FData[i][j] := FData[i][j] / av;
      end;
      // 第二矩阵
      for j:=0 to bMatrix.FColCount-1 do
        bMatrix.FData[i][j] := bMatrix.FData[i][j] / av;
      // 消元-
      for ii:=0 to FRowCount-1 do begin
        if ii<>i then begin
          ak := FData[ii][i];
          FData[ii][i] := 0;
          // 第一矩阵
          for j:=i+1 to FColCount-1 do
            FData[ii][j] := FData[ii][j] - FData[i][j] * ak;
          // 第二矩阵
          for j:=0 to bMatrix.FColCount-1 do
            bMatrix.FData[ii][j] := bMatrix.FData[ii][j] - bMatrix.FData[i][j] * ak;
        end;
      end;
    end;
end;

procedure TMatrix.MultiPlyMatrix(BMat, AResult: TMatrix);
Var
    M, L, N: integer;
    i,j,k: integer;
    v: double;
begin
    M := RowCount;
    L := ColCount;
    if L<>BMat.FRowCount then
      Raise
        EMathError.CreateFmt('矩阵的维数不匹配：A(%d,%d),B(%d,%d)',
          [M, L, BMat.RowCount, BMat.ColCount]);
    N := BMat.ColCount;
    if Aresult<>Nil then begin
      AResult.Resize(M,N);
      for i:=0 to M-1 do begin
        for j:=0 to N-1 do begin
          v := 0;
          for k:=0 to L-1 do v := v + FData[i][k] * BMat.Values[k][j];
          AResult.FData[i][j] := v;
        end;
      end;
    end;
end;

procedure TMatrix.LoadFromDatabase(ADataSet: TDataSet; FieldID, nFieldCount,
  Start, nCount: integer);
Var
    i,j,N: integer;
begin
    if nCount<0 then nCount := ADataSet.RecordCount;
    Resize(nCount, nFieldCount);
    N := nCount;
    ADataSet.FindFirst;
    for i:=0 to Start-1 do ADataSet.Next;
    if nCount<0 then begin
      for i:=0 to N-1 do
        if NOT ADataSet.Eof then begin
          for j:=0 to nFieldCount-1 do
            FData[i][j] := ADataSet.Fields[FieldId+j].AsFloat;
          ADataSet.Next;
        end else break;
    end else begin
      for i:=0 to nCount-1 do
        if NOT ADataSet.Eof then begin
          for j:=0 to nFieldCount-1 do
            FData[i][j] := ADataSet.Fields[FieldId+j].AsFloat;
          ADataSet.Next;
        end else break;
    end;
end;

procedure TMatrix.ShowAsSeries;
Var
    i: integer;
begin
    if FSeries<>Nil then begin
      FSeries.Clear;
      for i:=0 to FRowCount-1 do FSeries.AddXY(FData[i][FXCol], FData[i][FYCol]);
    end;
end;

procedure TMatrix.LoadFromDatabase;
Var
    i, j, M, N: integer;
begin
    if FDataSet<>Nil then begin
      M := FDataSet.RecordCount;
      if FFieldCount<1 then N := FColCount else N := fFieldCount;
      Resize(M, N);
      FDataSet.FindFirst;
      for i:=0 to M-1 do begin
        for j:=0 to N-1 do begin
          FData[i][j] := FDataSet.Fields[FStartField+j].AsFloat;
        end;
        FDataSet.Next;
      end;
    end;
end;

end.