File numerics/src/tools/SparseBlockMatrix.h

Go to the source code of this file

Structure definition and functions related to SparseBlockStructuredMatrix

Defines

NUMERICS_SBM_FREE_BLOCK
NUMERICS_SBM_FREE_SBM

Functions

SparseBlockStructuredMatrix *SBM_new(void)

Creation of an empty Sparse Block Matrix.

Returns

a pointer on allocated and initialized space

void SBM_null(SparseBlockStructuredMatrix *sbm)

set Sparse Block Matrix.

fields to NULL

Parameters

sbm – a matrix

void SBM_gemv(unsigned int sizeX, unsigned int sizeY, double alpha, const SparseBlockStructuredMatrix *const A, const double *x, double beta, double *y)

SparseMatrix - vector product y = alpha*A*x + beta*y.

Parameters

  • sizeX – [in] dim of the vectors x

  • sizeY – [in] dim of the vectors y

  • alpha – [in] coefficient

  • A – [in] the matrix to be multiplied

  • x – [in] the vector to be multiplied

  • beta – [in] coefficient

  • y – [inout] the resulting vector

void SBM_gemv_3x3(unsigned int sizeX, unsigned int sizeY, const SparseBlockStructuredMatrix *const A, double *const x, double *y)

SparseMatrix - vector product y = A*x + y for block of size 3x3.

Parameters

  • sizeX – [in] dim of the vectors x

  • sizeY – [in] dim of the vectors y

  • A – [in] the matrix to be multiplied

  • x – [in] the vector to be multiplied

  • y – [inout] the resulting vector

void SBM_gemm_without_allocation(double alpha, const SparseBlockStructuredMatrix *const A, const SparseBlockStructuredMatrix *const B, double beta, SparseBlockStructuredMatrix *C)

SparseBlockStructuredMatrix - SparseBlockStructuredMatrix product C = alpha*A*B + beta*C The routine has to be used with precaution.

The allocation of C is not done since we want to add beta*C. We assume that the structure and the allocation of the matrix C are right. Especially:

  • the blocks C(i,j) must exists

  • the sizes of blocks must be consistent

  • no extra block must be present in C

Parameters

  • alpha – [in] coefficient

  • A – [in] the matrix to be multiplied

  • B – [in] the matrix to be multiplied

  • beta – [in] coefficient

  • C – [inout] the resulting matrix

SparseBlockStructuredMatrix *SBM_multiply(const SparseBlockStructuredMatrix *const A, const SparseBlockStructuredMatrix *const B)

SparseBlockStructuredMatrix - SparseBlockStructuredMatrix multiplication C = A *B Correct allocation is performed.

Parameters

  • A – [in] the matrix to be multiplied

  • B – [in] the matrix to be multiplied

Returns

C the resulting matrix

SparseBlockStructuredMatrix *SBM_zero_matrix_for_multiply(const SparseBlockStructuredMatrix *const A, const SparseBlockStructuredMatrix *const B)

Perform the allocation of a zero matrix that is compatible qith multiplication.

SparseBlockStructuredMatrix *SBM_add(SparseBlockStructuredMatrix *A, SparseBlockStructuredMatrix *B, double alpha, double beta)

SparseBlockStructuredMatrix - SparseBlockStructuredMatrix addition C = alpha*A + beta*B.

Parameters

  • A – [in] the matrix to be added

  • B – [in] the matrix to be added

  • alpha – [in] coefficient

  • beta – [in] coefficient

Returns

C the resulting matrix

void SBM_add_without_allocation(SparseBlockStructuredMatrix *A, SparseBlockStructuredMatrix *B, double alpha, double beta, SparseBlockStructuredMatrix *C, double gamma)

SparseBlockStructuredMatrix - SparseBlockStructuredMatrix addition C = alpha*A + beta*B + gamma*C without allocation.

We assume that C has the correct structure

Parameters

  • A – [in] the matrix to be added

  • B – [in] the matrix to be added

  • alpha – [in] coefficient

  • beta – [in] coefficient

  • gamma – [in] coefficient

  • C – [inout] the resulting matrix

void SBM_scal(double alpha, SparseBlockStructuredMatrix *A)

Multiply a matrix with a double alpha*A –> A.

Parameters

  • alpha – the coefficient

  • A – the matrix

void SBM_row_prod(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix *const A, const double *const x, double *y, int init)

Row of a SparseMatrix - vector product y = rowA*x or y += rowA*x, rowA being a row of blocks of A.

Parameters

  • sizeX – [in] dim of the vector x

  • sizeY – [in] dim of the vector y

  • currentRowNumber – [in] number of the required row of blocks

  • A – [in] the matrix to be multiplied

  • x – [in] the vector to be multiplied

  • y – [inout] the resulting vector

  • init – [in] = 0 for y += Ax, =1 for y = Ax

void SBM_row_prod_no_diag(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix *const A, const double *const x, double *y, int init)

Row of a SparseMatrix - vector product y = rowA*x or y += rowA*x, rowA being a row of blocks of A.

Parameters

  • sizeX – [in] dim of the vector x

  • sizeY – [in] dim of the vector y

  • currentRowNumber – [in] number of the required row of blocks

  • A – [in] the matrix to be multiplied

  • x – [in] the vector to be multiplied

  • y – [inout] the resulting vector

  • init – [in] = 0 for y += Ax, =1 for y = Ax

void SBM_row_prod_no_diag_3x3(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix *const A, double *const x, double *y)

Row of a SparseMatrix - vector product y = rowA*x or y += rowA*x, rowA being a row of blocks of A of size 3x3.

Parameters

  • sizeX – [in] dim of the vector x

  • sizeY – [in] dim of the vector y

  • currentRowNumber – [in] number of the required row of blocks

  • A – [in] the matrix to be multiplied

  • x – [in] the vector to be multiplied

  • y – [inout] the resulting vector

void SBM_row_prod_no_diag_2x2(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix *const A, double *const x, double *y)
void SBM_row_prod_no_diag_1x1(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix *const A, double *const x, double *y)
void SBM_extract_component_3x3(const SparseBlockStructuredMatrix *const A, SparseBlockStructuredMatrix *B, unsigned int *row_components, unsigned int row_components_size, unsigned int *col_components, unsigned int col_components_size)
void SBM_clear(SparseBlockStructuredMatrix *blmat)

Destructor for SparseBlockStructuredMatrix objects.

Parameters

blmat – SparseBlockStructuredMatrix the matrix to be destroyed.

void SBMfree(SparseBlockStructuredMatrix *A, unsigned int level)

To free a SBM matrix (for example allocated by NM_new_from_file).

Parameters

  • A – [in] the SparseBlockStructuredMatrix that mus be de-allocated.

  • level – [in] use NUMERICS_SBM_FREE_BLOCK | NUMERICS_SBM_FREE_SBM

void SBM_print(const SparseBlockStructuredMatrix *const m)

Screen display of the matrix content.

Parameters

m – the matrix to be displayed

void SBM_write_in_file(const SparseBlockStructuredMatrix *const m, FILE *file)

print in file of the matrix content

Parameters

  • m – the matrix to be displayed

  • file – the corresponding file

void SBM_read_in_file(SparseBlockStructuredMatrix *const M, FILE *file)

read in file of the matrix content without performing memory allocation

Parameters

  • M – the matrix to be displayed

  • file – the corresponding name of the file

SparseBlockStructuredMatrix *SBM_new_from_file(FILE *file)

Create from file a SparseBlockStructuredMatrix with memory allocation.

Parameters

file – the corresponding name of the file

Returns

the matrix to be displayed

void SBM_write_in_fileForScilab(const SparseBlockStructuredMatrix *const M, FILE *file)

print in file of the matrix content in Scilab format for each block

Parameters

  • M – the matrix to be displayed

  • file – the corresponding file

void SBM_write_in_filename(const SparseBlockStructuredMatrix *const M, const char *filename)

print in file of the matrix content

Parameters

  • M – the matrix to be displayed

  • filename – the corresponding file

void SBM_read_in_filename(SparseBlockStructuredMatrix *const M, const char *filename)

read in file of the matrix content

Parameters

  • M – the matrix to be displayed

  • filename – the corresponding name of the file

void SBM_clear_pred(SparseBlockStructuredMatrixPred *blmatpred)

Destructor for SparseBlockStructuredMatrixPred objects.

Parameters

blmatpred – SparseBlockStructuredMatrix, the matrix to be destroyed.

unsigned int *SBM_diagonal_block_indices(SparseBlockStructuredMatrix *const M)

Compute the indices of blocks of the diagonal block.

Parameters

M – the SparseBlockStructuredMatrix matrix

Returns

the indices for all the rows

unsigned int SBM_diagonal_block_index(SparseBlockStructuredMatrix *const M, unsigned int row)

