# 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

• 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