Linear Algebra in Siconos

An overview of matrices and vectors in Siconos : storage, usage, operations … Based on boost (https://www.boost.org/) ublas linear algebra library.

Defined in :

  • numerics/src/tools
  • kernel/src/utils/SiconosAlgebra

Vectors

Vectors in siconos are defined with SiconosVector object.

  • Contains real numbers (type : double precision)
  • Can be dense or 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