简单的矩阵计算(结构体,静态数组)

有个需求需要用到矩阵工具..从十几年前的资料里找到了一个矩阵计算,但是是用类来实现的,改写了一下,改成结构体.只做了简单测试,发现问题欢迎告诉我...

  1 unit MatrixCal;
  2 
  3 interface
  4 
  5 uses
  6   Windows, SysUtils;
  7 
  8 const
  9   MaxRowCount = 8;
 10   MaxColCount = 8;
 11   
 12 type
 13   EMatrixException = class(Exception);
 14   
 15 type
 16   TMatrixDataAry = array[0..MaxRowCount - 1] of array[0..MaxColCount - 1] of Extended;
 17   TMatrixRecord = record
 18     RowCount, ColCount: Integer;
 19     Data: TMatrixDataAry;
 20   end;
 21 
 22 procedure Matrix_Init(const RowCount, ColCount: Integer; Value: Extended; var M: TMatrixRecord);
 23 
 24 
 25 procedure Matrix_MainBasis(Value: Extended; var M: TMatrixRecord);
 26 function  Matrix_SolveEx(var A, B: TMatrixRecord): Extended;
 27 
 28 
 29 function  Matrix_CalMultValue(M: TMatrixRecord; Value: Extended): TMatrixRecord;
 30 function  Matrix_CalMult(A, B: TMatrixRecord; var C: TMatrixRecord): Boolean;
 31 function  Matrix_CalAdd(A, B: TMatrixRecord; var C: TMatrixRecord): Boolean;
 32 function  Matrix_CalPower(M: TMatrixRecord; Pow: Integer): TMatrixRecord;
 33 function  Matrix_CalDeterminant(M: TMatrixRecord; var Value: Extended): Boolean;
 34 function  Matrix_CalInverse(M: TMatrixRecord; var B: TMatrixRecord): Boolean;
 35 
 36 implementation
 37 
 38 procedure Matrix_Init(const RowCount, ColCount: Integer; Value: Extended; var M: TMatrixRecord);
 39 var
 40   i, j: Integer;
 41 begin
 42   M.RowCount := RowCount;
 43   M.ColCount := ColCount;
 44   for i := 0 to RowCount - 1 do
 45     for j := 0 to ColCount - 1 do
 46       M.Data[i, j] := Value;
 47 end;
 48 
 49 
 50 procedure Matrix_MainBasis(Value: Extended; var M: TMatrixRecord);
 51 var
 52   i: Integer;
 53 begin
 54   if M.RowCount <> M.ColCount then
 55     raise EMatrixException.Create('This operation can be applied only to square matrix');
 56   Matrix_Init(M.RowCount, M.ColCount, 0, M);
 57   for i := 0 to M.RowCount - 1 do
 58     M.Data[i, i] := Value;
 59 end;
 60 
 61 procedure Matrix_ExChangeRows(i, j: Integer; var M: TMatrixRecord);
 62 var
 63   k: Integer;
 64   Temp: Extended;
 65 begin
 66   for k := 0 to M.ColCount - 1 do
 67   begin
 68     Temp := M.Data[i,k];
 69     M.Data[j,k] := M.Data[i,k];
 70     M.Data[j,k] := Temp;
 71   end;
 72 end;
 73 
 74 procedure Matrix_AddRows(i, j: integer; Factor: Extended; var M: TMatrixRecord);
 75 var
 76   k: Integer;
 77 begin
 78   for k := 0 to M.ColCount - 1 do
 79     M.Data[i,k] := M.Data[i,k] + M.Data[j,k] * Factor;
 80 end;
 81 
 82 procedure Matrix_MultLine(i: Integer; Factor: Extended; var M: TMatrixRecord);
 83 var
 84   k: Integer;
 85 begin
 86   for k := 0 to M.ColCount - 1 do
 87     M.Data[i,k] := M.Data[i,k] * Factor;
 88 end;
 89 
 90 function  Matrix_SolveEx(var A, B: TMatrixRecord): Extended;
 91 var
 92   i, j: Integer;
 93   Temp: Extended;
 94   Factor: Integer;
 95 begin
 96   Result := 1;
 97   { Forward pass, choose main element <> 0 }
 98   for i := 0 to A.RowCount - 1 do
 99   begin
100     j := 0;
101     while j < A.RowCount do
102       if A.Data[i, j] <> 0 then
103         Break
104       else
105         Inc(j);
106     if j = A.RowCount then
107     begin
108       Result := 0;
109       Exit;
110     end;
111     Matrix_ExChangeRows(i, j, A);
112     Matrix_ExChangeRows(i, j, B);
113     if not Odd(A.ColCount) then Result := -1 * Result;
114     for j := i + 1 to A.RowCount - 1 do
115     begin
116       Temp := A.Data[j, i]/A.Data[i, i];
117       Matrix_AddRows(j, i, -Temp, A);
118       Matrix_AddRows(j, i, -Temp, B);
119     end;
120   end;
121 
122   { Backward pass }
123   for i := A.RowCount - 1 downto 0 do
124     for j := i - 1 downto 0 do
125     begin
126       Temp := A.Data[j, i]/A.Data[i, i];
127       Matrix_AddRows(j, i, -Temp, A);
128       Matrix_AddRows(j, i, -Temp, B);
129     end;
130 
131   { Divide diagonal elements by themself, making identity matrix }
132   for i := 0 to A.RowCount - 1 do
133   begin
134     Result := Result * A.Data[i, i];
135     Temp := 1/A.Data[i, i];
136     Matrix_MultLine(i, Temp, B);
137     Matrix_MultLine(i, Temp, A);
138   end;
139 end;
140 
141 function  Matrix_CalMultValue(M: TMatrixRecord; Value: Extended): TMatrixRecord;
142 var
143   i, j: Integer;
144 begin
145   Result := M;
146   for i := 0 to Result.RowCount - 1 do
147     for j := 0 to Result.ColCount - 1 do
148       Result.Data[i,j] := Result.Data[i,j] * Value;
149 end;
150 
151 function  Matrix_CalMult(A, B: TMatrixRecord; var C: TMatrixRecord): Boolean;
152 var
153   i, j, k: Integer;
154   TmpC: TMatrixRecord;
155 begin
156   Result := False;
157   if A.ColCount <> B.RowCount then
158     Exit; //raise EMatrixException.Create('Can not Multiply these matrixes');
159   Matrix_Init(A.RowCount, B.ColCount, 0, TmpC);
160   for k := 0 to A.RowCount - 1 do
161     for i := 0 to B.ColCount - 1 do
162       for j := 0 to A.ColCount - 1 do
163         TmpC.Data[k, i] := TmpC.Data[k, i] + A.Data[k, j] * B.Data[j, i];
164   Result := True;
165   C := TmpC;
166 end;
167 
168 function  Matrix_CalAdd(A, B: TMatrixRecord; var C: TMatrixRecord): Boolean;
169 var
170   i, j, k: Integer;
171   TmpC: TMatrixRecord;
172 begin
173   Result := False;
174   if (A.ColCount <> B.ColCount) or (A.RowCount <> B.RowCount) then
175     Exit; //raise EMatrixException.Create('Can not add these matrixes');
176   TmpC := A;
177   for i := 0 to A.RowCount - 1 do
178     for j := 0 to A.ColCount - 1 do
179         TmpC.Data[i, j] := TmpC.Data[i, j] + B.Data[i, j];
180   Result := True;
181   C := TmpC;
182 end;
183 
184 function  Matrix_CalPower(M: TMatrixRecord; Pow: Integer): TMatrixRecord;
185 var
186   i: Integer;
187 begin
188   Result := M;
189   for i := 1 to Pow - 1 do
190     Matrix_CalMult(Result, M, Result);
191 end;
192 
193 function  Matrix_CalDeterminant(M: TMatrixRecord; var Value: Extended): Boolean;
194 var
195   B: TMatrixRecord;
196 begin
197   Result := False;
198   if M.RowCount <> M.ColCount then
199     Exit;  //raise EMatrixException.Create('This operation can be applied only to square matrix');
200   B := M;
201   Matrix_MainBasis(1, B);
202   Value := Matrix_SolveEx(M, B);
203   Result := True;
204 end;
205 
206 function  Matrix_CalInverse(M: TMatrixRecord; var B: TMatrixRecord): Boolean;
207 var
208   TmpB: TMatrixRecord;
209 begin
210   Result := False;
211   if M.RowCount <> M.ColCount then
212     Exit; //raise EMatrixException.Create('This operation can be applied only to square matrix');
213   TmpB := M;
214   Matrix_MainBasis(1, TmpB);
215   Result := Matrix_SolveEx(M, TmpB) <> 0;
216   if Result then
217     B := TmpB;
218 end;
219 
220 end.

View Code

同时附上用类实现的代码.作者不可考.

   1 unit MatLib;
   2 
   3 interface
   4 
   5 uses
   6   SysUtils, WinTypes, WinProcs, Messages, Classes, Graphics, Controls,
   7   Forms, Dialogs, StdCtrls;
   8 
   9 const
  10   { Size of Extended }
  11   szExtended = SizeOf(Extended);
  12   szPointer = SizeOf(Pointer);
  13 
  14 type
  15   { Array of Extended }
  16   PExtArray = ^TExtArray;
  17   TExtArray = array [0..MaxInt div szExtended - 1] of Extended;
  18   { Array of pointers }
  19   PPtrArray = ^TPtrArray;
  20   TPtrArray = array [0..MaxInt div szPointer - 1] of PExtArray;
  21 
  22 type
  23   EMatrixException = class(Exception);
  24 
  25   TMatrix = class;
  26   TVector = class;
  27   TMatrixClass = class of TMatrix;
  28   TVectorClass = class of TVector;
  29 
  30   { TMatrix - base abstract class for matrices }
  31   TMatrix = class(TComponent)
  32   protected
  33     FColCount: Integer;
  34     FRowCount: Integer;
  35     FClass: TMatrixClass;
  36     { Must be overrided in decendent class !}
  37     function  GetCells(ARow, ACol: Integer): Extended; virtual; abstract;
  38     { Must be overrided in decendent class !}
  39     procedure SetCells(ARow, ACol: Integer; AValue: Extended); virtual; abstract;
  40     { Must be overrided in decendent class !}
  41     procedure SetRanges(NewRow, NewCol: Integer); virtual; abstract;
  42     procedure SetCols(Value: Integer);
  43     procedure SetRows(Value: Integer);
  44   public
  45     constructor Create(AOwner: TComponent); override;
  46     { Assign matrix from Source }
  47     procedure   Assign(Source: TMatrix);
  48     { Assign Value to all cells of matrix }
  49     procedure   AssignValue(Value: Extended);
  50     { Assign Values array to NLine row of matrix  }
  51     procedure   AssignToLine(NLine: Integer; Values: array of Extended);
  52     { Assign Values array to NCol column of matrix  }
  53     procedure   AssignToCol (NCol : Integer; Values: array of Extended);
  54     { Add matrix AMatrix to Self (this matrix) }
  55     procedure   Add(AMatrix: TMatrix);
  56     { Multiply Self by AMatrix }
  57     procedure   MultAsLeft(AMatrix: TMatrix);
  58     { Multiply Self by AMatrix }
  59     procedure   MultAsRight(AMatrix: TMatrix);
  60     { Multiply Self by Value }
  61     procedure   MultValue(Value: Extended);
  62     { Assign Value to diagonal elements (only square matrix) }
  63     procedure   MainBasis(Value: Extended);
  64     { Calculate determinant }
  65     function    Determinant: Extended;
  66     { Inverses matrix  }
  67     procedure   Inverse;
  68     { Transpose matrix }
  69     procedure   Transpose;
  70     { ExChanges rows (Rows) number i and j }
  71     procedure   ExChangeRows(i, j: Integer);
  72     { Add to i-row j-row muliplied by Factor }
  73     procedure   AddRows(i, j: integer; Factor: Extended);
  74     { Multiply all elements of i-row by Factor }
  75     procedure   MultLine(i: Integer; Factor: Extended);
  76     { Raise to power Pow  }
  77     procedure   Power(Pow: Integer);
  78     property Cells[ARow, ACol: Integer]: Extended read GetCells write SetCells; default;
  79     property ColCount: Integer read FColCount write SetCols;
  80     property RowCount: Integer read FRowCount write SetRows;
  81   end;
  82 
  83   { Matrix of Extended [(MaxInt div 4)x(MaxInt div 10) elements maximum) }
  84   TRealMatrix = class(TMatrix)
  85   protected
  86     FRows: PPtrArray;
  87     destructor  Destroy; override;
  88     function  GetCells(ARow, ACol: Integer): Extended; override;
  89     procedure SetCells(ARow, ACol: Integer; AValue: Extended); override;
  90     procedure SetRanges(NewRow, NewCol: Integer); override;
  91   public
  92     property Cells;
  93   published
  94     property ColCount;
  95     property RowCount;
  96   end;
  97 
  98   { Sparse matrix cell }
  99   PSCell = ^TSCell;
 100   TSCell = record
 101     Col, Row: Integer;  { location of this cell in the matrix }
 102     Value: Extended;    { value to hold in the cell }
 103   end;
 104 
 105   { Sparse matrix row }
 106   TSRow = class(TList)
 107   private
 108     FRowNum: Integer;
 109     function Get(Index: Integer): PSCell;
 110     procedure Put(Index: Integer; ACell: PSCell);
 111   public
 112     constructor Create(ARowNum: Integer);
 113     destructor  Destroy; override;
 114     property RowNum: Integer read FRowNum write FRowNum;
 115     property Items[Index: Integer]: PSCell read Get write Put;
 116     function Find(ACellCol: Integer; var Index: Integer): Boolean;
 117     function Add(const ACellCol: Integer; const AValue: Extended): Integer;
 118   end;
 119 
 120   { Sparse matrix row list }
 121   TRowList = class(TList)
 122   private
 123     function Get(Index: Integer): TSRow;
 124     procedure Put(Index: Integer; ARow: TSRow);
 125   public
 126     destructor Destroy; override;
 127     property Items[Index: Integer]: TSRow read Get write Put;
 128     function Find(ARowNum: Integer; var Index: Integer): Boolean;
 129     function Add(const ARowNum: Integer): Integer;
 130   end;
 131 
 132   { Sparse matrix }
 133   TSparseMatrix = class(TMatrix)
 134   protected
 135     FRows: TRowList;
 136     destructor  Destroy; override;
 137     function  GetCells(ARow, ACol: Integer): Extended; override;
 138     procedure SetCells(ARow, ACol: Integer; AValue: Extended); override;
 139     procedure SetRanges(NewRow, NewCol: Integer); override;
 140   public
 141     property Cells;
 142   published
 143     property ColCount;
 144     property RowCount;
 145   end;
 146 
 147   { TVector - base abstract class for vectors }
 148   TVector = class(TComponent)
 149   protected
 150     FCount: Integer;
 151     FClass: TVectorClass;
 152     { Must be overrided in decendent class !}
 153     function  GetCells(APos: Integer): Extended; virtual; abstract;
 154     { Must be overrided in decendent class !}
 155     procedure SetCells(APos: Integer; AValue: Extended); virtual; abstract;
 156     { Must be overrided in decendent class !}
 157     procedure SetCount(NewCount: Integer); virtual; abstract;
 158   public
 159     constructor Create(AOwner: TComponent); override;
 160     { Assign Vector from Source }
 161     procedure   Assign(Source: TVector);
 162     { Assign Value to all cells of Vector }
 163     procedure   AssignValue(Value: Extended);
 164     { Add Vector AVector to Self (this Vector) }
 165     procedure   Add(AVector: TVector);
 166     {}
 167     procedure   Exchange(i, j: Integer);
 168     { Multiply Self by AVector }
 169     procedure   MultValue(Value: Extended);
 170     { Scalar product }
 171     function    ScalarProduct(AVector: TVector): Extended;
 172     { Transpose Vector }
 173     procedure   Transpose;
 174     property Cells[APos: Integer]: Extended read GetCells write SetCells; default;
 175     property Count: Integer read FCount write SetCount;
 176   end;
 177 
 178   { Vector of Extended [(MaxInt div 10) elements maximum) }
 179   TRealVector = class(TVector)
 180   protected
 181     FData: PExtArray;
 182     destructor  Destroy; override;
 183     function  GetCells(APos: Integer): Extended; override;
 184     procedure SetCells(APos: Integer; AValue: Extended); override;
 185     procedure SetCount(NewCount: Integer); override;
 186   public
 187     property Cells;
 188   published
 189     property Count;
 190   end;
 191 
 192   { Sparse Vector }
 193   TSparseVector = class(TVector)
 194   protected
 195     FData: TSRow;
 196     destructor  Destroy; override;
 197     function  GetCells(APos: Integer): Extended; override;
 198     procedure SetCells(APos: Integer; AValue: Extended); override;
 199     procedure SetCount(NewCount: Integer); override;
 200   public
 201     property Cells;
 202   published
 203     property Count;
 204   end;
 205 
 206 { Routines which returnes matrices. A and B matrices will not ExChange }
 207 { !!! Note that you are responsible for destroying returned matrix.  }
 208 { type of returned matrix is the same as A                           }
 209 
 210 { }
 211 function SolveMatrix(AMatrix: TMatrix; BVector: TVector): Extended;
 212 { }
 213 function SolveMatrixEx(AMatrix: TMatrix; BMatrix: TMatrix): Extended;
 214 { Inverses matrix A }
 215 function M_Inverse(A: TMatrix): TMatrix;
 216 { MultAsLeftiply A * B  }
 217 function M_Mult(A, B: TMatrix): TMatrix;
 218 { Add A to B, A + B }
 219 function M_Add(A, B: TMatrix): TMatrix;
 220 { Mulltiply A by Value }
 221 function M_MultValue(A: TMatrix; Value: Extended): TMatrix;
 222 { Raise A to power Pow  }
 223 function M_Power(A: TMatrix; Pow: Integer): TMatrix;
 224 
 225 procedure Register;
 226 
 227 implementation
 228 
 229 procedure Register;
 230 begin
 231   RegisterComponents('Math', [TRealMatrix, TSparseMatrix, TRealVector, TSparseVector]);
 232 end;
 233 
 234 { SolveMatrix uses Gaussian algorithm to transform AMatrix to identity matrix }
 235 { At the same time it perform same operation on BMatrix, as a result it }
 236 { returns determinant of AMatrix }
 237 function SolveMatrixEx(AMatrix: TMatrix; BMatrix: TMatrix): Extended;
 238 var
 239   i, j: Integer;
 240   Temp: Extended;
 241   Factor: Integer;
 242 begin
 243   Result := 1;
 244   { Forward pass, choose main element <> 0 }
 245   for i := 0 to AMatrix.RowCount - 1 do
 246   begin
 247     j := 0;
 248     while j < AMatrix.RowCount do
 249       if AMatrix.Cells[i, j] <> 0 then break
 250                                   else inc(j);
 251     if j = AMatrix.RowCount then
 252     begin
 253       Result := 0;
 254       Exit;
 255     end;
 256     AMatrix.ExChangeRows(i, j);
 257     BMatrix.ExChangeRows(i, j);
 258     if not Odd(AMatrix.ColCount) then Result := -1 * Result;
 259     for j := i + 1 to AMatrix.RowCount - 1 do
 260     begin
 261       Temp := AMatrix.Cells[j, i]/AMatrix.Cells[i, i];
 262       AMatrix.AddRows(j, i, -Temp);
 263       BMatrix.AddRows(j, i, -Temp);
 264     end;
 265   end;
 266 
 267   { Backward pass }
 268   for i := AMatrix.RowCount - 1 downto 0 do
 269     for j := i - 1 downto 0 do
 270     begin
 271       Temp := AMatrix.Cells[j, i]/AMatrix.Cells[i, i];
 272       AMatrix.AddRows(j, i, -Temp);
 273       BMatrix.AddRows(j, i, -Temp);
 274     end;
 275 
 276   { Divide diagonal elements by themself, making identity matrix }
 277   for i := 0 to AMatrix.RowCount - 1 do
 278   begin
 279     Result := Result * AMatrix.Cells[i, i];
 280     Temp := 1/AMatrix.Cells[i, i];
 281     BMatrix.MultLine(i, Temp);
 282     AMatrix.MultLine(i, Temp);
 283   end;
 284 end;
 285 
 286 function SolveMatrix(AMatrix: TMatrix; BVector: TVector): Extended;
 287 var
 288   i, j: Integer;
 289   Temp: Extended;
 290   Factor: Integer;
 291 begin
 292   Result := 1;
 293   { Forward pass, choose main element <> 0 }
 294   for i := 0 to AMatrix.RowCount - 1 do
 295   begin
 296     j := 0;
 297     while j < AMatrix.RowCount do
 298       if AMatrix.Cells[i, j] <> 0 then break
 299                                   else inc(j);
 300     if j = AMatrix.RowCount then
 301     begin
 302       Result := 0;
 303       Exit;
 304     end;
 305     AMatrix.ExChangeRows(i, j);
 306     BVector.ExChange(i, j);
 307     if not Odd(AMatrix.ColCount) then Result := -1 * Result;
 308     for j := i + 1 to AMatrix.RowCount - 1 do
 309     begin
 310       Temp := AMatrix.Cells[j, i]/AMatrix.Cells[i, i];
 311       AMatrix.AddRows(j, i, -Temp);
 312       BVector[j] := BVector[j] - Temp*BVector[i];
 313     end;
 314   end;
 315 
 316   { Backward pass }
 317   for i := AMatrix.RowCount - 1 downto 0 do
 318     for j := i - 1 downto 0 do
 319     begin
 320       Temp := AMatrix.Cells[j, i]/AMatrix.Cells[i, i];
 321       AMatrix.AddRows(j, i, -Temp);
 322       BVector[j] := BVector[j] - Temp*BVector[i];
 323     end;
 324 
 325   { Divide diagonal elements by themself, making identity matrix }
 326   for i := 0 to AMatrix.RowCount - 1 do
 327   begin
 328     Result := Result * AMatrix.Cells[i, i];
 329     Temp := 1/AMatrix.Cells[i, i];
 330     BVector[i] := BVector[i] * Temp;
 331     AMatrix.MultLine(i, Temp);
 332   end;
 333 end;
 334 
 335 function Min(Num1, Num2: Integer): Integer;
 336 begin
 337   if Num1 < Num2 then Result := Num1
 338                  else Result := Num2;
 339 end;
 340 
 341 { -- TMatrix -- }
 342 constructor TMatrix.Create;
 343 begin
 344   inherited Create(AOwner);
 345   FClass := TMatrixClass(ClassType);
 346   SetRanges(5, 5);
 347 end;
 348 
 349 procedure TMatrix.Assign(Source: TMatrix);
 350 var
 351   i, j: Integer;
 352 begin
 353   RowCount := Source.RowCount;
 354   ColCount := Source.ColCount;
 355   for i := 0 to RowCount - 1 do
 356     for j := 0 to ColCount - 1 do
 357       Cells[i,j] := Source[i, j];
 358 end;
 359 
 360 procedure TMatrix.AssignValue(Value: Extended);
 361 var
 362   i, j: Integer;
 363 begin
 364   for i := 0 to RowCount - 1 do
 365     for j := 0 to ColCount - 1 do
 366       Cells[i,j] := Value;
 367 end;
 368 
 369 procedure TMatrix.AssignToLine(NLine: Integer; Values: array of Extended);
 370 var
 371   j, Min: Integer;
 372 begin
 373   if (NLine >= RowCount) or (NLine < 0)
 374     then raise EMatrixException.Create('Invalid row number');
 375   if ColCount > High(Values)
 376     then Min := High(Values)
 377     else Min := ColCount - 1;
 378   for j := 0 to Min do
 379     Cells[NLine,j] := Values[j];
 380 end;
 381 
 382 procedure TMatrix.AssignToCol(NCol: Integer; Values: array of Extended);
 383 var
 384   j, Min: Integer;
 385 begin
 386   if (NCol >= ColCount) or (NCol < 0)
 387     then raise EMatrixException.Create('Invalid column number');
 388   if RowCount > High(Values)
 389     then Min := High(Values)
 390     else Min := RowCount - 1;
 391   for j := 0 to Min do
 392     Cells[j,NCol] := Values[j];
 393 end;
 394 
 395 procedure TMatrix.Add(AMatrix: TMatrix);
 396 var
 397   i, j: Integer;
 398 begin
 399   if (ColCount <> AMatrix.ColCount) or (RowCount <> AMatrix.RowCount)
 400     then raise EMatrixException.Create('Can not add these matrixes');
 401   for i := 0 to RowCount - 1 do
 402     for j := 0 to ColCount - 1 do
 403       Cells[i,j] := Cells[i,j] + AMatrix[i,j];
 404 end;
 405 
 406 procedure TMatrix.MultAsLeft(AMatrix: TMatrix);
 407 var
 408   A: TMatrix;
 409   i, j, k: Integer;
 410 begin
 411   if ColCount <> AMatrix.RowCount
 412     then raise EMatrixException.Create('Can not Multiply these matrixes');
 413   A := FClass.Create(nil);
 414   A.RowCount := RowCount;
 415   A.ColCount := AMatrix.ColCount;
 416   for k := 0 to RowCount - 1 do
 417     for i := 0 to AMatrix.ColCount - 1 do
 418       for j := 0 to ColCount - 1 do
 419         A[k, i] := A[k, i] + Cells[k, j] * AMatrix.Cells[j, i];
 420   Assign(A);
 421   A.Free;
 422 end;
 423 
 424 procedure TMatrix.MultAsRight(AMatrix: TMatrix);
 425 var
 426   A: TMatrix;
 427   i, j, k: Integer;
 428 begin
 429   if ColCount <> AMatrix.RowCount
 430     then raise EMatrixException.Create('Can not Multiply these matrixes');
 431   A := FClass.Create(nil);
 432   A.RowCount := AMatrix.RowCount;
 433   A.ColCount := ColCount;
 434   for k := 0 to AMatrix.RowCount - 1 do
 435     for i := 0 to ColCount - 1 do
 436       for j := 0 to AMatrix.ColCount - 1 do
 437         A[k, i] := A[k, i] + AMatrix[k, j] * Cells[j, i];
 438   Assign(A);
 439   A.Free;
 440 end;
 441 
 442 procedure TMatrix.MultValue(Value: Extended);
 443 var
 444   i, j: Integer;
 445 begin
 446   for i := 0 to RowCount - 1 do
 447     for j := 0 to ColCount - 1 do
 448       Cells[i,j] := Cells[i,j] * Value;
 449 end;
 450 
 451 procedure TMatrix.MainBasis(Value: Extended);
 452 var
 453   i: Integer;
 454 begin
 455   if RowCount <> ColCount
 456     then raise EMatrixException.Create('This operation can be applied only to square matrix');
 457   AssignValue(0);
 458   for i := 0 to RowCount - 1 do
 459     Cells[i,i] := Value;
 460 end;
 461 
 462 function TMatrix.Determinant: Extended;
 463 var
 464   A, B: TMatrix;
 465 begin
 466   if RowCount <> ColCount
 467     then raise EMatrixException.Create('This operation can be applied only to square matrix');
 468   A := FClass.Create(nil);
 469   A.Assign(Self);
 470 
 471   B := FClass.Create(nil);
 472   B.RowCount := RowCount;
 473   B.ColCount := ColCount;
 474   B.MainBasis(1);
 475 
 476   Result := SolveMatrixEx(A, B);
 477 end;
 478 
 479 procedure TMatrix.Inverse;
 480 var
 481   A: TMatrix;
 482 begin
 483   if RowCount <> ColCount
 484     then raise EMatrixException.Create('This operation can be applied only to' +
 485                                        ' square matrix');
 486   A := FClass.Create(nil);
 487   A.RowCount := RowCount;
 488   A.ColCount := ColCount;
 489   A.MainBasis(1);
 490 
 491   if SolveMatrixEx(Self, A) = 0 then
 492   begin
 493     raise EMatrixException.Create('Could not inverse this matrix');
 494     Exit;
 495   end;
 496   Assign(A);
 497 end;
 498 
 499 procedure TMatrix.Transpose;
 500 var
 501   A: TMatrix;
 502   T: Extended;
 503   i, j: Integer;
 504 begin
 505   for i := 0 to RowCount - 1 do
 506     for j := i + 1 to ColCount - 1 do
 507     begin
 508     end;
 509 end;
 510 
 511 procedure TMatrix.Power(Pow: Integer);
 512 var
 513   Temp: TMatrix;
 514   i: Integer;
 515 begin
 516   Temp := FClass.Create(nil);
 517   Temp.Assign(Self);
 518   for i := 1 to Pow - 1 do
 519     MultAsLeft(Temp);
 520   Temp.Free;
 521 end;
 522 
 523 procedure TMatrix.SetCols(Value: Integer);
 524 begin
 525   try
 526     SetRanges(RowCount, Value);
 527   except
 528   end;
 529 end;
 530 
 531 procedure TMatrix.SetRows(Value: Integer);
 532 begin
 533   try
 534     SetRanges(Value, ColCount);
 535   except
 536   end;
 537 end;
 538 
 539 procedure TMatrix.ExChangeRows(i, j: Integer);
 540 var
 541   k: Integer;
 542   Temp: Extended;
 543 begin
 544   for k := 0 to ColCount - 1 do
 545   begin
 546     Temp := Cells[i,k];
 547     Cells[j,k] := Cells[i,k];
 548     Cells[j,k] := Temp;
 549   end;
 550 end;
 551 
 552 procedure TMatrix.AddRows(i, j: integer; Factor: Extended);
 553 var
 554   k: Integer;
 555 begin
 556   for k := 0 to ColCount - 1 do
 557     Cells[i,k] := Cells[i,k] + Cells[j,k] * Factor;
 558 end;
 559 
 560 procedure TMatrix.MultLine(i: Integer; Factor: Extended);
 561 var
 562   k: Integer;
 563 begin
 564   for k := 0 to ColCount - 1 do
 565     Cells[i,k] := Cells[i,k] * Factor;
 566 end;
 567 
 568 function M_Inverse(A: TMatrix): TMatrix;
 569 var
 570   B: TMatrix;
 571 begin
 572   if A.RowCount <> A.ColCount
 573     then raise EMatrixException.Create('This operation can be applied only to' +
 574                                        ' square matrix');
 575   B := A.FClass.Create(nil);
 576   B.Assign(A);
 577 
 578   Result := A.FClass.Create(nil);
 579   Result.RowCount := A.RowCount;
 580   Result.ColCount := A.ColCount;
 581   Result.MainBasis(1);
 582 
 583   if SolveMatrixEx(B, Result) = 0 then
 584   begin
 585     raise EMatrixException.Create('Could not inverse this matrix');
 586     Result := nil;
 587   end;
 588 end;
 589 
 590 function M_Mult(A, B: TMatrix): TMatrix;
 591 begin
 592   Result := A.FClass.Create(nil);
 593   Result.Assign(A);
 594   Result.MultAsLeft(B);
 595 end;
 596 
 597 function M_Add(A, B: TMatrix): TMatrix;
 598 begin
 599   Result := A.FClass.Create(nil);
 600   Result.Assign(A);
 601   Result.Add(B);
 602 end;
 603 
 604 function M_MultValue(A: TMatrix; Value: Extended): TMatrix;
 605 begin
 606   Result := A.FClass.Create(nil);
 607   Result.Assign(A);
 608   Result.MultValue(Value);
 609 end;
 610 
 611 function M_Power(A: TMatrix; Pow: Integer): TMatrix;
 612 begin
 613   Result := A.FClass.Create(nil);
 614   Result.Assign(A);
 615   Result.Power(Pow);
 616 end;
 617 
 618 { TRealMatrix }
 619 procedure TRealMatrix.SetRanges(NewRow, NewCol: Integer);
 620 var
 621   OldCol, OldRow: Integer;
 622   i, j: Integer;
 623 begin
 624   if (NewRow < 1) or (NewCol < 1)
 625     then raise Exception.Create('Invalid dimensions...');
 626   if (RowCount <> NewRow) or (ColCount <> NewCol)
 627   then begin
 628     OldRow := RowCount;
 629     OldCol := ColCount;
 630     { reallocate memory }
 631     { if OldCol < NewCol then new elements will be equal to 0 }
 632     { if OldCol > NewCol then elements will be lost }
 633     for i := 0 to OldRow - 1 do
 634     begin
 635     {$IFDEF WIN32}
 636       ReAllocMem(FRows^[i], szExtended * NewCol);
 637     {$ELSE}
 638       FRows^[i] := ReAllocMem(FRows^[i], szExtended * OldCol, szExtended * NewCol);
 639     {$ENDIF}
 640       if OldCol < NewCol
 641         then FillChar(FRows^[i]^[OldCol], (NewCol - OldCol)*szExtended, 0);
 642     end;
 643     { if NewRow < OldRow, unnessesary Rows will be destroed }
 644     for i := OldRow - 1 downto NewRow do
 645       FreeMem(FRows^[i], szExtended * OldCol);
 646     { Resize FRows }
 647     {$IFDEF WIN32}
 648     ReAllocMem(FRows, szPointer * NewRow);
 649     {$ELSE}
 650     FRows := ReAllocMem(FRows, szPointer * OldRow, szPointer * NewRow);
 651     {$ENDIF}
 652     { if NewRow > OldRows, new Rows will be added }
 653     for i := OldRow to NewRow - 1 do
 654     begin
 655       FRows^[i] := AllocMem(szExtended * NewCol);
 656       FillChar(FRows^[i]^, NewCol*szExtended, 0);
 657     end;
 658     { update FRowCount }
 659     FRowCount := NewRow;
 660     { update FColCount }
 661     FColCount := NewCol;
 662   end;
 663 end;
 664 
 665 function  TRealMatrix.GetCells(ARow, ACol: Integer): Extended;
 666 begin
 667   { if index is invalid then raise exception }
 668   if (ACol < 0) or (ACol > ColCount - 1) or
 669      (ARow < 0) or (ARow > RowCount - 1)
 670   then  raise EMatrixException.Create('Index out of bounds' + IntTostr(ACol) + ' ' + IntToStr(ARow));
 671   Result := FRows^[ARow]^[ACol];
 672 end;
 673 
 674 procedure TRealMatrix.SetCells(ARow, ACol: Integer; AValue: Extended);
 675 begin
 676   { if index is invalid then raise exception }
 677   if (ACol < 0) or (ACol > ColCount - 1) or
 678      (ARow < 0) or (ARow > RowCount - 1)
 679   then  raise EMatrixException.Create('Index out of bounds');
 680   FRows^[ARow]^[ACol] := AValue;
 681 end;
 682 
 683 destructor TRealMatrix.Destroy;
 684 var
 685   i: Integer;
 686 begin
 687   { deallocate memory }
 688   for i := 0 to RowCount - 1 do
 689     FreeMem(FRows^[i], szExtended * ColCount);
 690   inherited Destroy;
 691 end;
 692 
 693 { -- TSRow methods -- }
 694 constructor TSRow.Create(ARowNum: Integer);
 695 begin
 696   inherited Create;
 697   FRowNum := ARowNum;
 698 end;
 699 
 700 destructor TSRow.Destroy;
 701 var
 702   i: Integer;
 703 begin
 704   for i := 0 to Count - 1 do
 705     Dispose(PSCell(Items[i]));
 706   inherited Destroy;
 707 end;
 708 
 709 function TSRow.Get(Index: Integer): PSCell;
 710 begin
 711   Result := inherited Get(Index);
 712 end;
 713 
 714 procedure TSRow.Put(Index: Integer; ACell: PSCell);
 715 begin
 716   inherited Put(Index, ACell);
 717 end;
 718 
 719 function TSRow.Find(ACellCol: Integer; var Index: Integer): Boolean;
 720 var
 721   L, H, C, I: Integer;
 722 begin
 723   Result := False;
 724   L := 0;
 725   H := Count - 1;
 726   while L <= H do
 727   begin
 728     I := (L + H) shr 1;
 729     if Items[I]^.Col < ACellCol
 730       then L := I + 1
 731       else begin
 732         H := I - 1;
 733         if Items[i]^.Col = ACellCol then
 734         begin
 735           Result := True;
 736           L := I;
 737         end
 738       end;
 739   end;
 740   Index := L;
 741 end;
 742 
 743 function TSRow.Add(const ACellCol: Integer; const AValue: Extended): Integer;
 744 
 745 function NewCell(ARow, ACol: Integer; AValue: Extended): PSCell;
 746 begin
 747   New(Result);
 748   with Result^ do
 749   begin
 750     Row := ARow;
 751     Col := ACol;
 752     Value := AValue;
 753   end;
 754 end;
 755 
 756 begin
 757   Find(ACellCol, Result);
 758   Insert(Result, NewCell(FRowNum, ACellCol, AValue));
 759 end;
 760 
 761 { -- TRowList methods -- }
 762 destructor TRowList.Destroy;
 763 var
 764   i: Integer;
 765 begin
 766   for i := 0 to Count - 1 do
 767     Items[i].Free;
 768   inherited Destroy;
 769 end;
 770 
 771 function TRowList.Get(Index: Integer): TSRow;
 772 begin
 773   Result := inherited Get(Index);
 774 end;
 775 
 776 procedure TRowList.Put(Index: Integer; ARow: TSRow);
 777 begin
 778   inherited Put(Index, ARow);
 779 end;
 780 
 781 function TRowList.Find(ARowNum: Integer; var Index: Integer): Boolean;
 782 var
 783   L, H, C, I: Integer;
 784 begin
 785   Result := False;
 786   L := 0;
 787   H := Count - 1;
 788   while L <= H do
 789   begin
 790     I := (L + H) shr 1;
 791     if Items[i].RowNum < ARowNum
 792       then L := I + 1
 793       else begin
 794         H := I - 1;
 795         if Items[i].RowNum = ARowNum then
 796         begin
 797           Result := True;
 798           L := I;
 799         end
 800       end;
 801   end;
 802   Index := L;
 803 end;
 804 
 805 function TRowList.Add(const ARowNum: Integer): Integer;
 806 begin
 807   Find(ARowNum, Result);
 808   Insert(Result, TSRow.Create(ARowNum));
 809 end;
 810 
 811 {-- TSparseMatrix --}
 812 destructor TSparseMatrix.Destroy;
 813 begin
 814   FRows.Free;
 815   inherited Destroy;
 816 end;
 817 
 818 function  TSparseMatrix.GetCells(ARow, ACol: Integer): Extended;
 819 var
 820   Index: Integer;
 821   tRow: TSRow;
 822 begin
 823   { if index is invalid then raise exception }
 824   if (ACol < 0) or (ACol >= ColCount) or
 825      (ARow < 0) or (ARow >= RowCount)
 826   then  raise EMatrixException.Create('Index out of bounds');
 827   with FRows do
 828     if not Find(ARow, Index)
 829       then Result := 0
 830       else begin
 831         tRow := Items[Index];
 832         if not tRow.Find(ACol, Index)
 833           then Result := 0
 834           else Result := tRow.Items[Index]^.Value;
 835       end;
 836 end;
 837 
 838 procedure TSparseMatrix.SetCells(ARow, ACol: Integer; AValue: Extended);
 839 var
 840   Index: Integer;
 841   tRow: TSRow;
 842 begin
 843   { if index is invalid then raise exception }
 844   if (ACol < 0) or (ACol >= ColCount) or
 845      (ARow < 0) or (ARow >= RowCount)
 846   then  raise EMatrixException.Create('Index out of bounds');
 847   with FRows do
 848     if not Find(ARow, Index)
 849       then begin
 850         Index := Add(ARow);
 851         Items[Index].Add(ACol, AValue);
 852       end
 853       else begin
 854         tRow := Items[Index];
 855         if not tRow.Find(ACol, Index)
 856           then begin
 857             if AValue <> 0 then tRow.Add(ACol, AValue);
 858           end
 859           else begin
 860             if AValue <> 0
 861               then tRow.Items[Index]^.Value := AValue
 862               else begin
 863                 Dispose(PSCell(tRow.Items[Index]));
 864                 tRow.Delete(Index);
 865               end;
 866           end;
 867       end;
 868 end;
 869 
 870 procedure TSparseMatrix.SetRanges(NewRow, NewCol: Integer);
 871 var
 872   i, j, Index: Integer;
 873   tRow: TSRow;
 874 begin
 875   if (NewRow < 1) or (NewCol < 1)
 876     then raise Exception.Create('Invalid dimensions...');
 877   if FRows = nil then FRows := TRowList.Create;
 878   if RowCount > NewRow then
 879     with FRows do
 880       for i := Count - 1 downto 0 do
 881         if Items[i].RowNum >= NewRow then
 882         begin
 883           Items[i].Free;
 884           Delete(i);
 885         end
 886         else break;
 887   if ColCount > NewCol then
 888     with FRows do
 889       for i := 0 to Count - 1 do
 890       begin
 891         tRow := Items[i];
 892         for j := tRow.Count - 1 downto 0 do
 893           if tRow.Items[j]^.Col > NewCol then
 894           begin
 895             Dispose(tRow.Items[j]);
 896             tRow.Delete(j);
 897           end
 898           else break;
 899       end;
 900    FRowCount := NewRow;
 901    FColCount := NewCol;
 902 end;
 903 
 904 {-- TVector --}
 905 constructor TVector.Create(AOwner: TComponent);
 906 begin
 907   inherited Create(AOwner);
 908   SetCount(5);
 909 end;
 910 
 911 procedure TVector.Assign(Source: TVector);
 912 var
 913   i: Integer;
 914 begin
 915   Count := Source.Count;
 916   for i := 0 to Count - 1 do
 917     Cells[i] := Source[i];
 918 end;
 919 
 920 procedure TVector.AssignValue(Value: Extended);
 921 var
 922   i: Integer;
 923 begin
 924   for i := 0 to Count - 1 do
 925     Cells[i] := Value;
 926 end;
 927 
 928 procedure TVector.Add(AVector: TVector);
 929 var
 930   i: Integer;
 931 begin
 932   if Count <> AVector.Count
 933     then raise EMatrixException.Create('Can not add these Vectors');
 934   for i := 0 to Count - 1 do
 935     Cells[i] := Cells[i] + AVector[i];
 936 end;
 937 
 938 procedure TVector.MultValue(Value: Extended);
 939 var
 940   i: Integer;
 941 begin
 942   for i := 0 to Count - 1 do
 943     Cells[i] := Cells[i] * Value;
 944 end;
 945 
 946 procedure TVector.Transpose;
 947 var
 948   T: Extended;
 949   i: Integer;
 950 begin
 951   for i := 0 to (Count - 1) div 2 do
 952     begin
 953       T := Cells[i];
 954       Cells[i] := Cells[Count - i];
 955       Cells[Count - i] := T;
 956     end;
 957 end;
 958 
 959 function TVector.ScalarProduct(AVector: TVector): Extended;
 960 var
 961   i: Integer;
 962 begin
 963   if Count <> AVector.Count
 964     then raise EMatrixException.Create('Can not process these Vectors');
 965   Result := 0;
 966   for i := 0 to Count - 1 do
 967     Result := Cells[i] * AVector[i];
 968 end;
 969 
 970 procedure TVector.Exchange(i, j: Integer);
 971 var
 972   T: Extended;
 973 begin
 974   T := Cells[i];
 975   Cells[i] := Cells[j];
 976   Cells[j] := T;
 977 end;
 978 
 979 {-- TRealVector --}
 980 destructor  TRealVector.Destroy;
 981 begin
 982   FreeMem(FData, szExtended * Count);
 983   inherited Destroy;
 984 end;
 985 
 986 function  TRealVector.GetCells(APos: Integer): Extended;
 987 begin
 988   { if index is invalid then raise exception }
 989   if (APos < 0) or (APos > Count - 1)
 990     then  raise EMatrixException.Create('Index out of bounds');
 991   Result := FData^[APos];
 992 end;
 993 
 994 procedure TRealVector.SetCells(APos: Integer; AValue: Extended);
 995 begin
 996   { if index is invalid then raise exception }
 997   if (APos < 0) or (APos > Count - 1)
 998     then  raise EMatrixException.Create('Index out of bounds');
 999   FData^[APos] := AValue;
1000 end;
1001 
1002 procedure TRealVector.SetCount(NewCount: Integer);
1003 var
1004   OldCount: Integer;
1005 begin
1006   if (NewCount < 1)
1007     then raise Exception.Create('Invalid dimension...');
1008   if Count <> NewCount
1009   then begin
1010     OldCount := Count;
1011     { reallocate memory }
1012     {$IFDEF WIN32}
1013       ReAllocMem(FData, szExtended * NewCount);
1014     {$ELSE}
1015       FData := ReAllocMem(FData, szExtended * OldCount, szExtended * NewCount);
1016     {$ENDIF}
1017       if OldCount < NewCount
1018         then FillChar(FData^[OldCount], (NewCount - OldCount)*szExtended, 0);
1019     FCount := NewCount;
1020   end;
1021 end;
1022 
1023 {-- TSparseVector --}
1024 destructor  TSparseVector.Destroy;
1025 begin
1026   FData.Free;
1027   inherited Destroy;
1028 end;
1029 
1030 function  TSparseVector.GetCells(APos: Integer): Extended;
1031 var
1032   Index: Integer;
1033 begin
1034   { if index is invalid then raise exception }
1035   if (APos < 0) or (APos >= Count)
1036     then  raise EMatrixException.Create('Index out of bounds');
1037   if not FData.Find(APos, Index)
1038     then Result := 0
1039     else Result := FData.Items[Index]^.Value;
1040 end;
1041 
1042 procedure TSparseVector.SetCells(APos: Integer; AValue: Extended);
1043 var
1044   Index: Integer;
1045 begin
1046   { if index is invalid then raise exception }
1047   if (APos < 0) or (APos >= Count)
1048     then  raise EMatrixException.Create('Index out of bounds');
1049   if not FData.Find(APos, Index)
1050     then begin
1051       if AValue <> 0 then FData.Add(APos, AValue);
1052     end
1053     else begin
1054       if AValue <> 0
1055         then FData.Items[Index]^.Value := AValue
1056         else begin
1057           Dispose(PSCell(FData.Items[Index]));
1058           FData.Delete(Index);
1059         end;
1060     end;
1061 end;
1062 
1063 procedure TSparseVector.SetCount(NewCount: Integer);
1064 var
1065   i: Integer;
1066 begin
1067   if NewCount < 1
1068     then raise Exception.Create('Invalid dimension...');
1069   if FData = nil then FData := TSRow.Create(1);
1070   if Count > NewCount then
1071     for i := FData.Count - 1 downto 0 do
1072       if FData.Items[i]^.Col > NewCount then
1073       begin
1074         Dispose(FData.Items[i]);
1075         FData.Delete(i);
1076       end
1077       else break;
1078   FCount := NewCount;
1079 end;
1080 
1081 end.

View Code