#include "geoslib_old_f.h"
#include "geoslib_enum.h"
#include "Enum/ECst.hpp"
#include "Basic/Utilities.hpp"
#include "Basic/VectorHelper.hpp"
#include "Basic/OptCst.hpp"
#include <math.h>
#include <cmath>

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

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

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, double v1, double w1)

void	matrix_int_transpose (int n1, int n2, int v1, int w1)

void	matrix_transpose_in_place (int n1, int n2, double *v1)

void	matrix_int_transpose_in_place (int n1, int n2, int *v1)

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

int	matrix_invert_copy (const double a, int neq, double b)

int	matrix_solve (int mode, const double a, const double b, double x, int neq, int nrhs, int pivot)

int	is_matrix_symmetric (int neq, const double *a, int verbose)

int	is_matrix_definite_positive (int neq, const double a, double valpro, double *vecpro, int verbose)

int	is_matrix_non_negative (int nrow, int ncol, double *a, int verbose)

double	matrix_determinant (int neq, const double *b)

int	is_matrix_rotation (int neq, const double *a, int verbose)

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_cholesky_norme (int mode, int neq, const double tl, const double a, double *b)

int	matrix_invgen (double a, int neq, double tabout, double *cond)

double	matrix_normA (double b, double a, int neq, int subneq)

void	matrix_triangle_to_square (int mode, int neq, const double tl, double a)

void	matrix_produit_cholesky (int neq, const double tl, double a)

VectorDouble	matrix_produit_cholesky_VD (int neq, const double *tl)

void	matrix_svd_inverse (int neq, double s, double u, double v, double tabout)

int	matrix_invsvdsym (double *mat, int neq, int rank)

void	matrix_manage (int nrows, int ncols, int nr, int nc, int rowsel, int colsel, double v1, double v2)

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

double	matrix_norminf (int nval, double *tab)

void	matrix_product_by_diag (int mode, int neq, double a, double c, double *b)

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

int	matrix_qo (int neq, double hmat, double gmat, double *xmat)

int	matrix_qoc (int flag_invert, int neq, double hmat, double gmat, int na, double amat, double bmat, double xmat, double lambda)

int	matrix_qoci (int neq, double hmat, double gmat, int nae, double aemat, double bemat, int nai, double aimat, double bimat, double *xmat)

double *	matrix_bind (int mode, int n11, int n12, double a1, int n21, int n22, double a2, int n31, int n32)

int	matrix_get_extreme (int mode, int ntab, double *tab)

int	matrix_invreal (double *mat, int neq)

int	matrix_invert_triangle (int neq, double *tl, int rank)

void	matrix_tl2tu (int neq, const double tl, double tu)

int	matrix_LU_decompose (int neq, const double a, double tls, double *tus)

int	matrix_LU_solve (int neq, const double tus, const double tls, const double b, double x)

int	matrix_LU_invert (int neq, double *a)

void	matrix_fill_symmetry (int neq, double *a)

Function Documentation

int is_matrix_definite_positive	(	int	neq,
		const double *	a,
		double *	valpro,
		double *	vecpro,
		int	verbose
	)

Check if a matrix is definite positive

Returns: 1 if the matrix is definite positive; 0 otherwise

Parameters

[in]	neq	Size of the matrix
[in]	a	Symmetric square matrix to be checked
[in]	verbose	1 for the verbose option
[out]	valpro	Array of eigen values (Dimension: neq)
[out]	vecpro	Array of eigen vectors (Dimension: neq * neq)

int is_matrix_non_negative	(	int	nrow,
		int	ncol,
		double *	a,
		int	verbose
	)

Check if all the elements of a matrix are non-negative

Returns: 1 if the matrix is non-negative; 0 otherwise

Parameters

[in]	nrow	Number of rows
[in]	ncol	Number of columns
[in]	a	Array to be checked
[in]	verbose	1 for the verbose option

int is_matrix_rotation	(	int	neq,
		const double *	a,
		int	verbose
	)

Check if a matrix is a rotation matrix

Returns: 1 if the matrix is a rotation matrix; 0 otherwise

Parameters

[in]	neq	Size of the matrix
[in]	a	Square matrix to be checked
[in]	verbose	1 for the verbose option

