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:

    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.

    Also popular now: