Linear Algebra in Siconos

Linear Algebra in Siconos#

An overview of matrices and vectors in Siconos : storage, usage, operations …

Defined in :

  • numerics/src/tools

  • kernel/src/utils/SiconosAlgebra

Siconos vectors and matrices in the kernel component are mostly an interface to boost ublas objects (https://www.boost.org/doc/libs/1_65_1/libs/numeric/ublas/doc/index.html)

Vectors#

Vectors in siconos are defined with SiconosVector object.

  • Contains real numbers (type : double precision)

  • Can be dense (type: Siconos::DENSE) or sparse (type Siconos::SPARSE).

Matrices#

Matrices are represented with SimpleMatrix class, which is a wrapper around boost ublas matrix types.

  • Contains real numbers (type : double precision)

  • Can be dense, triangular, symmetric, sparse, banded, zero or identity

Notes about SimpleMatrix#

  • Concerning sparse matrices, see http://freenet-homepage.de/guwi17/ublas/matrix_sparse_usage.html#Q2, for operations improvments.

  • Comparison between ublas (direct call) and Siconos are available in sandbox (developers only) gitlab project.

  • Different ways to compute matrix-vector or matrix-matrix products are proposed (prod, axpy_prod, gem…) based either on ublas or boost numeric bindings.

  • axpy_prod is only efficient for sparse or for large objects. For small matrices and vectors it is slower.

Matrix Storage in numerics component#

Numerics component proposes different ways to store ‘matrix-like’ objects, all handled through a C structure, NumericsMatrix .

Numerics component proposes different ways to store ‘matrix-like’ objects, all handled through a C structure, NumericsMatrix .

A number ( NumericsMatrix.storageType ) identify the type of storage while only one pointer among NumericsMatrix.matrixX, X = storageType = 0, 1 or 2, is not NULL and hold the values of the matrix.

At the time, the following storages are available:

  • “classical” (i.e. dense) column-major storage in a double*, NumericsMatrix.matrix0

  • sparse block storage, in a structure of type SparseBlockStructuredMatrix (warning: only for square matrices!!), NumericsMatrix.matrix1

  • sparse storage (csc, csr or triplet), based on CSparse (from T.Davis), NumericsMatrix.matrix2

As an example, consider the following matrix A of size 8X8:

\[\begin{split}\begin{equation*} A=\left[\begin{array}{cccc|cc|cc} 1 & 2 & 0 & 4 & 3 &-1 & 0 & 0\\ 2 & 1 & 0 & 0 & 4 & 1 & 0 & 0\\ 0 & 0 & 1 &-1 & 0 & 0 & 0 & 0\\ 5 & 0 &-1 & 6 & 0 & 6 & 0 & 0\\ \hline 0 & 0 & 0 & 0 & 1 & 0 & 0 & 5\\ 0 & 0 & 0 & 0 & 0 & 2 & 0 & 2\\ \hline 0 & 0 & 2 & 1 & 0 & 0 & 2 & 2\\ 0 & 0 & 2 & 2 & 0 & 0 & -1& 2\\ \end{array}\right] \end{equation*}\end{split}\]

How can we store this matrix ?

The first classical storage results in:

  • M.storageType = 0

  • M.size0 = 8, M.size1 = 8

  • M.matrix0 = [1 2 0 5 0 0 0 0 2 1 0 0 …]

    matrix0 being a double* of size 64.

For the second way of storage, SparseBlockStructuredMatrix we have:

  • M.storageType = 1

  • M.size0 = 8, M.size1 = 8

  • M.matrix1 a SparseBlockStructuredMatrix in which we save:

    • the number of non null blocks, 6 (matrix1->nbblocks) and the number of diagonal blocks, 3 (matrix1->size).

    • the vector matrix1->blocksize which collects the sum of diagonal blocks sizes until the present one, is equal to [4,6,8],

      blocksize[i] = blocksize[i-1] + ni, ni being the size of the diagonal block at row(block) i.

      Note that the last element of blocksize corresponds to the real size of the matrix.

    • the list of positions of non null blocks in vectors matrix1->ColumnIndex and matrix1->RowIndex, equal to [0,1,1,2,0,2] and [0,0,1,1,2,2]

    • the list of non null blocks, in matrix1->block, stored in Fortran order (column-major) as

      matrix1->block[0] = [1,2,0,5,2,1,0,0,0,0,1,-1,4,0,-1,6]

      matrix1->block[1] = [3,4,0,0,-1,1,0,6]

      …

      matrix1->block[5] = [2,-1,2,2]

Todo write proper doc for CSparse storage and complete the example above.

Functions on NumericsMatrix:

Create, fill and delete NumericsMatrix functions:

  • NM_create() : allocation without initial values

  • NM_create_from_data() : allocation and set default values from external data

  • NM_fill() : needs a pre-defined NumericsMatrix , set default values from external data

  • NM_free() : free a NumericsMatrix

These last two functions accept a data parameter, which if non-NULL contains the matrix data.

Linear Algebra:

The following linear algebra operation are supported:

  • BLAS-like functions:

    • product matrix - vector: NM_gemv() and NM_tgemv() (transpose)

    • product matrix - matrix: NM_gemm()

    • partial product matrix - vector: NM_row_prod()

-LAPACK-like functions -NM_gesv(): solve a linear system Ax = b

Input / Output:

  • NM_display() : display a NumericsMatrix

  • NM_display_row_by_row() : display a NumericsMatrix row by row

  • NM_write_in_filename() , NM_write_in_file() : save to filesystem

  • NM_read_in_filename() , NM_read_in_file() : fill a NumericsMatrix from a file

  • NM_new_from_file() : create new NumericsMatrix from a file

Contents