MeteoIODoc 20241221.207bde49
mio::Matrix Class Reference

Detailed Description

This class implements the basic operations on matrices. Elements are access in matrix notation: that is A(1,2) represents the second element of the first line. Index go from 1 to nrows/ncols.

It might not be the best ever such implementation, but the goal is to provide a standalone matrix class. It might be later possible to chose between using the embedded implementation or to act as a front end to BLAS for those who have it installed on their system.

If the compilation flag NOSAFECHECKS is used, bounds check is turned off (leading to increased performances).

Author
Mathias Bavay

#include <Matrix.h>

Public Member Functions

 Matrix ()
 
 Matrix (const int &rows, const int &cols)
 A constructor that creates a matrix of a given size. More...
 
 Matrix (const size_t &rows, const size_t &cols)
 
 Matrix (const size_t &rows, const size_t &cols, const double &init)
 A constructor that creates a matrix filled with constant values. More...
 
 Matrix (const size_t &n, const double &init)
 A constructor that creates a diagonal matrix of size n. More...
 
 Matrix (const size_t &rows, const size_t &cols, const std::vector< double > data)
 A constructor that takes a data vector to fill the matrix. More...
 
void identity (const size_t &n, const double &init=1.)
 Convert the current matrix to a identity matrix of size n. More...
 
void resize (const size_t &rows, const size_t &cols)
 
void resize (const size_t &rows, const size_t &cols, const double &init)
 
void resize (const size_t &rows, const size_t &cols, const std::vector< double > &data)
 
void size (size_t &rows, size_t &cols) const
 get the dimensions of the current object More...
 
size_t getNx () const
 
size_t getNy () const
 
void clear ()
 free the memory and set the matrix dimensions to (0,0) More...
 
void random (const double &range)
 fill the matrix with random numbers. More...
 
double & operator() (const size_t &x, const size_t &y)
 
double operator() (const size_t &x, const size_t &y) const
 
double scalar () const
 Converts a 1x1 matrix to a scalar. More...
 
Matrix getT () const
 matrix transpose. More...
 
void T ()
 
Matrix getInv () const
 matrix invert. It first performs LU decomposition (Dolittle algorithm) and then computes the inverse by backward and forward solving of LU * A-1 = I. This inversion is in O(n³). see Press, William H.; Flannery, Brian P.; Teukolsky, Saul A.; Vetterling, William T. (1992), "LU Decomposition and Its Applications", Numerical Recipes in FORTRAN: The Art of Scientific Computing (2nd ed.), Cambridge University Press, pp. 34–42 More...
 
bool inv ()
 
Matrix getRow (const size_t i) const
 Copy a matrix row to a new 1-row-matrix. More...
 
void setRow (const size_t i, const Matrix &row)
 Set a matrix row from a 1-row-matrix. More...
 
Matrix getCol (const size_t j) const
 Copy a matrix column to a new 1-column-matrix. More...
 
void setCol (const size_t j, const Matrix &col)
 Set a matrix column from a 1-column-matrix. More...
 
Matrix getDiagonal () const
 Copy the matrix diagonal to a new 1-column-matrix. More...
 
Matrix extract (size_t r_low, size_t r_high, size_t c_low, size_t c_high) const
 Extract a submatrix to a new matrix. IOUtils::npos for an index means "from beginning" (i. e. 1) or "until end" (i. e. ncols resp. nrows). More...
 
double maxCoeff (size_t &max_row, size_t &max_col) const
 Find the first maximum value in the matrix. More...
 
double det () const
 matrix determinant More...
 
bool LU (Matrix &L, Matrix &U) const
 matrix LU decomposition. Perform LU decomposition by the Dolittle algorithm, (cf http://math.fullerton.edu/mathews/numerical/linear/dol/dol.html) More...
 
void partialPivoting (std::vector< size_t > &pivot_idx)
 matrix partial pivoting. This reorders the rows so that each diagonal element is the maximum in its column (see https://secure.wikimedia.org/wikipedia/en/wiki/Pivot_element) More...
 
void partialPivoting ()
 
void maximalPivoting ()
 
void swapCols (const size_t &j1, const size_t &j2)
 
const std::string toString (const int &precision=2, const bool &prettify=true) const
 
Matrixoperator+= (const Matrix &rhs)
 
const Matrix operator+ (const Matrix &rhs) const
 
Matrixoperator+= (const double &rhs)
 
const Matrix operator+ (const double &rhs) const
 
Matrixoperator-= (const Matrix &rhs)
 
const Matrix operator- (const Matrix &rhs) const
 
Matrixoperator-= (const double &rhs)
 
const Matrix operator- (const double &rhs) const
 
Matrixoperator*= (const Matrix &rhs)
 
const Matrix operator* (const Matrix &rhs) const
 
Matrixoperator*= (const double &rhs)
 
const Matrix operator* (const double &rhs) const
 
Matrixoperator/= (const double &rhs)
 
const Matrix operator/ (const double &rhs) const
 
bool operator== (const Matrix &) const
 Operator that tests for equality. More...
 
bool operator!= (const Matrix &) const
 Operator that tests for inequality. More...
 
bool isIdentity () const
 check if a matrix is the identity matrix More...
 

Static Public Member Functions

static double scalar (const Matrix &m)
 
static double dot (const Matrix &A, const Matrix &B)
 Dot product. More...
 
static Matrix T (const Matrix &m)
 
static Matrix solve (const Matrix &A, const Matrix &B)
 matrix solving for A·X=B. It first performs LU decomposition and then solves A·X=B by backward and forward solving of LU * X = B More...
 
static bool solve (const Matrix &A, const Matrix &B, Matrix &X)
 matrix solving for A·X=B. It first performs LU decomposition and then solves A·X=B by backward and forward solving of LU * X = B More...
 
static Matrix TDMA_solve (const Matrix &A, const Matrix &B)
 Solving system of equations using Thomas Algorithm The following function solves a A·X=B with X and B being vectors and A a tridiagonal matrix, using Thomas Algorithm (see http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm) More...
 
static bool TDMA_solve (const Matrix &A, const Matrix &B, Matrix &X)
 Solving system of equations using Thomas Algorithm The following function solves a A·X=B with X and B being vectors and A a tridiagonal matrix, using Thomas Algorithm (see http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm) More...
 
static void gaussElimination (Matrix &M, std::vector< size_t > &p)
 Gaussian elimimination with partial pivoting. This function performs row reduction with (virtual) row-swapping. Attention: The original matrix is destroyed! More...
 
static bool gaussSolve (Matrix &M, Matrix &A, Matrix &X)
 Equation system solver via Gaussian elimination. This function solves a set of equations M·X=A via Gauss elimination with partial pivoting. Attention: The original matrices are destroyed! More...
 
static bool gaussSolve (const Matrix &M, const Matrix &A, Matrix &X)
 Convenience call to the solver with input matrix preservation. This function makes a copy of the input matrices and invokes the solver on that, i. e. they are not changed in the process. More...
 
static bool gaussInverse (Matrix &M)
 Matrix inversion via Gaussian elimination. Attention: The original matrix is destroyed! More...
 
static bool gaussInverse (const Matrix &M, Matrix &Inv)
 Convenience call for matrix inversion with input matrix preservation. This function makes a copy of the input matrix, i. e. it is left unchanged by the solver. More...
 
static double gaussDet (Matrix &M)
 Get the matrix determinant via Gaussian elimination. Attention: The original matrix is destroyed! More...
 
static double normEuclid (const Matrix &vv)
 Calculate Euclidean norm (l2) for a vector. The matrix must either be a row or a column vector. More...
 
static unsigned int eigenvaluesJacobi (Matrix &AA, Matrix &DD)
 matrix bidiagonalization This uses Householder's reduction, see Golub, 1970. (see https://secure.wikimedia.org/wikipedia/en/wiki/Bidiagonalization) More...
 
static double jacobiEpsilon (Matrix &AA)
 
static void svdJacobi (const Matrix &MM, Matrix &UU, Matrix &SS, Matrix &VV)
 Singular Value Decomposition. This function calculates the SVD of a matrix. It uses an eigenvector search on M·M^T, which means problems may be squared. More...
 
static void sortEigenvalues (Matrix &EE, Matrix &VV)
 
static bool isIdentity (const Matrix &A)
 

Static Public Attributes

static const double epsilon = 1e-9
 
static const double epsilon_mtr = 1e-6
 

Protected Member Functions

size_t findMaxInCol (const size_t &col)
 
size_t findMaxInRow (const size_t &row)
 
void swapRows (const size_t &i1, const size_t &i2)
 

Protected Attributes

std::vector< double > vecData
 
size_t ncols
 
size_t nrows
 

Constructor & Destructor Documentation

◆ Matrix() [1/6]

mio::Matrix::Matrix ( )
inline

◆ Matrix() [2/6]

mio::Matrix::Matrix ( const int &  rows,
const int &  cols 
)

A constructor that creates a matrix of a given size.

Parameters
rowsnumber of rows of the new matrix
colsnumber of columns of the new matrix

◆ Matrix() [3/6]

mio::Matrix::Matrix ( const size_t &  rows,
const size_t &  cols 
)
inline

◆ Matrix() [4/6]

mio::Matrix::Matrix ( const size_t &  rows,
const size_t &  cols,
const double &  init 
)
inline

A constructor that creates a matrix filled with constant values.

Parameters
rowsnumber of rows of the new matrix
colsnumber of columns of the new matrix
initinitial value to fill the matrix with

◆ Matrix() [5/6]

mio::Matrix::Matrix ( const size_t &  n,
const double &  init 
)

A constructor that creates a diagonal matrix of size n.

Parameters
ndimension of the new square matrix
initinitial value to fill the matrix with

◆ Matrix() [6/6]

mio::Matrix::Matrix ( const size_t &  rows,
const size_t &  cols,
const std::vector< double >  data 
)
inline

A constructor that takes a data vector to fill the matrix.

Parameters
rowsnumber of rows of the new matrix
colsnumber of columns of the new matrix
datadata vector to init the matrix (must be of size rows*cols)

Member Function Documentation

◆ clear()

void mio::Matrix::clear ( )

free the memory and set the matrix dimensions to (0,0)

◆ det()

double mio::Matrix::det ( ) const

matrix determinant

Returns
determinant

◆ dot()

double mio::Matrix::dot ( const Matrix A,
const Matrix B 
)
static

Dot product.

Returns
scalar value

◆ eigenvaluesJacobi()

unsigned int mio::Matrix::eigenvaluesJacobi ( Matrix AA,
Matrix DD 
)
static

matrix bidiagonalization This uses Householder's reduction, see Golub, 1970. (see https://secure.wikimedia.org/wikipedia/en/wiki/Bidiagonalization)

Eigenvalue and eigenvector computation of symmetrical matrices. This function performs the iterative Jacobi method, cf. https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm Attention: The given matrix must be symmetrical on input!

Author
Michael Reisecker
Parameters
[in,out]AAThe matrix to search eigenvalues for. They will be found on its diagonal after transformation (unsorted).
[out]DDEigenvectors
Returns
Number of iterations that were needed

◆ extract()

Matrix mio::Matrix::extract ( size_t  r_low,
size_t  r_high,
size_t  c_low,
size_t  c_high 
) const

Extract a submatrix to a new matrix. IOUtils::npos for an index means "from beginning" (i. e. 1) or "until end" (i. e. ncols resp. nrows).

Parameters
r_lowStart from this row
r_highEnd at this row
c_lowStart from this column
c_highEnd at this column
Returns
The extracted submatrix

◆ findMaxInCol()

size_t mio::Matrix::findMaxInCol ( const size_t &  col)
protected

◆ findMaxInRow()

size_t mio::Matrix::findMaxInRow ( const size_t &  row)
protected

◆ gaussDet()

double mio::Matrix::gaussDet ( Matrix M)
static

Get the matrix determinant via Gaussian elimination. Attention: The original matrix is destroyed!

Author
Michael Reisecker
Parameters
[in]MMatrix to calculate the determinant for
Returns
The determinant

◆ gaussElimination()

void mio::Matrix::gaussElimination ( Matrix M,
std::vector< size_t > &  p 
)
static

Gaussian elimimination with partial pivoting. This function performs row reduction with (virtual) row-swapping. Attention: The original matrix is destroyed!

Author
Michael Reisecker
Parameters
[in,out]MThe matrix to reduce
[out]pUsed permutation vector (ignore p[0])

◆ gaussInverse() [1/2]

bool mio::Matrix::gaussInverse ( const Matrix M,
Matrix Inv 
)
static

Convenience call for matrix inversion with input matrix preservation. This function makes a copy of the input matrix, i. e. it is left unchanged by the solver.

Parameters
[in]MThe matrix to invert
[out]InvThe inverse matrix
Returns
False if the matrix is singular

◆ gaussInverse() [2/2]

bool mio::Matrix::gaussInverse ( Matrix M)
static

Matrix inversion via Gaussian elimination. Attention: The original matrix is destroyed!

Author
Michael Reisecker
Parameters
[in,out]MThe matrix to invert
Returns
False if the matrix is singular

◆ gaussSolve() [1/2]

bool mio::Matrix::gaussSolve ( const Matrix M,
const Matrix A,
Matrix X 
)
static

Convenience call to the solver with input matrix preservation. This function makes a copy of the input matrices and invokes the solver on that, i. e. they are not changed in the process.

Parameters
[in]MThe matrix M in M·X=A
[in]AThe solution vector (or matrix for multiple equations) A in M·X=A
[out]XThe calculated solution X in M·X=A
Returns
False if the equations are linearly dependent

◆ gaussSolve() [2/2]

bool mio::Matrix::gaussSolve ( Matrix M,
Matrix A,
Matrix X 
)
static

Equation system solver via Gaussian elimination. This function solves a set of equations M·X=A via Gauss elimination with partial pivoting. Attention: The original matrices are destroyed!

Author
Michael Reisecker
Parameters
[in]MThe matrix M in M·X=A
[in]AThe solution vector (or matrix for multiple equations) A in M·X=A
[out]XThe calculated solution X in M·X=A
Returns
False if the equations are linearly dependent

◆ getCol()

Matrix mio::Matrix::getCol ( const size_t  j) const

Copy a matrix column to a new 1-column-matrix.

Parameters
jColumn index
Returns
The column as a new matrix

◆ getDiagonal()

Matrix mio::Matrix::getDiagonal ( ) const

Copy the matrix diagonal to a new 1-column-matrix.

Returns
The diagonal as a new matrix

◆ getInv()

Matrix mio::Matrix::getInv ( ) const

matrix invert. It first performs LU decomposition (Dolittle algorithm) and then computes the inverse by backward and forward solving of LU * A-1 = I. This inversion is in O(n³). see Press, William H.; Flannery, Brian P.; Teukolsky, Saul A.; Vetterling, William T. (1992), "LU Decomposition and Its Applications", Numerical Recipes in FORTRAN: The Art of Scientific Computing (2nd ed.), Cambridge University Press, pp. 34–42

matrix inversion performance with the matrix dimension (running on a single core of a Intel i7-9750H)
Returns
inversed matrix

◆ getNx()

size_t mio::Matrix::getNx ( ) const
inline

◆ getNy()

size_t mio::Matrix::getNy ( ) const
inline

◆ getRow()

Matrix mio::Matrix::getRow ( const size_t  i) const

Copy a matrix row to a new 1-row-matrix.

Parameters
iRow index
Returns
The row as a new matrix

◆ getT()

Matrix mio::Matrix::getT ( ) const

matrix transpose.

Returns
transposed matrix

◆ identity()

void mio::Matrix::identity ( const size_t &  n,
const double &  init = 1. 
)

Convert the current matrix to a identity matrix of size n.

Parameters
ndimension of the new square matrix
initinitial value to fill the matrix diagonal with

◆ inv()

bool mio::Matrix::inv ( )

◆ isIdentity() [1/2]

bool mio::Matrix::isIdentity ( ) const

check if a matrix is the identity matrix

Returns
true if it is I

◆ isIdentity() [2/2]

bool mio::Matrix::isIdentity ( const Matrix A)
static

◆ jacobiEpsilon()

double mio::Matrix::jacobiEpsilon ( Matrix AA)
static

◆ LU()

bool mio::Matrix::LU ( Matrix L,
Matrix U 
) const

matrix LU decomposition. Perform LU decomposition by the Dolittle algorithm, (cf http://math.fullerton.edu/mathews/numerical/linear/dol/dol.html)

Parameters
Llower diagonal matrix
Uupper diagonal matrix
Returns
false if the decomposition can not be performed (division by zero)

◆ maxCoeff()

double mio::Matrix::maxCoeff ( size_t &  max_row,
size_t &  max_col 
) const

Find the first maximum value in the matrix.

Parameters
max_rowRow index of the maximum
max_colColumn index of the maximum
Returns
The maximum value

◆ maximalPivoting()

void mio::Matrix::maximalPivoting ( )

◆ normEuclid()

double mio::Matrix::normEuclid ( const Matrix vv)
static

Calculate Euclidean norm (l2) for a vector. The matrix must either be a row or a column vector.

Parameters
[in]vvVector to calculate norm for
Returns
L2 norm

◆ operator!=()

bool mio::Matrix::operator!= ( const Matrix in) const

Operator that tests for inequality.

◆ operator()() [1/2]

double & mio::Matrix::operator() ( const size_t &  x,
const size_t &  y 
)

◆ operator()() [2/2]

double mio::Matrix::operator() ( const size_t &  x,
const size_t &  y 
) const

◆ operator*() [1/2]

const Matrix mio::Matrix::operator* ( const double &  rhs) const

◆ operator*() [2/2]

const Matrix mio::Matrix::operator* ( const Matrix rhs) const

◆ operator*=() [1/2]

Matrix & mio::Matrix::operator*= ( const double &  rhs)

◆ operator*=() [2/2]

Matrix & mio::Matrix::operator*= ( const Matrix rhs)

◆ operator+() [1/2]

const Matrix mio::Matrix::operator+ ( const double &  rhs) const

◆ operator+() [2/2]

const Matrix mio::Matrix::operator+ ( const Matrix rhs) const

◆ operator+=() [1/2]

Matrix & mio::Matrix::operator+= ( const double &  rhs)

◆ operator+=() [2/2]

Matrix & mio::Matrix::operator+= ( const Matrix rhs)

◆ operator-() [1/2]

const Matrix mio::Matrix::operator- ( const double &  rhs) const

◆ operator-() [2/2]

const Matrix mio::Matrix::operator- ( const Matrix rhs) const

◆ operator-=() [1/2]

Matrix & mio::Matrix::operator-= ( const double &  rhs)

◆ operator-=() [2/2]

Matrix & mio::Matrix::operator-= ( const Matrix rhs)

◆ operator/()

const Matrix mio::Matrix::operator/ ( const double &  rhs) const

◆ operator/=()

Matrix & mio::Matrix::operator/= ( const double &  rhs)

◆ operator==()

bool mio::Matrix::operator== ( const Matrix in) const

Operator that tests for equality.

◆ partialPivoting() [1/2]

void mio::Matrix::partialPivoting ( )

◆ partialPivoting() [2/2]

void mio::Matrix::partialPivoting ( std::vector< size_t > &  pivot_idx)

matrix partial pivoting. This reorders the rows so that each diagonal element is the maximum in its column (see https://secure.wikimedia.org/wikipedia/en/wiki/Pivot_element)

Parameters
pivot_idxnew indices (to apply when solving A * X = B, for example

◆ random()

void mio::Matrix::random ( const double &  range)

fill the matrix with random numbers.

Parameters
rangerange of the randoms numbers (they will be between -range and +range)

◆ resize() [1/3]

void mio::Matrix::resize ( const size_t &  rows,
const size_t &  cols 
)

◆ resize() [2/3]

void mio::Matrix::resize ( const size_t &  rows,
const size_t &  cols,
const double &  init 
)

◆ resize() [3/3]

void mio::Matrix::resize ( const size_t &  rows,
const size_t &  cols,
const std::vector< double > &  data 
)

◆ scalar() [1/2]

double mio::Matrix::scalar ( ) const

Converts a 1x1 matrix to a scalar.

Returns
scalar value

◆ scalar() [2/2]

double mio::Matrix::scalar ( const Matrix m)
static

◆ setCol()

void mio::Matrix::setCol ( const size_t  j,
const Matrix col 
)

Set a matrix column from a 1-column-matrix.

Parameters
jColumn index
colThe new column

◆ setRow()

void mio::Matrix::setRow ( const size_t  i,
const Matrix row 
)

Set a matrix row from a 1-row-matrix.

Parameters
iRow index
rowThe new row

◆ size()

void mio::Matrix::size ( size_t &  rows,
size_t &  cols 
) const

get the dimensions of the current object

Parameters
rowsnumber of rows of the matrix
colsnumber of columns of the matrix

◆ solve() [1/2]

Matrix mio::Matrix::solve ( const Matrix A,
const Matrix B 
)
static

matrix solving for A·X=B. It first performs LU decomposition and then solves A·X=B by backward and forward solving of LU * X = B

Parameters
AA matrix
BB matrix
Returns
solution matrix

◆ solve() [2/2]

bool mio::Matrix::solve ( const Matrix A,
const Matrix B,
Matrix X 
)
static

matrix solving for A·X=B. It first performs LU decomposition and then solves A·X=B by backward and forward solving of LU * X = B

Parameters
AA matrix
BB matrix
Xsolution matrix
Returns
true is success

◆ sortEigenvalues()

void mio::Matrix::sortEigenvalues ( Matrix EE,
Matrix VV 
)
static

◆ svdJacobi()

void mio::Matrix::svdJacobi ( const Matrix MM,
Matrix UU,
Matrix SS,
Matrix VV 
)
static

Singular Value Decomposition. This function calculates the SVD of a matrix. It uses an eigenvector search on M·M^T, which means problems may be squared.

Author
Michael Reisecker
Parameters
[in]MMThe matrix to factorize.
[out]UULeft singular vectors
[out]SSDiagonal matrix with the singular values
[out]VVRight singular vectors

◆ swapCols()

void mio::Matrix::swapCols ( const size_t &  j1,
const size_t &  j2 
)

◆ swapRows()

void mio::Matrix::swapRows ( const size_t &  i1,
const size_t &  i2 
)
protected

◆ T() [1/2]

void mio::Matrix::T ( )

◆ T() [2/2]

Matrix mio::Matrix::T ( const Matrix m)
static

◆ TDMA_solve() [1/2]

Matrix mio::Matrix::TDMA_solve ( const Matrix A,
const Matrix B 
)
static

Solving system of equations using Thomas Algorithm The following function solves a A·X=B with X and B being vectors and A a tridiagonal matrix, using Thomas Algorithm (see http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)

Author
Nander Wever, Mathias Bavay
Parameters
AA matrix
BB matrix
Returns
solution matrix

◆ TDMA_solve() [2/2]

bool mio::Matrix::TDMA_solve ( const Matrix A,
const Matrix B,
Matrix X 
)
static

Solving system of equations using Thomas Algorithm The following function solves a A·X=B with X and B being vectors and A a tridiagonal matrix, using Thomas Algorithm (see http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)

Author
Nander Wever, Mathias Bavay
Parameters
AA matrix
BB matrix
Xsolution matrix
Returns
true is success

◆ toString()

const std::string mio::Matrix::toString ( const int &  precision = 2,
const bool &  prettify = true 
) const

Member Data Documentation

◆ epsilon

const double mio::Matrix::epsilon = 1e-9
static

◆ epsilon_mtr

const double mio::Matrix::epsilon_mtr = 1e-6
static

◆ ncols

size_t mio::Matrix::ncols
protected

◆ nrows

size_t mio::Matrix::nrows
protected

◆ vecData

std::vector<double> mio::Matrix::vecData
protected

The documentation for this class was generated from the following files: