List of available 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 rarely 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 three properties: General (default), Symmetric and Hermitian .

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

ColMajor: dense matrix stored by columns (as in Fortran);
RowMajor: dense matrix stored by rows (recommended in C++);
ColUpTriang: upper triangular matrix stored in a dense matrix and by columns;
ColLoTriang: lower triangular matrix stored in a dense matrix and by columns;
RowUpTriang: upper triangular matrix stored in a dense matrix and by rows;
RowLoTriang: lower triangular matrix stored in a dense 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 dense matrix and by columns;
RowSym: symmetric matrix stored in a dense 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 dense matrix and by columns;
RowHerm: hermitian matrix stored in a dense 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;
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.
BandedCol : banded matrix stored by columns as prescribed in Lapack documentation
ArrowCol : banded matrix stored by columns with additional rows and columns at the end

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. Additional storages dedicated to complex matrices are available when including SeldonComplexMatrix.hxx :

ColComplexSparse: complex sparse matrix (Harwell-Boeing form) stored by columns and for which the real part and the imaginary part are split 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 split 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 split 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 split into two sparse matrices.
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.

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;

Using 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);

However, the operator operator(int, int) can be used to modify the matrix only for storages RowMajor, ColMajor, RowSymPacked and ColSymPacked (dense symmetric or non-symmetric matrices). In order to propose a coherent and generic way to handle matrices, the four following access methods are available

operator (int i, int j) : returns always the value of A(i, j), by reference when it is possible.
method Val(int i, int j) : returns always a reference to A(i, j), as fast as possible without checking if the asked entry is on the upper part for symmetric matrices. An exception is raised when the entry A(i, j) does not correspond to a non-zero entry for sparse matrices
method Get(int i, int j) : returns always a reference to A(i, j), can modify the matrix for any (i, j) even for sparse matrices. In this last case, the matrix may be resized to insert a new non-zero entry.
method Set(int i, int j, T val) : sets A(i, j) to a given value, creates a new non-zero entry if needed for sparse matrices.

Here an exemple of use for hermitian matrices and sparse matrices

// hermitian matrix
Matrix<complex<double>, Hermitian, RowHermPacked> H(5, 5);
// use of Val to modify upper part of H
H.Val(0, 1) = complex<double>(-3.1, 0.9);
H.Val(0, 0) = 1.2;
H.Val(1, 2) = H(0,1)*4.1;

// operator () will always return the correct value of H(i, j)
cout << "H(2, 1) = " << H(2, 1) << endl;

// for Hermitian matrices, Get is equivalent to Val,
// and will fail when you will try to modify lower part
H.Get(2, 1) = complex<double>(0.5, 0.7); // => exception will be raised

// so if you want to set a value of the lower part, the only possible method is Set
H.Set(2, 1, complex<double>(0.5, 0.7)); // => CORRECT

// sparse matrix
Matrix<double, General, ArrayRowSparse> A(4, 4);
// use of Get(i, j) to modify A(i, j)
A.Get(2, 3) = 8.1;
A.Get(1, 0) = -1.5; 
A.Get(2, 3) += 0.5;
// you can use Set as well
A.Set(1, 4, 2.3);

// operator(int i, int j) can only be used in a read-only mode
cout << "A(3, 3) = " << A(3, 3) << endl;

// Val(int i, int j) to modify an already existing non-zero entry
A.Val(1, 4) = -2.7;

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/matrix_test.cpp (unitary tests are present in test/unit/matrix_test.cc)..

Sparse matrices

Sparse matrices - Harwell-Boeing form

Sparse matrices are not managed as dense 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_. The methods Get(i, j) or Set(i, j, val) can modify the matrix, but they are rather slow, since the arrays are resized for each new non-zero entry.

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 sometimes functions are implemented only for storage by rows. Therefore the user is strongly encourage to use only storages by rows. These storages are accessible if you have included SeldonSolver.hxx after the inclusion of Seldon.hxx (or if you have included SeldonLib.hxx):

// including Seldon and SeldonSolver
#include "Seldon.hxx"
#include "SeldonSolver.hxx"

using namespace Seldon;

// then you can use ArrayRowSparse storage
Matrix<double, General, ArrayRowSparse> A;

In these storages, each row is stored as a sparse vector, allowing fast insertions of entries. However, the access operator cannot be used as for dense matrices, since it doesn't allow modification, you should use function Get(i, j) instead.

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;
// format available only if you have included SeldonSolver.hxx
Matrix<double, General, ArrayRowSparse> A(m, n);

// you can use AddInteraction
A.AddInteraction(3, 4, 3.0);
// or directly use the method Get(i, j)
// if the entry doesn't exist, it is added
A.Get(2, 3) = 2.0;
// the following instruction also adds an entry :
double y = A(3, 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/matrix_sparse.cpp (unitary tests are present in test/unit/matrix_sparse_test.cc).