Remarks: A rotation matrix must be orthogonal with determinant equal to 1

int is_matrix_symmetric	(	int	neq,
		const double *	a,
		int	verbose
	)

Check if a matrix is symmetric

Returns: 1 if the matrix is symmetric; 0 otherwise

Parameters

[in]	neq	Size of the matrix
[in]	a	Symmetric square matrix to be checked
[in]	verbose	1 for the verbose option

double* matrix_bind	(	int	mode,
		int	n11,
		int	n12,
		double *	a1,
		int	n21,
		int	n22,
		double *	a2,
		int *	n31,
		int *	n32
	)

Concatenate two matrices

Returns: Pointer on the newly created concatenated matrix (or NULL)

Parameters

[in]	mode	Concatenation type: 1 : by row (n31=n11+n21; n32=n12=n22) 2 : by column (n31=n11=n21; n32=n12+n22)
[in]	n11	First dimension of the first matrix
[in]	n12	Second dimension of the first matrix
[in]	a1	Pointer to the first matrix
[in]	n21	First dimension of the second matrix
[in]	n22	Second dimension of the second matrix
[in]	a2	Pointer to the second matrix
[out]	n31	First dimension of the output matrix
[out]	n32	Second dimension of the output matrix

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

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

void matrix_cholesky_norme	(	int	mode,
		int	neq,
		const double *	tl,
		const double *	a,
		double *	b
	)

Performs the product B = TL * A * TU or TU * A * TL where TL,TU is a triangular matrix and A a square symmetric matrix

Parameters

[in]	mode	0: TL * A * TU; 1: TU * A * TL
[in]	neq	number of equations in the system
[in]	tl	Triangular matrix defined by column
[in]	a	Square symmetric matrix (optional)
[out]	b	Square symmetric matrix

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)

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

Perform a linear combinaison 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[]