Find index of the diagonal block in a row.

Parameters
Returns

pos the position of the block

int SBM_entry(SparseBlockStructuredMatrix *M, unsigned int row, unsigned int col, double val)

insert an entry into a SparseBlockStructuredMatrix.

This method is expensive in terms of memory management. For a lot of entries, use an alternative

Parameters
double SBM_get_value(const SparseBlockStructuredMatrix *const M, unsigned int row, unsigned int col)

get the element of row i and column j of the matrix M

Parameters
Returns

the value

int SBM_copy(const SparseBlockStructuredMatrix *const A, SparseBlockStructuredMatrix *B, unsigned int copyBlock)

Copy of a SBM A into B.

Parameters
Returns

0 if ok

int SBM_transpose(const SparseBlockStructuredMatrix *const A, SparseBlockStructuredMatrix *B)

Transpose by copy of a SBM A into B.

Parameters
Returns

0 if ok

int SBM_inverse_diagonal_block_matrix_in_place(const SparseBlockStructuredMatrix *M, int *ipiv)

Inverse (in place) a square diagonal block matrix.

Parameters
Returns

0 ik ok

void SBM_to_dense(const SparseBlockStructuredMatrix *const A, double *denseMat)

Copy a SBM into a Dense Matrix.

Parameters
int SBM_to_sparse(const SparseBlockStructuredMatrix *const A, CSparseMatrix *outSparseMat)

Copy a SBM into a Sparse (CSR) Matrix.

Parameters
Returns

0 if ok

int SBM_to_sparse_init_memory(const SparseBlockStructuredMatrix *const A, CSparseMatrix *sparseMat)

initMemory of a Sparse (CSR) Matrix form a SBM matrix

Parameters
Returns

0 if ok

void SBM_row_to_dense(const SparseBlockStructuredMatrix *const A, int row, double *denseMat, int rowPos, int rowNb)

Copy a block row of the SBM into a Dense Matrix.

Parameters

  • A – [in] the SparseBlockStructuredMatrix matrix to be inversed.

  • row – [in] the block row index copied.

  • denseMat – [in] pointer on the filled dense Matrix.

  • rowPos – [in] line pos in the dense matrix.

  • rowNb – [in] total number of line of the dense matrix. The number of line copied is contained in M.

void SBM_row_permutation(unsigned int *rowIndex, SparseBlockStructuredMatrix *A, SparseBlockStructuredMatrix *C)
Parameters

  • rowIndex – [in] permutation: the row numC of C is the row rowIndex[numC] of A.

  • A – [in] The source SBM.

  • C – [out] The target SBM. It assumes the structure SBM has been allocated. The memory allocation for its menber is done inside. NB : The blocks are not copied.

void SBM_column_permutation(unsigned int *colIndex, SparseBlockStructuredMatrix *A, SparseBlockStructuredMatrix *C)
Parameters

  • colIndex – [in] permutation: the col numC of C is the col colIndex[numC] of A.

  • A – [in] The source SBM.

  • C – [out] The target SBM. It assumes the structure SBM has been allocated. The memory allocation for its menber is done inside. NB : The blocks are not copied.

void SBCM_null(SparseBlockCoordinateMatrix *MC)
SparseBlockCoordinateMatrix *SBCM_new(void)
SparseBlockCoordinateMatrix *SBCM_new_3x3(unsigned int m, unsigned int n, unsigned int nbblocks, unsigned int *row, unsigned int *column, double *block)

allocate a SparseBlockCoordinateMatrix from a list of 3x3 blocks

Parameters

  • m – [in] the number of rows

  • n – [in] the number of colums

  • nbblocks – [in] the number of blocks

  • row – [in] a pointer to row of each block

  • column – [in] a pointer to column of each block

  • block – [in] a pointer to each block

Returns

a pointer to a SparseBlockCoordinateMatrix structure

void SBCM_free_3x3(SparseBlockCoordinateMatrix *MC)

free allocated memory in newSparseBlockCoordinateMatrix functions

Parameters

MC – [in] matrix pointer

SparseBlockStructuredMatrix *SBCM_to_SBM(SparseBlockCoordinateMatrix *MC)

copy a SparseBlockCoordinateMatrix to a SparseBlockStructuredMatrix

Parameters

MC – [in] the SparseBlockCoordinateMatrix matrix

Returns

a pointer to a SparseBlockCoordinateMatrix structure

void SBM_free_from_SBCM(SparseBlockStructuredMatrix *M)

free a SparseBlockStructuredMatrix created with SBCM_to_SBM

Parameters

M – [inout] a SparseBlockStructuredMatrix to free

int SBM_from_csparse(int blocksize, const CSparseMatrix *const sparseMat, SparseBlockStructuredMatrix *outSBM)

Copy a Sparse Matrix into a SBM, with fixed blocksize.

Parameters

  • blocksize – [in] the blocksize

  • sparseMat – [in] pointer on the Sparse Matrix

  • outSBM – [inout] pointer on an empty SparseBlockStructuredMatrix

Returns

0 in ok

struct SparseBlockStructuredMatrix
#include <>

Structure to store sparse block matrices with square diagonal blocks.

Related functions: SBM_gemv(), SBM_row_prod(), SBM_clear(), SBM_print, SBM_diagonal_block_index()

Note: the sparse format is the same as the one used by Boost C++ library to store compressed sparse row matrices. The same member names have been adopted in order to simplify usage from Siconos Kernel : filled1, filled2, index1_data, index2_data. Reference : http://ublas.sourceforge.net/refdoc/classboost_1_1numeric_1_1ublas_1_1compressed__matrix.html

If we consider the matrix M and the right-hand-side q defined as

\[\begin{split} M=\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] \quad, q=\left[\begin{array}{c}-1 \\ -1 \\ 0 \\ -1 \\ \hline 1 \\ 0 \\ \hline -1 \\ 2 \end{array}\right]. \end{split}\]

then

  • the number of non null blocks is 6 (nbblocks=6)

  • the number of rows of blocks is 3 (blocknumber0 =3) and the number of columns of blocks is 3 (blocknumber1 =3)

  • the vector blocksize0 is equal to {4,6,8} and the vector blocksize1 is equal to {4,6,8}

  • the integer filled1 is equal to 4

  • the integer filled2 is equal to 6

  • the vector index1_data is equal to {0,2,4,6}

  • the vector index2_data is equal to {0,1,1,2,0,2}

  • the block contains all non null block matrices stored in Fortran order (column by column) as block[0] = {1,2,0,5,2,1,0,0,0,0,1,-1,4,0,-1,6} block[1] = {3,4,0,0,-1,1,0,6} … block[5] = {2,-1,2,2}

Param nbblocks

the total number of non null blocks

Param **block

: *block contains the double values of one block in Fortran storage (column by column) **block is the list of non null blocks

Param blocknumber0

the first dimension of the block matrix (number of block rows)

Param blocknumber1

the second dimension of the block matrix (number of block columns)

Param *blocksize0

the list of sums of the number of rows of the first column of blocks of M: blocksize0[i] = blocksize0[i-1] + ni, ni being the number of rows of the block at row i

Param *blocksize1

the list of sums of the number of columns of the first row of blocks of M: blocksize1[i] = blocksize1[i-1] + ni, ni being the number of columns of the block at column i

Param filled1

index of the last non empty line + 1

Param filled2

number of non null blocks

Param index1_data

index1_data is of size equal to number of non empty lines + 1. A block with number blockNumber inside a row numbered rowNumber verify index1_data[rowNumber]<= blockNumber <index1_data[rowNumber+1]`

Param index2_data

index2_data is of size filled2 index2_data[blockNumber] -> columnNumber.

Public Members

unsigned int nbblocks
double **block
unsigned int blocknumber0
unsigned int blocknumber1
unsigned int *blocksize0
unsigned int *blocksize1
size_t filled1
size_t filled2
size_t *index1_data
size_t *index2_data
unsigned int *diagonal_blocks
NumericsDataVersion version

version of storage

struct SparseBlockCoordinateMatrix

Public Members

unsigned int nbblocks

number of blocks

unsigned int blocknumber0

number of rows

unsigned int blocknumber1

number of columns

double **block

block pointers

unsigned int *blocksize0

cumulative number of rows in blocks

unsigned int *blocksize1

cumulative number of columns in blocks

unsigned int *row

row indices

unsigned int *column

column indices

struct SparseBlockStructuredMatrixPred

Public Members

int nbbldiag
int **indic
int **indicop
double **submatlcp
double **submatlcpop
int **ipiv
int *sizesublcp
int *sizesublcpop
double **subq
double **bufz
double **newz
double **workspace
struct SBM_index_by_column

Public Members

size_t filled3
size_t filled4
size_t *index3_data
size_t *index4_data
size_t *blockMap