Classes of matrices and vectors in Delphi
This article discusses the design of types for working with linear algebra objects: vectors, matrices, quaternions. The classical application of the mechanism of overloading standard operations, the use of the “Copy On Write” technique and annotations is shown.
When working in the field of mathematical modeling, I often have to deal with computational algorithms that use operations on matrices, vectors, and quaternions. Surprisingly, I found that, despite the capabilities of the modern development environment, fellow programmers often use a procedural approach to solving such problems. So, to calculate the product of a matrix and a vector, types and functions like these are described:
It is proposed to use the object approach, which, in turn, implies the joint placement of data and methods of their processing. Let's see what possibilities Delphi provides for solving a similar class of tasks in object form. When designing object structures, we will proceed from the following requirements:
For convenience, we define auxiliary types:
Now we define the structures of the quaternion, vector and matrix:
We use the record , because operation overloading for a class construct in Delphi is not allowed. In addition, record objects have a useful property - their data is expanded in memory at the place of declaration, in other words, the record object is not a reference to an instance in dynamic memory.
However, in our case, the elements of vectors and matrices will be stored in a dynamic array, the object of which is a link. Therefore, it will be convenient to use explicit constructors. They initialize the internal fields, allocating memory for the required number of elements:
Quaternion at this stage does not require a constructor, because it stores data in a static array and is expanded in memory at the place of its declaration.
To access the elements here are indexer properties, it is convenient to make them default to omit the name. Access to the requested element occurs after checking its index for valid values. Shown implementation for TVector:
At this stage, to create our objects, you will have to use the following code:
Practice has shown that it is useful to have tools that allow you to use a more concise syntax to create a vector or matrix. To do this, add additional constructors, as well as implement the implicit cast operation, which allows you to overload ": =".
And implementation:
Now, to create and initialize a vector or matrix, it is enough to write:
Here, for example, only the overload of operation * will be implemented to multiply the matrix by a vector. Other operations can be found in the file attached to the article . A complete list of overload options is here .
The first argument to the Multiply () method is the matrix to the left of the * sign, the second argument is the column vector to the right of the * sign. The result of the work is a new vector whose object is created in the process of calculation. In the case of a discrepancy between the number of columns of the matrix and the number of vector elements, an exception is raised. Here's what the use of this operation in the program looks like:
It is convenient to use wrapper functions to construct anonymous vectors and matrices from literals of arrays "on the fly":
Using wrappers is as follows. The equivalent to the expression from the previous example is shown:
In addition to standard operations, it is useful to add specific methods to the types of our objects, for example, transposition or inversion. The following is an example of a matrix inversion (inversion) method. Despite its size, it is the fastest of all that I have seen (in high-level languages).
Using dynamic arrays to store elements of vectors and matrices leads to the fact that when you try to copy them in their entirety in the receiver object (the one to the left of ": ="), a copy of the link to this dynamic array is created.
For example, an attempt to save the value of the matrix M after evaluating the expression will invert the MStore matrix as well.
In order to correctly implement copying by value, we use the fact that, at a negative offset from the address of the first element of the dynamic array, along with the length value, a reference counter to this array is stored. If the counter value is 0, then the memory manager frees this array. If the counter value is 1, this means that there is only one reference to the array instance in memory.
Therefore, when copying, we must analyze the value of the counter and, if it is more than 1, then create a full copy of the array by copying it into the destination object element -by- element . Below is the code of a function that returns True only if the value of the reference count passed in the dynamic array input parameter exceeds 1.
At what point should you perform a full copy? This operation is quite expensive in time, so it makes no sense to perform it when accessing the element of the vector \ matrix for reading. If we have at least a thousand references to the original, if it is not subject to any changes, then all of them remain the same. Therefore, you need to copy only when accessing an item for recording. To do this, we modify the SetElement () methods for vectors and matrices, adding at the beginning a check for uniqueness of an instance of the FData array:
Thus, when you try to change the value of an element, a check is made for the uniqueness of the link, and if it is not confirmed, an element-by-element copy will be created in which the change will be made.
Elements of vectors and matrices should be accessed only after allocating memory for them. Element values are stored in a dynamic array, the dimensions of which are set in the constructor of the object. An implicit constructor call can occur when the object is initialized, or in the process of evaluating an expression.
Using implicit constructors can lead to errors when, sooner or later, access to an element of an un-created object is allowed. According to the rules of good form, the constructor should be called explicitly.
But what if the vectors and matrices in our program are hundreds and thousands? Consider a description of a class that uses many vectors and matrices as its fields.
It is required to develop a method for the automated initialization of all fields of the type TVector and TMatrix: allocate memory for vectors and matrices in accordance with the required number of elements and sizes. The annotation mechanism (or attributes, in terms of Delphi) will help us with this - a tool that allows types to be supplemented with arbitrary metadata. So, for each vector, the number of its elements must be known in advance, and for the matrix, the number of rows and columns.
Let's create a class that encapsulates the data on dimensions according to the rules for creating attribute classes.
The constructor gets the number of rows and columns, and in the case of a vector, we can do with only the number of rows. Now we supplement the definition of types from the previous listing with new annotations:
Below is the code that initializes objects of the type TVector and TMatrix based on information taken from the attributes.
The Init () procedure receives the input of the container object and its RTTI data. Next, a recursive traversal of all fields of the container occurs, and for all counter fields with type names “TVector” and “TMatrix” their constructors will be explicitly called.
We refine the TMovementScheme class using the Init () procedure:
Variant of calling Init () for an arbitrary record:
By default, Init () creates vectors with three elements, and 3x3 matrices, so in the declaration of TMovement and TMovementScheme types the attributes [TDim (3)] and [TDim (3, 3)] can be omitted, leaving only [TDim (4) ].
The file is attached to the article , in which the implementation of the described ideas is given in full.
When working in the field of mathematical modeling, I often have to deal with computational algorithms that use operations on matrices, vectors, and quaternions. Surprisingly, I found that, despite the capabilities of the modern development environment, fellow programmers often use a procedural approach to solving such problems. So, to calculate the product of a matrix and a vector, types and functions like these are described:
TVec3 = array[1..3] of Extended;
TMatrix3x3 = array[1..3, 1..3] of Extended;
function MVMult(M: TMatrix3x3; V: TVec3): TVec3;
It is proposed to use the object approach, which, in turn, implies the joint placement of data and methods of their processing. Let's see what possibilities Delphi provides for solving a similar class of tasks in object form. When designing object structures, we will proceed from the following requirements:
- Delphi XE7 Development Environment Version
- For numeric data, use the Extended type as the most accurate;
- Use dynamic arrays for data storage, as the sizes of vectors and matrices can be any and it is convenient to watch them in the debugger;
- Elements of vectors and matrices are numbered from 1, elements of quaternions from 0;
- For calculations, use the operations +, -, *, /;
- Provide the ability to transfer by value and copy vectors and matrices with the operation:: =;
- Provide the ability to automatically initialize vectors and matrices at predetermined sizes.
Design
For convenience, we define auxiliary types:
TAbstractVector = array of Extended;
TAbstractMatrix = array of array of Extended;
Now we define the structures of the quaternion, vector and matrix:
TQuaternion = record
private
FData: array[0..3] of Extended;
procedure SetElement(Index: Byte; Value: Extended);
function GetElement(Index: Byte): Extended;
public
property Element[Index: Byte]: Extended read GetElement write SetElement; default;
end;
TVector = record
private
FData: TAbstractVector;
FCount: Word;
procedure SetElement(Index: Word; Value: Extended);
function GetElement(Index: Word): Extended;
public
constructor Create(ElementsCount: Word);
property Count: Word read FCount;
property Elements[Index: Word]: Extended read GetElement write SetElement; default;
end;
TMatrix = record
private
FData: TAbstractMatrix;
FRowsCount: Word;
FColsCount: Word;
procedure SetElement(Row, Col: Word; Value: Extended);
function GetElement(Row, Col: Word): Extended;
public
constructor Create(RowsCount, ColsCount: Word);
property RowCount: Word read FRowsCount;
property ColCount: Word read FColsCount;
property Elements[Row, Col: Word]: Extended read GetElement write SetElement; default;
end;
We use the record , because operation overloading for a class construct in Delphi is not allowed. In addition, record objects have a useful property - their data is expanded in memory at the place of declaration, in other words, the record object is not a reference to an instance in dynamic memory.
However, in our case, the elements of vectors and matrices will be stored in a dynamic array, the object of which is a link. Therefore, it will be convenient to use explicit constructors. They initialize the internal fields, allocating memory for the required number of elements:
constructor TVector.Create(ElementsCount: Word);
begin
FCount := ElementsCount;
FData := nil;
SetLength(FData, FCount);
end;
constructor TMatrix.Create(RowsCount, ColsCount: Word);
begin
FRowsCount := RowsCount;
FColsCount := ColsCount;
FData := nil;
SetLength(FData, FRowsCount, FColsCount);
end;
Quaternion at this stage does not require a constructor, because it stores data in a static array and is expanded in memory at the place of its declaration.
To access the elements here are indexer properties, it is convenient to make them default to omit the name. Access to the requested element occurs after checking its index for valid values. Shown implementation for TVector:
function TVector.GetElement(Index: Word): Extended;
begin
{$R+}
Result := FData[Pred(Index)];
end;
procedure TVector.SetElement(Index: Word; Value: Extended);
begin
{$R+}
FData[Pred(Index)] := Value;
end;
At this stage, to create our objects, you will have to use the following code:
var
V: TVector;
. . .
V := TVector.Create(3);
V[1] := 1;
V[2] := 2;
V[3] := 3;
Practice has shown that it is useful to have tools that allow you to use a more concise syntax to create a vector or matrix. To do this, add additional constructors, as well as implement the implicit cast operation, which allows you to overload ": =".
TQuaternion = record
public
. . .
constructor Create(Q: TAbstractVector);
class operator Implicit(V: TAbstractVector): TQuaternion;
end;
TVector = record
public
. . .
constructor Create(V: TAbstractVector); overload;
class operator Implicit(V: TAbstractVector): TVector;
end;
TMatrix = record
public
. . .
constructor Create(M: TAbstractMatrix); overload;
class operator Implicit(M: TAbstractMatrix): TMatrix;
end;
And implementation:
constructor TQuaternion.Create(Q: TAbstractVector);
begin
if Length(Q) <> 4 then
raise EMathError.Create(WRONG_SIZE);
Move(Q[0], FData[0], SizeOf(FData));
end;
class operator TQuaternion.Implicit(V: TAbstractVector): TQuaternion;
begin
Result.Create(V);
end;
constructor TVector.Create(V: TAbstractVector);
begin
FCount := Length(V);
FData := Copy(V);
end;
class operator TVector.Implicit(V: TAbstractVector): TVector;
begin
Result.Create(V);
end;
constructor TMatrix.Create(M: TAbstractMatrix);
var
I: Integer;
begin
FRowsCount := Length(M);
FColsCount := Length(M[0]);
FData := nil;
SetLength(FData, FRowsCount, FColsCount);
for I := 0 to Pred(FRowsCount) do
FData[I] := Copy(M[I]);
end;
class operator TMatrix.Implicit(M: TAbstractMatrix): TMatrix;
begin
Result.Create(M);
end;
Now, to create and initialize a vector or matrix, it is enough to write:
var
V: TVector;
M: TMatrix;
. . .
V := [4, 5, 6];
//
M := [[1, 2, 3],
[4, 5, 6],
[7, 8, 9]];
Operations Overload
Here, for example, only the overload of operation * will be implemented to multiply the matrix by a vector. Other operations can be found in the file attached to the article . A complete list of overload options is here .
TMatrix = record
public
. . .
class operator Multiply(M: TMatrix; V: TVector): TVector;
end;
class operator TMatrix.Multiply(M: TMatrix; V: TVector): TVector;
var
I, J: Integer;
begin
if (M.FColsCount <> V.FCount) then
raise EMathError.Create(WRONG_SIZE);
Result.Create(M.FRowsCount);
for I := 0 to M.FRowsCount - 1 do
for J := 0 to M.FColsCount - 1 do
Result.FData[I] := Result.FData[I] + M.FData[I, J] * V.FData[J];
end;
The first argument to the Multiply () method is the matrix to the left of the * sign, the second argument is the column vector to the right of the * sign. The result of the work is a new vector whose object is created in the process of calculation. In the case of a discrepancy between the number of columns of the matrix and the number of vector elements, an exception is raised. Here's what the use of this operation in the program looks like:
var
V, VResult: TVector;
M: TMatrix;
. . .
VResult := M * V;
It is convenient to use wrapper functions to construct anonymous vectors and matrices from literals of arrays "on the fly":
function TVec(V: TAbstractVector): TVector;
begin
Result.Create(V);
end;
function TMat(M: TAbstractMatrix): TMatrix;
begin
Result.Create(M);
end;
function TQuat(Q: TAbstractVector): TQuaternion;
begin
Result.Create(Q);
end;
Using wrappers is as follows. The equivalent to the expression from the previous example is shown:
V := TMat([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]]) * TVec([4, 5, 6]);
In addition to standard operations, it is useful to add specific methods to the types of our objects, for example, transposition or inversion. The following is an example of a matrix inversion (inversion) method. Despite its size, it is the fastest of all that I have seen (in high-level languages).
TMatrix = record
public
. . .
function Inv: TMatrix;
end;
function TMatrix.Inv: TMatrix;
var
Ipiv, Indxr, Indxc: array of Integer;
DimMat, I, J, K, L, N, ICol, IRow: Integer;
Big, Dum, Pivinv: Extended;
begin
// Алгоритм Жордана.
if (FRowsCount <> FColsCount) then
raise EMathError.Create(NOT_QUAD);
Result := Self;
DimMat := FRowsCount;
SetLength(Ipiv, DimMat);
SetLength(Indxr, DimMat);
SetLength(Indxc, DimMat);
IRow := 1;
ICol := 1;
for I := 1 to DimMat do
begin
Big := 0;
for J := 1 to DimMat do
if (Ipiv[J - 1] <> 1) then
for K := 1 to DimMat do
if (Ipiv[K - 1] = 0) then
if (Abs(Result[J, K]) >= Big) then
begin
Big := Abs(Result[J, K]);
IRow := J;
ICol := K;
end;
Ipiv[ICol - 1] := Ipiv[ICol - 1] + 1;
if (IRow <> ICol) then
for L := 1 to DimMat do
begin
Dum := Result[IRow, L];
Result[IRow, L] := Result[ICol, L];
Result[ICol, L] := Dum;
end;
Indxr[I - 1] := IRow;
Indxc[I - 1] := ICol;
if Result[ICol, ICol] = 0 then
raise EMathError.Create(SINGULAR);
Pivinv := 1.0 / Result[ICol, ICol];
Result[ICol, ICol] := 1.0;
for L := 1 to DimMat do
Result[ICol, L] := Result[ICol, L] * Pivinv;
for N := 1 to DimMat do
if (N <> ICol) then
begin
Dum := Result[N, ICol];
Result[N, ICol] := 0.0;
for L := 1 to DimMat do
Result[N, L] := Result[N, L] - Result[ICol, L] * Dum;
end;
end;
for L := DimMat downto 1 do
if (Indxr[L - 1] <> Indxc[L - 1]) then
for K := 1 to DimMat do
begin
Dum := Result[K, Indxr[L - 1]];
Result[K, Indxr[L - 1]] := Result[K, Indxc[L - 1]];
Result[K, Indxc[L - 1]] := Dum;
end;
end;
Copy by value
Using dynamic arrays to store elements of vectors and matrices leads to the fact that when you try to copy them in their entirety in the receiver object (the one to the left of ": ="), a copy of the link to this dynamic array is created.
For example, an attempt to save the value of the matrix M after evaluating the expression will invert the MStore matrix as well.
var
M, MStore: TMatrix;
. . .
MStore := M;
M := M.Inv;
In order to correctly implement copying by value, we use the fact that, at a negative offset from the address of the first element of the dynamic array, along with the length value, a reference counter to this array is stored. If the counter value is 0, then the memory manager frees this array. If the counter value is 1, this means that there is only one reference to the array instance in memory.
Therefore, when copying, we must analyze the value of the counter and, if it is more than 1, then create a full copy of the array by copying it into the destination object element -by- element . Below is the code of a function that returns True only if the value of the reference count passed in the dynamic array input parameter exceeds 1.
{$POINTERMATH ON}
function NotUnique(var Arr): Boolean;
begin
Result := (PCardinal(Arr) - 2)^ > 1;
end;
At what point should you perform a full copy? This operation is quite expensive in time, so it makes no sense to perform it when accessing the element of the vector \ matrix for reading. If we have at least a thousand references to the original, if it is not subject to any changes, then all of them remain the same. Therefore, you need to copy only when accessing an item for recording. To do this, we modify the SetElement () methods for vectors and matrices, adding at the beginning a check for uniqueness of an instance of the FData array:
procedure TVector.SetElement(Index: Word; Value: Extended);
begin
{$R+}
CheckUnique;
FData[Pred(Index)] := Value;
end;
procedure TVector.CheckUnique;
begin
if NotUnique(FData) then
FData := Copy(FData);
end;
procedure TMatrix.SetElement(Row, Col: Word; Value: Extended);
begin
{$R+}
CheckUnique;
FData[Pred(Row), Pred(Col)] := Value;
end;
procedure TMatrix.CheckUnique;
var
I: Integer;
begin
if NotUnique(FData) then
begin
FData := Copy(FData);
for I := 0 to Pred(FRowsCount) do
FData[i] := Copy(FData[i]);
end;
end;
Thus, when you try to change the value of an element, a check is made for the uniqueness of the link, and if it is not confirmed, an element-by-element copy will be created in which the change will be made.
Annotations and automatic initialization
Elements of vectors and matrices should be accessed only after allocating memory for them. Element values are stored in a dynamic array, the dimensions of which are set in the constructor of the object. An implicit constructor call can occur when the object is initialized, or in the process of evaluating an expression.
var
V: TVector;
M: TMatrix;
begin
// V[1] := 1; // Ошибка: объект не создан
V := TVector.Create(4); // Явный вызов конструктора
M := TMatrix.Create(4, 4); // Явный вызов конструктора
// V := [1, 0, 0, 0]; // Неявный вызов конструктора
// V := M * TVec([1, 0, 0, 0]); // Неявный вызов конструктора
V[1] := 1; // Корректное обращение к элементу: объект создан
Using implicit constructors can lead to errors when, sooner or later, access to an element of an un-created object is allowed. According to the rules of good form, the constructor should be called explicitly.
But what if the vectors and matrices in our program are hundreds and thousands? Consider a description of a class that uses many vectors and matrices as its fields.
TMovement = record
R: TVector;
V: TVector;
W: TVector;
Color: TVector;
end;
TMovementScheme = class
private
FMovement: array[1..100] of TMovement;
FOrientation: TMatrix;
end;
It is required to develop a method for the automated initialization of all fields of the type TVector and TMatrix: allocate memory for vectors and matrices in accordance with the required number of elements and sizes. The annotation mechanism (or attributes, in terms of Delphi) will help us with this - a tool that allows types to be supplemented with arbitrary metadata. So, for each vector, the number of its elements must be known in advance, and for the matrix, the number of rows and columns.
Let's create a class that encapsulates the data on dimensions according to the rules for creating attribute classes.
TDim = class(TCustomAttribute)
private
FRowCount: Integer;
FColCount: Integer;
public
constructor Create(ARowCount: Integer; AColCount: Integer = 0); overload;
property RowCount: Integer read FRowCount;
property ColCount: Integer read FColCount;
end;
constructor TDim.Create(ARowCount: Integer; AColCount: Integer = 0);
begin
FRowCount := ARowCount;
FColCount := AColCount;
end;
The constructor gets the number of rows and columns, and in the case of a vector, we can do with only the number of rows. Now we supplement the definition of types from the previous listing with new annotations:
TMovement = record
[TDim(3)] R: TVector;
[TDim(3)] V: TVector;
[TDim(3)] W: TVector;
[TDim(4)] Color: TVector;
end;
TMovementScheme = class
private
FMovement: array[1..100] of TMovement;
[TDim(3, 3)] FOrientation: TMatrix;
end;
Below is the code that initializes objects of the type TVector and TMatrix based on information taken from the attributes.
procedure Init(Obj, TypeInfoOfObj: Pointer; Offset: Integer = 0);
const
DefaultRowCount = 3;
DefaultColCount = 3;
VectorTypeName = 'TVector';
MatrixTypeName = 'TMatrix';
var
RTTIContext: TRttiContext;
Field : TRttiField;
ArrFld: TRttiArrayType;
I: Integer;
Dim: TCustomAttribute;
RowCount, ColCount: Integer;
OffsetFromArray: Integer;
begin
for Field in RTTIContext.GetType(TypeInfoOfObj).GetFields do
begin
if Field.FieldType <> nil then
begin
RowCount := DefaultRowCount;
ColCount := DefaultColCount;
for Dim in Field.GetAttributes do
begin
RowCount := (Dim as TDim).RowCount;
ColCount := (Dim as TDim).ColCount;
end;
if Field.FieldType.TypeKind = tkArray then
begin
ArrFld := TRttiArrayType(Field.FieldType);
if ArrFld.ElementType.TypeKind = tkRecord then
begin
for I := 0 to ArrFld.TotalElementCount - 1 do
begin
OffsetFromArray := I * ArrFld.ElementType.TypeSize;
if ArrFld.ElementType.Name = VectorTypeName then
PVector(Integer(Obj) +
Field.Offset +
OffsetFromArray +
Offset)^ := TVector.Create(RowCount)
else if ArrFld.ElementType.Name = MatrixTypeName then
PMatrix(Integer(Obj) +
Field.Offset +
OffsetFromArray +
Offset)^ := TMatrix.Create(RowCount, ColCount)
else
Init(Obj, ArrFld.ElementType.Handle, Field.Offset + OffsetFromArray);
end;
end;
end
else if Field.FieldType.TypeKind = tkRecord then
begin
if Field.FieldType.Name = VectorTypeName then
PVector(Integer(Obj) +
Field.Offset +
Offset)^ := TVector.Create(RowCount)
else if Field.FieldType.Name = MatrixTypeName then
PMatrix(Integer(Obj) +
Field.Offset +
Offset)^ := TMatrix.Create(RowCount, ColCount)
else
Init(Obj, Field.FieldType.Handle, Field.Offset)
end;
end;
end;
end;
The Init () procedure receives the input of the container object and its RTTI data. Next, a recursive traversal of all fields of the container occurs, and for all counter fields with type names “TVector” and “TMatrix” their constructors will be explicitly called.
We refine the TMovementScheme class using the Init () procedure:
TMovementScheme = class
. . .
public
constructor Create;
end;
constructor TMovementScheme.Create;
begin
Init(Self, Self.ClassInfo);
end;
Variant of calling Init () for an arbitrary record:
var
Movement: TMovement;
. . .
Init(@Movement, TypeInfo(TMovement));
By default, Init () creates vectors with three elements, and 3x3 matrices, so in the declaration of TMovement and TMovementScheme types the attributes [TDim (3)] and [TDim (3, 3)] can be omitted, leaving only [TDim (4) ].
The file is attached to the article , in which the implementation of the described ideas is given in full.