double matrix_determinant	(	int	neq,
		const double *	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

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

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

void matrix_fill_symmetry	(	int	neq,
		double *	a
	)

Fill the matrix to take the symmetry into account

Parameters

[in]	neq	matrix dimension
[in,out]	a	square symmetric matrix (dimension = neq*neq)

Remarks: We assume that, in the input matrix, the elements A[i,j] where; i >= j have already been filled

int matrix_get_extreme	(	int	mode,
		int	ntab,
		double *	tab
	)

Find the minimum or maximum value within an array

Returns: Rank of the minimum or maximum value

Parameters

[in]	mode	1 for minimum and 2 for maximum
[in]	ntab	Number of samples in the array
[in]	tab	Array of values

void matrix_int_transpose	(	int	n1,
		int	n2,
		int *	v1,
		int *	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[]

void matrix_int_transpose_in_place	(	int	n1,
		int	n2,
		int *	v1
	)

Transpose an integer matrix in place

Parameters

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

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

int matrix_invert_copy	(	const double *	a,
		int	neq,
		double *	b
	)

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

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

Parameters

[in]	a	input matrix
[in]	neq	number of equations in the matrix 'a'
[out]	b	output matrix

Remarks: The difference with matrix_invert() is that the output; matrix is different from input matrix

int matrix_invert_triangle	(	int	neq,
		double *	tl,
		int	rank
	)

Invert a symmetric square matrix (stored as triangular)

Returns: Error returned code

Parameters

[in,out]	tl	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

int matrix_invgen	(	double *	a,
		int	neq,
		double *	tabout,
		double *	cond
	)

Calculate the generalized inverse of a square symmetric matrix

Returns: Error returned code

Parameters

[in]	a	Symmetric matrix to be inverted
[in]	neq	Number of equations
[out]	tabout	Inverted matrix
[out]	cond	Condition number (MAX(ABS(eigval))/MIN(ABS(eigval)))

Remarks: The input and output matrices can matchprint

int matrix_invreal	(	double *	mat,
		int	neq
	)

Invert a square real full matrix

Returns: Error return code

Parameters

[in]	mat	input matrix, destroyed in computation and replaced by resultant inverse
[in]	neq	number of equations in the matrix 'a'

int matrix_invsvdsym	(	double *	mat,
		int	neq,
		int	rank
	)

Invert a symmetric square matrix (SVD method)

Returns: Error return code

Parameters

[in]	mat	input matrix, destroyed in computation and replaced by resultant inverse
[in]	neq	number of equations in the matrix 'a'
[out]	rank	instability rank (if > 0)

Remarks: This routine has been adapted from a Pascal implementation; (c) 1988 J. C. Nash, "Compact numerical methods for computers",; Hilger 1990.

int matrix_LU_decompose	(	int	neq,
		const double *	a,
		double *	tls,
		double *	tus
	)

LU Decomposition of a square matrix (not necessarily symmetric)

Parameters

neq	Size of the matrix (number of rows or columns)
a	Input square matrix (stored column-wise)
tls	Output square matrix containing lower triangle (stored columnwise)
tus	Output square matrix containing upper triangle (stored linewise)

Remarks: The output matrices 'tus' and 'tls' must be allocated beforehand

int matrix_LU_invert	(	int	neq,
		double *	a
	)

int matrix_LU_solve	(	int	neq,
		const double *	tus,
		const double *	tls,
		const double *	b,
		double *	x
	)

void matrix_manage	(	int	nrows,
		int	ncols,
		int	nr,
		int	nc,
		int *	rowsel,
		int *	colsel,
		double *	v1,
		double *	v2
	)

Manage a matrix and derive a submatrix

Parameters

[in]	nrows	Number of rows of the input matrix
[in]	ncols	Number of columns of the input matrix
[in]	nr	Number of rows of interest >0 : for extraction 0 : no action <0 : for suppression
[in]	nc	Number of columns of interest >0 : for extraction 0 : no action <0 : for suppression
[in]	rowsel	Array of rows indices of interest (starting from 0) (Dimension: ABS(nr)
[in]	colsel	Array of columns indices of interest (starting from 0) (Dimension: ABS(nr)
[in]	v1	Input rectangular matrix (Dimension: nrows * ncols)
[out]	v2	Output rectangular matrix

double matrix_normA	(	double *	b,
		double *	a,
		int	neq,
		int	subneq
	)

Calculate the square norm of a vector issued from the inner product relative to a matrix A

Returns: Value of the norm

Parameters

[in]	b	Vector (Dimension: neq)
[in]	a	Matrix (Symmetric, Dimension neq x neq)
[in]	neq	Space dimension
[in]	subneq	Dimension of the subspace to be considered

double matrix_norminf	(	int	nval,
		double *	tab
	)

Calculate the maximum of the absolute values of a vector

Parameters

[in]	nval	vector dimension
[in]	tab	vector

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)

void matrix_product	(	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 coincide with one of the two initial ones

void matrix_product_by_diag	(	int	mode,
		int	neq,
		double *	a,
		double *	c,
		double *	b
	)

Performs the A %*% diag(c) where A is a square matrix and c a vector

Parameters

[in]	mode	0: c as is; 1: sqrt(c); 2: 1/c; 3: 1/sqrt(c)
[in]	neq	matrix dimension for A
[in]	a	square matrix
[in]	c	vector
[out]	b	square matrix

Remarks: Matrices a() and b() may coincide

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

void matrix_produit_cholesky	(	int	neq,
		const double *	tl,
		double *	a
	)

Calculate the product of 'tl' (lower triangle) by its transpose

Parameters

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

VectorDouble matrix_produit_cholesky_VD	(	int	neq,
		const double *	tl
	)

Calculate the product of 'tl' (lower triangle) by its transpose and store the result in a VectorDouble

Returns: The resulting VectorDouble

Parameters

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

int matrix_qo	(	int	neq,
		double *	hmat,
		double *	gmat,
		double *	xmat
	)

Solve a linear system: H %*% g = x

Returns: Error return code

Parameters

[in]	neq	number of equations in the system
[in]	hmat	symmetric square matrix (Dimension: neq * neq)
[in]	gmat	right-hand side vector (Dimension: neq)
[out]	xmat	solution vector (Dimension: neq)

Remarks: In output, hmat contains the inverse matrix

int matrix_qoc	(	int	flag_invert,
		int	neq,
		double *	hmat,
		double *	gmat,
		int	na,
		double *	amat,
		double *	bmat,
		double *	xmat,
		double *	lambda
	)

Minimize 1/2 t(x) %*% H %*% x + t(g) %*% x under the constraints t(A) %*% x = b

Returns: Error return code

Parameters

[in]	flag_invert	Tells if the inverse has already been calculated
[in]	neq	number of equations in the system
[in]	hmat	symmetric square matrix (Dimension: neq * neq)
[in]	gmat	right-hand side vector (Dimension: neq)
[in]	na	Number of equalities
[in]	amat	matrix for inequalities (Dimension: neq * na)
[in]	bmat	inequality vector (Dimension: na)
[in]	xmat	solution of the linear system with no constraint. On return, solution with constraints (Dimension: neq)
[out]	lambda	working vector (Dimension: na)

Remarks: In input:; If flag_invert== 1, H is provided as the generalized inverse; and x contains the solution of the linear system with no constraint; If flag_invert==0, H is the primal matrix; In output, H contains the inverse matrix

int matrix_qoci	(	int	neq,
		double *	hmat,
		double *	gmat,
		int	nae,
		double *	aemat,
		double *	bemat,
		int	nai,
		double *	aimat,
		double *	bimat,
		double *	xmat
	)

Minimize 1/2 t(x) %*% H %*% x + t(g) %*% x under the constraints t(Ae) %*% x = be and t(Ai) %*% x = bi

Returns: Error return code

Parameters

[in]	neq	number of equations in the system
[in]	hmat	symmetric square matrix (Dimension: neq * neq)
[in]	gmat	right-hand side vector (Dimension: neq)
[in]	nae	Number of equalities
[in]	aemat	matrix for equalities (Dimension: neq * nae)
[in]	bemat	right-hand side for equalities (Dimension: nae)
[in]	nai	Number of inequalities
[in]	aimat	matrix for inequalities (Dimension: neq * nai)
[in]	bimat	right-hand side for inequalities (Dimension: nai)
[in,out]	xmat	solution of the linear system with constraints (neq)

REMARKS: The initial xmat has to satisfied all the constraints.

int matrix_solve	(	int	mode,
		const double *	a,
		const double *	b,
		double *	x,
		int	neq,
		int	nrhs,
		int *	pivot
	)

Solve a system of linear equations with symmetric coefficient matrix upper triangular part of which is stored columnwise. Use is restricted to matrices whose leading principal minors have non-zero determinants.

Returns: Error return code : 1 for core problem; 0 otherwise

Parameters

[in]	mode	Storage: 0 for upper triangular, 1 for lower triangular
[in]	a	triangular matrix (dimension = neq*(neq+1)/2)
[in]	b	right-hand side matrix (dimension = neq*nrhs)
[in]	neq	number of equations in the system
[in]	nrhs	number of right-hand side vectors
[out]	x	matrix of solutions (dimension = neq*nrhs)
[out]	pivot	rank of the pivoting error (0 none)

void matrix_svd_inverse	(	int	neq,
		double *	s,
		double *	u,
		double *	v,
		double *	tabout
	)

Invert a symmetric square matrix (SVD method)

Parameters

[in]	neq	number of equations in the matrix 'a'
[in]	s	input matrix (square)
[out]	u	first decomposition matrix
[out]	v	second decomposition matrix
[out]	tabout	inverse matrix

Remarks: This routine has been adapted from a Pascal implementation; (c) 1988 J. C. Nash, "Compact numerical methods for computers",; Hilger 1990.

void matrix_tl2tu	(	int	neq,
		const double *	tl,
		double *	tu
	)

Change the storage of triangular matrices (stored columnwise)

Parameters

[in]	neq	Matrix dimension
[in]	tl	Lower Triangular Matrix
[in]	tu	Upper Triangular Matrix

void matrix_transpose	(	int	n1,
		int	n2,
		double *	v1,
		double *	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[]

void matrix_transpose_in_place	(	int	n1,
		int	n2,
		double *	v1
	)

Transpose a matrix in place

Parameters

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

void matrix_triangle_to_square	(	int	mode,
		int	neq,
		const double *	tl,
		double *	a
	)

Fill a square matrix with a triangular matrix

Parameters

[in]	mode	0: TL (upper); 1: TL (lower)
[in]	neq	number of equations in the system
[in]	tl	Triangular matrix (lower part)
[out]	a	Resulting square matrix

Functions

Function Documentation