LinearAlgebra.md

Linear Algebra in ROOT

The linear algebra package is supposed to give a complete environment in ROOT to perform calculations like equation solving and eigenvalue decompositions. Most calculations are performed in double precision. For backward compatibility, some classes are also provided in single precision like TMatrixF, TMatrixFSym and TVectorF. Copy constructors exist to transform these into their double precision equivalent, thereby allowing easy access to decomposition and eigenvalue classes, only available in double precision.

The choice was made not to provide the less frequently used complex matrix classes. If necessary, users can always reformulate the calculation in 2 parts, a real one and an imaginary part. Although, a linear equation involving complex numbers will take about a factor of 8 more computations, the alternative of introducing a set of complex classes in this non-template library would create a major maintenance challenge.

Another choice was to fill in both the upper-right corner and the bottom-left corner of a symmetric matrix. Although most algorithms use only the upper-right corner, implementation of the different matrix views was more straightforward this way. When stored only the upper-right part is written to file.

For a detailed description of the interface, the user should look at the root reference guide at: http://root.cern.ch/root/Reference.html

Overview of Matrix Classes

The figure below shows an overview of the classes available in the linear algebra library, libMatrix.so. At the center is the base class TMatrixDBase from which three different matrix classes, TMatrixD, TMatrixDSym and TMatrixDFSparse derive. The user can define customized matrix operations through the classes TElementActionD and TElementsPosActionD.

Reference to different views of the matrix can be created through the classes on the right-hand side, see "Matrix Views". These references provide a natural connection to vectors.

Matrix decompositions (used in equation solving and matrix inversion) are available through the classes on the left-hand side (see "Matrix Decompositions"). They inherit from the TDecompBase class. The Eigen Analysis is performed through the classes at the top, see "Matrix Eigen Analysis". In both cases, only some matrix types can be analyzed. For instance, TDecompChol will only accept symmetric matrices as defined TMatrixDSym. The assignment operator behaves somewhat different than of most other classes. The following lines will result in an error:

   TMatrixD a(3,4);
   TMatrixD b(5,6);
   b = a;

It required to first resize matrix b to the shape of a.

   TMatrixD a(3,4);
   TMatrixD b(5,6);
   b.ResizeTo(a);
   b = a;

Matrix Properties

A matrix has five properties, which are all set in the constructor:

• precision - float or double. In the first case you will use the TMatrixF class family, in the latter case the TMatrixD one;

• type - general (TMatrixD), symmetric (TMatrixDSym) or sparse (TMatrixDSparse);

• size - number of rows and columns;

• index - range start of row and column index. By default these start at zero;

• sparse map - this property is only relevant for a sparse matrix. It indicates where elements are unequal zero.

Accessing Properties

The following table shows the methods to access the information about the relevant matrix property:

+-----------------------------+----------------------------------------------+ | Method | Descriptions | +-----------------------------+----------------------------------------------+ | Int_t GetRowLwb () | row lower-bound index | +-----------------------------+----------------------------------------------+ | Int_t GetRowUpb () | row upper-bound index | +-----------------------------+----------------------------------------------+ | Int_t GetNrows () | number of rows | +-----------------------------+----------------------------------------------+ | Int_t GetColLwb () | column lower-bound index | +-----------------------------+----------------------------------------------+ | Int_t GetColUpb () | column upper-bound index | +-----------------------------+----------------------------------------------+ | Int_t GetNcols () | number of columns | +-----------------------------+----------------------------------------------+ | Int_t GetNoElements () | number of elements, for a dense matrix this | | | equals: fNrows x fNcols | +-----------------------------+----------------------------------------------+ | Double_t GetTol () | tolerance number which is used in | | | decomposition operations | +-----------------------------+----------------------------------------------+ | Int_t *GetRowIndexArray | for sparse matrices, access to the row index | | () | of fNrows+1 entries | +-----------------------------+----------------------------------------------+ | Int_t *GetColIndexArray | for sparse matrices, access to the column | | () | index of fNelems entries | +-----------------------------+----------------------------------------------+

The last two methods in this table are specific to the sparse matrix, which is implemented according to the Harwell-Boeing format. Here, besides the usual shape/size descriptors of the matrix like fNrows, fRowLwb, fNcols and fColLwb, we also store a row index, fRowIndex and column index, fColIndex for the elements unequal zero: