Matrices

Definition

Matrices are instances of the class Matrix. Class Matrix is a template class: Matrix<T, Prop, Storage, Allocator>. As for vectors, T is the type of the elements to be stored (e.g. double).

Prop indicates that the matrix has given properties (symmetric, hermitian, positive definite or whatever). This template parameter is never used by Seldon; so the user may define its own properties. Thanks to this template parameter, the user may overload some functions based on the properties of the matrix. Seldon defines two properties: General (default) and Symmetric.

Storage defines how the matrix is stored. Matrices may be stored in several ways:

  • ColMajor: full matrix stored by columns (as in Fortran);
  • RowMajor: full matrix stored by rows (recommended in C++);
  • ColUpTriang: upper triangular matrix stored in a full matrix and by columns;
  • ColLoTriang: lower triangular matrix stored in a full matrix and by columns;
  • RowUpTriang: upper triangular matrix stored in a full matrix and by rows;
  • RowLoTriang: lower triangular matrix stored in a full matrix and by rows;
  • ColUpTriangPacked: upper triangular matrix stored in packed form (Blas format) and by columns;
  • ColLoTriangPacked: lower triangular matrix stored in packed form (Blas format) and by columns;
  • RowUpTriangPacked: upper triangular matrix stored in packed form (Blas format) and by rows;
  • RowLoTriangPacked: lower triangular matrix stored in packed form (Blas format) and by rows;
  • ColSym: symmetric matrix stored in a full matrix and by columns;
  • RowSym: symmetric matrix stored in a full matrix and by rows;
  • ColSymPacked: symmetric matrix stored in packed form (Blas format) and by columns;
  • RowSymPacked: symmetric matrix stored in packed form (Blas format) and by rows;
  • ColHerm: hermitian matrix stored in a full matrix and by columns;
  • RowHerm: hermitian matrix stored in a full matrix and by rows;
  • ColHermPacked: hermitian matrix stored in packed form (Blas format) and by columns;
  • RowHermPacked: hermitian matrix stored in packed form (Blas format) and by rows;
  • ColSparse: sparse matrix (Harwell-Boeing form) stored by columns;
  • RowSparse: sparse matrix (Harwell-Boeing form) stored by rows;
  • ColSymSparse: symmetric sparse matrix (Harwell-Boeing form, storing only the upper part) stored by columns;
  • RowSymSparse: symmetric sparse matrix (Harwell-Boeing form, storing only the upper part) stored by rows;
  • ColComplexSparse: complex sparse matrix (Harwell-Boeing form) stored by columns and for which the real part and the imaginary part are splitted into two sparse matrices;
  • RowComplexSparse: complex sparse matrix (Harwell-Boeing form) stored by rows and for which the real part and the imaginary part are splitted into two sparse matrices;
  • ColSymComplexSparse: symmetric sparse matrix (Harwell-Boeing form, storing only the upper part) stored by columns and for which the real part and the imaginary part are splitted into two sparse matrices;
  • RowSymComplexSparse: symmetric sparse matrix (Harwell-Boeing form, storing only the upper part) stored by rows and for which the real part and the imaginary part are splitted into two sparse matrices.
  • ArrayRowSparse: sparse matrix, each row being stored as a sparse vector.
  • ArrayColSparse: sparse matrix, each column being stored as a sparse vector.
  • ArrayRowSymSparse: symmetric sparse matrix, each row being stored as a sparse vector. Only upper part of the matrix is stored.
  • ArrayColSymSparse: symmetric sparse matrix, each column being stored as a sparse vector. Only upper part of the matrix is stored.
  • ArrayRowComplexSparse: sparse matrix, each row being stored as a sparse vector. Real part and imaginary part are stored separately.
  • ArrayRowSymComplexSparse: symmetric sparse matrix, each row being stored as a sparse vector. Real part and imaginary part are stored separately.

RowMajor is the default storage.

Finally, Allocator defines the way memory is managed. It is close to STL allocators. See the section "Allocators" for further details.

Declaration

There is a default Allocator (see the section "Allocators"), a default Storage (RowMajor) and a default property Prop (General). It means that the three last template parameters may be omitted. Then a matrix of integers may be declared thanks to the line: 

Matrix<int> M;

This defines a matrix of size 0 x 0, that is to say an empty matrix. To define a matrix of size 5 x 3, one may write: 

Matrix<int> M(5, 3);

Other declarations may be: 

Matrix<int, Symmetric> M(20, 30);
Matrix<int, General, ColMajor> M(20, 20);
Matrix<int, General, RowSparse, NewAlloc<int> > M;

Use of matrices

Access to elements is achieved through the operator(int, int), and indices start at 0: 

Matrix<double> M(5, 6);
M(0, 1) = -3.1;
M(0, 0) = 1.2 * M(0, 1);

To display matrices, there are two convenient options: 

M.Print();
cout << M << endl;

There are lots of methods that are described in the documentation. One may point out:

  • GetM() and GetN() return respectively the number of rows and the number of columns.
  • Fill fills with 1, 2, 3, etc. or fills the matrix with a given value.
  • Reallocate resizes the matrix (warning, data may be lost, depending on the allocator).
  • Resize resizes the matrix while keeping previous entries
  • Read, ReadText, Write, WriteText are useful methods for input/ouput operations.

Symmetric matrices in packed form are managed as square full matrices. There is a difference. If one affects the element (i, j) in the matrix, the element (j, i) will be affected at the same time, without any cost. A comprehensive test of full matrices is achieved in file test/program/matrices_test.cpp.

Sparse matrices - Harwell-Boeing form

Sparse matrices are not managed as full matrices. The access operator operator(int, int) is still available, but it doesn't allow any affectation, because it may return a (non-stored) zero of the matrix. In practice, one should deal with the underlying vectors that define elements and their indices: the vector of values data_, the vector of "start indices" ptr_ (i.e. indices - in the vector data_ - of the first element of each row or column) and the vector of row or column indices ind_.

Six types of storage are available : RowSparse, ColSparse, RowSymSparse, ColSymSparse, RowComplexSparse, RowSymComplexSparse. To deal with sparse matrices, the following methods should be useful:

  • SetData (or the constructor with the same arguments) to set the three vectors data_, ptr_ and ind_.
  • GetData, GetPtr and GetInd which return respectively data_, ptr_ and ind_ described above.
  • GetDataSize, GetPtrSize and GetIndSize which return respectively sizes of data_, ptr_ and ind_.

Sparse matrices - array of sparse vectors

Since the Harwell-Boeing form is difficult to handle, a more flexible form can be used in Seldon. Four types of storage are available : ArrayRowSparse, ArrayRowSymSparse, ArrayRowComplexSparse, ArrayRowSymComplexSparse. Their equivalents with a storage of columns : ArrayColSparse, ArrayColSymSparse, ArrayColComplexSparse, ArrayColSymComplexSparse are available as well, but often functions are implemented only for storage by rows. Therefore the user is strongly encourage to use only storages by rows. In this form, each row is stored as a sparse vector, allowing fast insertions of entries. Moreover, the access operator has the same functionnality as for dense matrices, because it allows affections.

The following methods should be useful:

  • Reallocate initialization of the matrix with a number of rows and columns. Previous entries are cleared.
  • Resize The number of rows and columns are modified, previous entries are kept.
  • GetDataSize returns number of stored values.
  • AddInteraction adds/inserts an entry in the matrix
  • AddInteractionRow adds/inserts severals entries in a row of the matrix
  • ReadText/WriteText reads/writes matrix in ascii.
The methods AddInteraction and AddInteractionRow inserts or adds entries in the right position, so that each row contains sorted elements (column numbers are sorted in ascending order). The function Copy converts matrices in any form, especially Harwell-Boeing form, since this last form induces faster matrix-vector product. 
int m = 10, n = 5;
Matrix<double, General, ArrayRowSparse> A(m, n);

// you can use AddInteraction
A.AddInteraction(3, 4, 3.0);
// or directly use the operator ()
// if the entry doesn't exist, it is added
A(2, 3) = 2.0;
// or add several entries at the same time
int nb_elt = 3;
IVect col_number(nb_elt);
DVect values(nb_elt);
col_number(0) = 1; col_number(0) = 3; col_number(0) = 2;
values(0) = 1.0; values(1) = -1.5; values(2) = 3.0;
A.AddInteractionRow(2, col_number, values);

// you can read/write the matrix
A.WriteText("Ah.dat"); A.ReadText("Ah.dat");
// and convert the matrix if you want a faster matrix-vector product
Matrix<double, General, RowSparse> Acsr;
Copy(A, Acsr);
A.Clear();
DVect x(n), b(m);
x.Fill();
// b = A*x
MltAdd(1.0, A, x, 0, b);
A comprehensive test of sparse matrices is achieved in file test/program/sparse_matrices_test.cpp.
Email Contact