matrix.cpp File Reference

#include "geoslib_old_f.h"
#include "Basic/Utilities.hpp"
#include "Basic/Memory.hpp"
#include <math.h>
#include <cmath>

Functions
static double	_getTolInvert ()
static double	_getEpsMatrix ()
int	matrix_eigen (const double *a_in, int neq, double *value, double *vector)
void	matrix_product_safe (int n1, int n2, int n3, const double *v1, const double *v2, double *v3)
int	matrix_prod_norme (int transpose, int n1, int n2, const double *v1, const double *a, double *w)
void	matrix_transpose (int n1, int n2, VectorDouble &v1, VectorDouble &w1)
int	matrix_invert (double *a, int neq, int rank)
double	matrix_determinant (int neq, const VectorDouble &b)
int	matrix_cholesky_decompose (const double *a, double *tl, int neq)
void	matrix_cholesky_product (int mode, int neq, int nrhs, const double *tl, const double *a, double *x)
void	matrix_cholesky_invert (int neq, const double *tl, double *xl)
void	matrix_combine (int nval, double coeffa, const double *a, double coeffb, const double *b, double *c)
int	matrix_eigen_tridiagonal (const double *vecdiag, const double *vecinf, const double *vecsup, int neq, double *eigvec, double *eigval)

Function Documentation

_getEpsMatrix()

static double _getEpsMatrix ( )

static

_getTolInvert()

static double _getTolInvert ( )

static

matrix_cholesky_decompose()

int matrix_cholesky_decompose	(	const double *	a,
		double *	tl,
		int	neq
	)

Performs the Cholesky triangular decomposition of a definite positive symmetric matrix A = t(TL) * TL

Returns

Return code: >0 rank of zero pivot or 0 if no error

Parameters

[in]	a	symmetric matrix
[in]	neq	number of equations in the system
[out]	tl	Lower triangular matrix defined by column

Remarks

the matrix a[] is destroyed during the calculations

matrix_cholesky_invert()

void matrix_cholesky_invert	(	int	neq,
		const double *	tl,
		double *	xl
	)

Invert the Cholesky matrix

Parameters

[in]	neq	number of equations in the system
[in]	tl	lower triangular matrix defined by column
[out]	xl	lower triangular inverted matrix defined by column

matrix_cholesky_product()

void matrix_cholesky_product	(	int	mode,
		int	neq,
		int	nrhs,
		const double *	tl,
		const double *	a,
		double *	x
	)

Performs the product between a triangular and a square matrix TL is the lower triangular matrix and X is a square matrix

Parameters

[in]	mode	Type of calculations: 0 : X=TU%*A 1 : X=TL%*A 2 : X=A%*TU 3 : X=A%*TL 4 : X=t(A)%*TU 5 : X=t(A)%*TL
[in]	neq	number of equations in the system
[in]	nrhs	number of columns in x
[in]	tl	Triangular matrix defined by column (dimension: neq * neq)
[in]	a	matrix (dimension neq * nrhs)
[out]	x	resulting matrix (dimension neq * nrhs)

matrix_combine()

void matrix_combine	(	int	nval,
		double	coeffa,
		const double *	a,
		double	coeffb,
		const double *	b,
		double *	c
	)

Perform a linear combination of matrices or vectors [C] = 'coeffa' * [A] + 'coeffb' * [B]

Parameters

[in]	nval	Number of elements of the matrices or vectors
[in]	coeffa	Coefficient applied to the first matrix or vector
[in]	a	First matrix or vector (not used if NULL)
[in]	coeffb	Coefficient applied to the second matrix or vector
[in]	b	Second matrix or vector (not used if NULL)
[out]	c	Resulting matrix or vector

Remarks

No test is performed to check that the three matrices or vectors

have the same dimensions

Matrix c[] can coincide with matrices a[] or b[]

matrix_determinant()

double matrix_determinant	(	int	neq,
		const VectorDouble &	b
	)

Calculate the determinant of the square matrix (full storage)

Returns

Value of the determinant

Parameters

[in]	neq	Size of the matrix
[in]	b	Square matrix to be checked

matrix_eigen()

int matrix_eigen	(	const double *	a_in,
		int	neq,
		double *	value,
		double *	vector
	)

Calculates the eigen values and eigen vectors of a symmetric square matrix A X = X l

Returns

Return code:

0 no error

1 convergence problem

Parameters

[in]	a_in	square symmetric matrix (dimension = neq*neq)
[in]	neq	matrix dimension
[out]	value	matrix of the eigen values (dimension: neq)
[out]	vector	matrix of the eigen vectors (dimension: neq*neq)

Remarks

a_in is protected

matrix_eigen_tridiagonal()

int matrix_eigen_tridiagonal	(	const double *	vecdiag,
		const double *	vecinf,
		const double *	vecsup,
		int	neq,
		double *	eigvec,
		double *	eigval
	)

Transform a tridiagonal non-symmetric matrix into a symmetric one

Returns

Error return code (see remarks)

Parameters

[in]	vecdiag	Vector for the main diagonal
[in]	vecinf	Vector for the subdiagonal (in the last neq-1 positions)
[in]	vecsup	Vector for the superdiagonal (in the first neq positions)
[in]	neq	matrix dimension
[out]	eigvec	square symmetric matrix (dimension: neq * neq)
[out]	eigval	vector (dimensionL neq)

Remarks

Given the nonsymmetric tridiagonal matrix, we must have the products

of corresponding pairs of off-diagonal elements are all non-negative,

and zero only when both factors are zero

matrix_invert()

int matrix_invert	(	double *	a,
		int	neq,
		int	rank
	)

Invert a symmetric square matrix Pivots are assumed to be located on the diagonal

Returns

Return code: 0 no error; k if the k-th pivot is zero

Parameters

[in,out]	a	input matrix, destroyed in computation and replaced by resultant inverse
[in]	neq	number of equations in the matrix 'a'
[in]	rank	Type of message when inversion problem is encountered >=0: message involves 'rank+1' -1: neutral message -2: no message

Remarks

It is unnecessary to edit a message if inversion problem occurs

matrix_prod_norme()

int matrix_prod_norme	(	int	transpose,
		int	n1,
		int	n2,
		const double *	v1,
		const double *	a,
		double *	w
	)

Performs the product t(G) %*% A %*% G or G %*% A %*% t(G)

Returns

Error return code

Parameters

[in]	transpose	transposition mode -1 : transpose the first term +1 : transpose the last term
[in]	n1	matrix dimension
[in]	n2	matrix dimension
[in]	v1	rectangular matrix (n1,n2)
[in]	a	square matrix (optional)
[out]	w	square matrix

Remarks

According to the value of 'transpose':

-1: the output array has dimension (n2,n2)

+1: the output array has dimension (n1,n1)

According to the value of 'transpose':

-1: the optional array A has dimension (n1,n1)

+1: the optional array A has dimension (n2,n2)

matrix_product_safe()

void matrix_product_safe	(	int	n1,
		int	n2,
		int	n3,
		const double *	v1,
		const double *	v2,
		double *	v3
	)

Performs the product of two matrices

Parameters

[in]	n1	matrix dimension
[in]	n2	matrix dimension
[in]	n3	matrix dimension
[in]	v1	rectangular matrix (n1,n2)
[in]	v2	rectangular matrix (n2,n3)
[out]	v3	rectangular matrix (n1,n3)

Remarks

The matrix v3[] may NOT coincide with one of the two initial ones

matrix_transpose()

void matrix_transpose	(	int	n1,
		int	n2,
		VectorDouble &	v1,
		VectorDouble &	w1
	)

Transpose a (square or rectangular) matrix

Parameters

[in]	n1	matrix dimension
[in]	n2	matrix dimension
[in]	v1	rectangular matrix (n1,n2)
[out]	w1	rectangular matrix (n2,n1)

Remarks

The matrix w1[] may NOT coincide with v1[]