File numerics/src/tools/NumericsMatrix.h

Go to the source code of this file

Structure definition and functions related to matrix storage in Numerics.

Typedefs

typedef enum NumericsMatrix_types NM_types

Available types of storage for NumericsMatrix.

Enums

enum NM_gesv_opts

option for gesv factorization

Values:

NM_NONE

keep nothing

NM_KEEP_FACTORS

keep all the factorization data (useful to reuse the factorization)

NM_PRESERVE

keep the matrix as-is (useful for the dense case)

enum NumericsMatrix_types

Available types of storage for NumericsMatrix.

Values:

NM_DENSE

dense format

NM_SPARSE_BLOCK

sparse block format

NM_SPARSE

compressed column format

Functions

NumericsMatrix *NM_add(double alpha, NumericsMatrix *A, double beta, NumericsMatrix *B)

Add two matrices with coefficients C = alpha*A + beta*B.

Return
C a new NumericsMatrix
Parameters
  • alpha: the first coefficient
  • A: the first matrix
  • beta: the second coefficient
  • B: the second matrix

void NM_add_to_diag3(NumericsMatrix *M, double alpha)

Add a constant term to the diagonal elements, when the block of the SBM are 3x3.

Parameters
  • M: the matrix
  • alpha: the term to add

static void NM_assert(const int type, NumericsMatrix *M)

assert that a NumericsMatrix has the right structure given its type

Parameters
  • type: expected type
  • M: the matrix to check

int NM_check(const NumericsMatrix *const A)

Check the matrix (the sparse format for now)

Return
0 if the matrix storage is fine, 1 if not
Parameters
  • A: the matrix to check

void NM_clearCSC(NumericsMatrix *A)

Clear compressed column storage, if it is existent.

Parameters
  • A: a Numericsmatrix

void NM_clearCSCTranspose(NumericsMatrix *A)

Clear transposed compressed column storage, if it is existent.

Parameters
  • A: a Numericsmatrix

void NM_clearCSR(NumericsMatrix *A)

Clear compressed row storage, if it is existent.

Parameters
  • A: a Numericsmatrix

void NM_clearDense(NumericsMatrix *A)

Clear dense storage, if it is existent.

Parameters
  • A: a Numericsmatrix

void NM_clearSparse(NumericsMatrix *A)

Clear sparse data, if it is existent.

The linear solver parameters are also cleared.

Parameters
  • A: a Numericsmatrix

void NM_clearSparseBlock(NumericsMatrix *A)

Clear sparse block storage, if it is existent.

Parameters
  • A: a Numericsmatrix

void NM_clearSparseStorage(NumericsMatrix *A)

Clear triplet, csc, csc transposed storage, if they are existent.

Linear solver parameters are preserved.

Parameters
  • A: a Numericsmatrix

void NM_clearTriplet(NumericsMatrix *A)

Clear triplet storage, if it is existent.

Parameters
  • A: a Numericsmatrix

bool NM_compare(NumericsMatrix *A, NumericsMatrix *B, double tol)

compare to NumericsMatrix up to a given tolerance

Parameters

void NM_copy(const NumericsMatrix *const A, NumericsMatrix *B)

Copy a NumericsMatrix inside another NumericsMatrix (deep).

Reallocations are performed if B cannot hold a copy of A

Parameters

void NM_copy_diag_block3(NumericsMatrix *M, int block_row_nb, double **Block)

get a 3x3 diagonal block of a NumericsMatrix.

No allocation is done.

Parameters
  • M: a NumericsMatrix
  • block_row_nb: the number of the block row
  • Block: the target. In all cases (dense, sbm, and sparse) (*Block) must be allocated by caller. A copy is always performed

void NM_copy_to_sparse(const NumericsMatrix *const A, NumericsMatrix *B)

Copy a NumericsMatrix to s sparse one.

Allocation or reallocation are performed on B

Warning
It is assumed that B has been properly initialized: its storageType must be set to NM_SPARSE.
Parameters

NumericsMatrix *NM_create(int storageType, int size0, int size1)

create a NumericsMatrix and allocate the memory according to the matrix type

Return
a pointer to a NumericsMatrix
Parameters
  • storageType: the type of storage
  • size0: number of rows
  • size1: number of columns

NumericsMatrix *NM_create_from_data(int storageType, int size0, int size1, void *data)

create a NumericsMatrix and possibly set the data

Return
a pointer to a NumericsMatrix
Parameters
  • storageType: the type of storage
  • size0: number of rows
  • size1: number of columns
  • data: pointer to the matrix data. If NULL, all matrixX fields are set to NULL

NumericsMatrix *NM_create_from_file(FILE *file)
NumericsMatrix *NM_create_from_filename(char *filename)
CSparseMatrix *NM_csc(NumericsMatrix *A)

Creation, if needed, of compress column storage of a NumericsMatrix.

Return
the compressed column CSparseMatrix created in A.
Parameters

void NM_csc_alloc(NumericsMatrix *A, CS_INT nzmax)

Allocate a csc matrix in A.

Parameters
  • A: the matrix
  • nzmax: number of non-zero elements

void NM_csc_empty_alloc(NumericsMatrix *A, CS_INT nzmax)

Allocate a csc matrix in A and set the vector of column pointers to 0 such that the matrix is empty.

Parameters
  • A: the matrix
  • nzmax: number of non-zero elements

CSparseMatrix *NM_csc_trans(NumericsMatrix *A)

Creation, if needed, of the transposed compress column storage from compress column storage.

Return
the transposed compressed column matrix created in A.
Parameters

CSparseMatrix *NM_csr(NumericsMatrix *A)

Creation, if needed, of compress row storage of a NumericsMatrix.

Warning
This rely on the MKL
Return
the compressed row CSparseMatrix created in A.
Parameters

void NM_csr_alloc(NumericsMatrix *A, CS_INT nzmax)

Allocate a csr matrix in A.

Parameters
  • A: the matrix
  • nzmax: number of non-zero elements

void NM_dense_display(double *m, int nRow, int nCol, int lDim)

Screen display of the matrix content stored as a double * array in Fortran style.

Parameters
  • m: the matrix to be displayed
  • nRow: the number of rows
  • nCol: the number of columns
  • lDim: the leading dimension of M

void NM_dense_display_matlab(double *m, int nRow, int nCol, int lDim)

Screen display of the matrix content stored as a double * array in Fortran style.

Parameters
  • m: the matrix to be displayed
  • nRow: the number of rows
  • nCol: the number of columns
  • lDim: the leading dimension of M

void NM_dense_to_sparse(const NumericsMatrix *const A, NumericsMatrix *B)

matrix conversion display

matrix and vector display

void NM_display(const NumericsMatrix *const M)

Screen display of the matrix content.

Parameters
  • M: the matrix to be displayed

void NM_display_row_by_row(const NumericsMatrix *const m)

Screen display raw by raw of the matrix content.

Parameters
  • m: the matrix to be displayed

NumericsMatrix *NM_duplicate(NumericsMatrix *mat)

create a NumericsMatrix similar to the another one.

The structure is the same

Return
a pointer to a NumericsMatrix
Parameters
  • mat: the model matrix

double *NM_dWork(NumericsMatrix *A, int size)

Double workspace initialization, if needed.

Return
pointer on A->dWork allocated space of with the right size
Parameters

bool NM_equal(NumericsMatrix *A, NumericsMatrix *B)

compare to NumericsMatrix up to machine accuracy (DBL_EPSILON)

Parameters

void NM_extract_diag_block(NumericsMatrix *M, int block_row_nb, size_t start_row, int size, double **Block)

get the (square) diagonal block of a NumericsMatrix.

No allocation is done.

Parameters
  • M: a NumericsMatrix
  • block_row_nb: the number of the block Row. Useful only in sparse case
  • start_row: the starting row. Useful only in dense case.
  • size: of the diag block. Only useful in dense case.
  • Block: the target. In the dense and sparse case (*Block) must be allocated by caller. In case of SBM case **Bout contains the resulting block (from the SBM).

void NM_extract_diag_block3(NumericsMatrix *M, int block_row_nb, double **Block)

get a 3x3 diagonal block of a NumericsMatrix.

No allocation is done.

Parameters
  • M: a NumericsMatrix
  • block_row_nb: the number of the block row
  • Block: the target. In the dense and sparse case (*Block) must be allocated by caller. In case of SBM case **Bout contains the resulting block (from the SBM).

NumericsMatrix *NM_eye(int size)
void NM_fill(NumericsMatrix *M, int storageType, int size0, int size1, void *data)

fill an existing NumericsMatrix struct

Parameters
  • M: the struct to fill
  • storageType: the type of storage
  • size0: number of rows
  • size1: number of columns
  • data: pointer to the matrix data. If NULL, all matrixX fields are set to NULL

void NM_free(NumericsMatrix *m)

Free memory for a NumericsMatrix.

Warning: call this function only if you are sure that memory has been allocated for the structure in Numerics. This function is assumed that the memory is “owned” by this structure. Note that this function does not free m.

Parameters
  • m: the matrix to be deleted.

void NM_gemm(const double alpha, NumericsMatrix *A, NumericsMatrix *B, const double beta, NumericsMatrix *C)

Matrix matrix multiplication : C = alpha A B + beta C.

Parameters

void NM_gemv(const double alpha, NumericsMatrix *A, const double *x, const double beta, double *y)

Matrix vector multiplication : y = alpha A x + beta y.

Parameters
  • alpha: scalar
  • A: a NumericsMatrix
  • x: pointer on a dense vector of size A->size1
  • beta: scalar
  • y: pointer on a dense vector of size A->size1

static int NM_gesv(NumericsMatrix *A, double *b, bool preserve)

Direct computation of the solution of a real system of linear equations: A x = b.

Return
0 if successful, else the error is specific to the backend solver used
Parameters
  • A: a NumericsMatrix. On a dense factorisation A.iWork is initialized.
  • b: pointer on a dense vector of size A->size1
  • preserve: preserve the original matrix data. Only useful in the dense case, where the LU factorization is done in-place.

int NM_gesv_expert(NumericsMatrix *A, double *b, unsigned keep)

Direct computation of the solution of a real system of linear equations: A x = b.

The factorized matrix A is kept for future solve. If A is already factorized, the solve the linear system from it

Warning
this is not enable for all the solvers, your mileage may vary
Return
0 if successful, else the error is specific to the backend solver used
Parameters
  • A: a NumericsMatrix. On a dense factorisation A.iWork is initialized.
  • b: pointer on a dense vector of size A->size1
  • keep: if set to NM_KEEP_FACTORS, keep all the info related to the factorization to allow for future solves. If A is already factorized, just solve the linear system. If set to NM_PRESERVE, preserve the original matrix (just used in the dense case). if NM_NONE, discard everything.

int NM_gesv_expert_multiple_rhs(NumericsMatrix *A, double *b, unsigned int n_rhs, unsigned keep)
double NM_get_value(NumericsMatrix *M, int i, int j)

get the value of a NumericsMatrix.

Return
the value to be inserted.
Parameters

NumericsMatrixInternalData *NM_internalData(NumericsMatrix *A)

Get Matrix internal data with initialization if needed.

Return
a pointer on internal data.
Parameters

void NM_internalData_copy(const NumericsMatrix *const A, NumericsMatrix *B)

Copy the internalData structure.

Parameters
  • M: the matrix to modify

static void NM_internalData_new(NumericsMatrix *M)

Allocate the internalData structure (but not its content!)

Parameters
  • M: the matrix to modify

int NM_inv(NumericsMatrix *A, NumericsMatrix *Ainv)

Computation of the inverse of a NumericsMatrix A usinf NM_gesv_expert.

Parameters

int NM_inverse_diagonal_block_matrix_in_place(NumericsMatrix *A)
int NM_is_symmetric(NumericsMatrix *A)
void *NM_iWork(NumericsMatrix *A, size_t size, size_t sizeof_elt)

Integer work vector initialization, if needed.

Return
pointer on A->iWork allocated space of with the right size
Parameters
  • A: pointer on a NumericsMatrix.
  • size: number of element to allocate
  • sizeof_elt: of an element in bytes (result of sizeof for instance)

NumericsMatrix *NM_multiply(NumericsMatrix *A, NumericsMatrix *B)

Matrix matrix multiplication : C = A B.

Parameters

NumericsMatrix *NM_new(void)

Constructors and destructors.

Creation of an empty NumericsMatrix.

Return
a pointer to allocated space

NumericsMatrix *NM_new_from_file(FILE *file)

Create from file a NumericsMatrix with memory allocation.

Return
0 if the matrix
Parameters
  • file: the corresponding file

NumericsMatrix *NM_new_from_filename(char *filename)
NumericsMatrix *NM_new_SBM(int size0, int size1, SparseBlockStructuredMatrix *m1)

new NumericsMatrix with sparse storage from minimal set of data

Return
a pointer to a NumericsMatrix
Parameters

size_t NM_nnz(const NumericsMatrix *M)

return the number of non-zero element.

For a dense matrix, it is the product of the dimensions (e.g. an upper bound). For a sparse matrix, it is the true number

Return
the number (or an upper bound) of non-zero elements in the matrix
Parameters
  • M: the matrix

double NM_norm_1(NumericsMatrix *const A)

Compute the 1-norm of a sparse matrix = max (sum (abs (A))), largest column sum of a matrix (the sparse format for now)

Return
the norm
Parameters
  • A: the matrix

double NM_norm_inf(NumericsMatrix *const A)

Compute the inf-norm of a sparse matrix = max (sum (abs (A^T))), largest row sum of a matrix (the sparse format for now)

Return
the norm
Parameters
  • A: the matrix

static void NM_null(NumericsMatrix *A)

set NumericsMatrix fields to NULL

Parameters
  • A: a matrix

void NM_prod_mv_3x3(int sizeX, int sizeY, NumericsMatrix *A, double *const x, double *y)

Matrix - vector product.

Matrix - vector product y = A*x + y

Parameters
  • sizeX: dim of the vector x
  • sizeY: dim of the vector y
  • A: the matrix to be multiplied
  • x: the vector to be multiplied
  • y: the resulting vector

void NM_read_in_file(NumericsMatrix *const M, FILE *file)

Read in file of the matrix content without performing memory allocation.

Parameters
  • M: the matrix to be read
  • file: the corresponding file

void NM_read_in_file_scilab(NumericsMatrix *const M, FILE *file)

Read in file for scilab of the matrix content.

Parameters
  • M: the matrix to be read
  • file: the corresponding file

void NM_read_in_filename(NumericsMatrix *const M, const char *filename)

Read in file of the matrix content.

Parameters
  • M: the matrix to be read
  • filename: the corresponding name of the file

void NM_row_prod(int sizeX, int sizeY, int currentRowNumber, const NumericsMatrix *const A, const double *const x, double *y, int init)

Row of a Matrix - vector product y = rowA*x or y += rowA*x, rowA being a submatrix of A (sizeY rows and sizeX columns)

Parameters
  • sizeX: dim of the vector x
  • sizeY: dim of the vector y
  • currentRowNumber: position of the first row of rowA in A (warning: real row if A is a double*, block-row if A is a SparseBlockStructuredMatrix)
  • A: the matrix to be multiplied
  • x: the vector to be multiplied
  • y: the resulting vector
  • init: = 0 for y += Ax, =1 for y = Ax

void NM_row_prod_no_diag(size_t sizeX, size_t sizeY, int block_start, size_t row_start, NumericsMatrix *A, double *x, double *y, double *xsave, bool init)

Row of a Matrix - vector product y = rowA*x or y += rowA*x, rowA being a submatrix of A (sizeY rows and sizeX columns)

Parameters
  • sizeX: dim of the vector x
  • sizeY: dim of the vector y
  • block_start: block number (only used for SBM)
  • row_start: position of the first row of A (unused if A is SBM)
  • A: the matrix to be multiplied
  • x: the vector to be multiplied
  • y: the resulting vector
  • xsave: storage for saving the part of x set to 0
  • init: if True y = Ax, else y += Ax

void NM_row_prod_no_diag1x1(size_t sizeX, int block_start, size_t row_start, NumericsMatrix *A, double *x, double *y, bool init)
void NM_row_prod_no_diag3(size_t sizeX, int block_start, size_t row_start, NumericsMatrix *A, double *x, double *y, bool init)

Row of a Matrix - vector product y = rowA*x or y += rowA*x, rowA being a submatrix of A (3 rows and sizeX columns)

Parameters
  • sizeX: dim of the vector x
  • block_start: block number (only used for SBM)
  • row_start: position of the first row of A (unused if A is SBM)
  • A: the matrix to be multiplied
  • x: the vector to be multiplied
  • y: the resulting vector
  • init: if True y = Ax, else y += Ax

void NM_setSparseSolver(NumericsMatrix *A, unsigned solver_id)

Set the linear solver.

Parameters
  • A: the matrix
  • solver_id: the solver

double NM_symmetry_discrepancy(NumericsMatrix *A)
void NM_tgemv(const double alpha, NumericsMatrix *A, const double *x, const double beta, double *y)

Transposed matrix multiplication : y += alpha transpose(A) x + y.

Parameters
  • alpha: scalar
  • A: a NumericsMatrix
  • x: pointer on a dense vector of size A->size1
  • beta: scalar
  • y: pointer on a dense vector of size A->size1

int NM_to_dense(const NumericsMatrix *const A, NumericsMatrix *B)
NumericsMatrix *NM_transpose(NumericsMatrix *A)

new NumericsMatrix equal to the transpose of a given matrix

Return
a pointer to a NumericsMatrix
Parameters
  • A:

CSparseMatrix *NM_triplet(NumericsMatrix *A)

Creation, if needed, of triplet storage from sparse block storage.

Return
the triplet sparse Matrix created in A.
Parameters

void NM_triplet_alloc(NumericsMatrix *A, CS_INT nzmax)

Allocate a triplet matrix in A.

Parameters
  • A: the matrix
  • nzmax: maximum number of non-zero elements

void NM_update_size(NumericsMatrix *A)

update the size of the matrix based on the matrix data

Parameters
  • A: the matrix which size is updated

void NM_vector_display(double *m, int nRow)

Screen display of the vector content stored as a double * array.

Parameters
  • m: the vector to be displayed
  • nRow: the number of rows

void NM_write_in_file(const NumericsMatrix *const M, FILE *file)

PrintInFile of the matrix content.

Parameters
  • M: the matrix to be printed
  • file: filename the corresponding file

void NM_write_in_file_python(const NumericsMatrix *const M, FILE *file)

NM_write_in_file_python of the matrix content.

Parameters
  • M: the matrix to be printed
  • file: the corresponding file

void NM_write_in_file_scilab(const NumericsMatrix *const M, FILE *file)

NM_write_in_file_scilab of the matrix content.

Parameters
  • M: the matrix to be printed
  • file: the corresponding file

void NM_write_in_filename(const NumericsMatrix *const M, const char *filename)

matrix I/O

PrintInFile of the matrix content

Parameters
  • M: the matrix to be printed
  • filename: the corresponding name of the file

void NM_zentry(NumericsMatrix *M, int i, int j, double val)

setters and getters

insert an non zero entry into a NumericsMatrix. for storageType = NM_SPARSE, a conversion to triplet is done for performing the entry in the matrix. This method is expensice in terms of memory management. For a lot of entries, use preferably a triplet matrix.

Parameters
  • M: the NumericsMatrix
  • i: row index
  • j: column index
  • val: the value to be inserted.

NumericsSparseMatrix *numericsSparseMatrix(NumericsMatrix *A)

Creation, if needed, of sparse matrix storage.

Return
a pointer on the sparse matrix storage
Parameters

struct NumericsMatrix
#include <NumericsMatrix.h>

Interface to different type of matrices in numerics component.

See NM_* functions for linear algebra operations on dense, sparse block and sparse storage.

Public Members

NumericsMatrixInternalData *internalData

internal storage, used for workspace among other things

double *matrix0

dense storage

SparseBlockStructuredMatrix *matrix1

sparse block storage

NumericsSparseMatrix *matrix2

csc, csr or triplet storage

int size0

number of rows

int size1

number of columns

int storageType

the type of storage: 0: dense (double*), 1: SparseBlockStructuredMatrix, 2: classical sparse (csc, csr or triplet) via CSparse (from T.

Davis)

struct NumericsMatrixInternalData
#include <NumericsMatrix.h>

Structure for simple workspaces.

Public Members

double *dWork

double workspace

size_t dWorkSize

sizeof_elt of an element in bytes (result of sizeof for instance)

size of dWork

bool isInversed

true if the matrix containes its inverse (in place inversion)

bool isLUfactorized

true if the matrix has already been LU-factorized

void *iWork

integer workspace

size_t iWorkSize

size of iWork

size_t sizeof_elt