cs_math.c File Reference

#include "cs_defs.h"
#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <float.h>
#include <mpi.h>
#include "bft_error.h"
#include "bft_mem.h"
#include "cs_math.h"

Include dependency graph for cs_math.c:

Functions
void	cs_f_math_sym_33_inv_cramer (const cs_real_t s[6], cs_real_t sout[6])
void	cs_f_math_sym_33_product (const cs_real_t s1[6], const cs_real_t s2[6], cs_real_t sout[6])
void	cs_f_math_reduce_sym_prod_33_to_66 (const cs_real_t s[3][3], cs_real_t sout[6][6])
void	cs_f_math_3_normalize (const cs_real_t vin[3], cs_real_t vout[3])
void	cs_math_sym_33_eigen (const cs_real_t m[6], cs_real_t eig_vals[3])
	Compute all eigenvalues of a 3x3 symmetric matrix with symmetric storage. More...
void	cs_math_33_eigen (const cs_real_t m[3][3], cs_real_t *eig_ratio, cs_real_t *eig_max)
	Compute max/min eigenvalues ratio and max. eigenvalue of a 3x3 symmetric matrix with non-symmetric storage. More...
void	cs_math_3_length_unitv (const cs_real_t xa[3], const cs_real_t xb[3], cs_real_t *len, cs_real_3_t unitv)
	Compute the length (Euclidean norm) between two points xa and xb in a Cartesian coordinate system of dimension 3. More...
double	cs_math_surftri (const cs_real_t xv[3], const cs_real_t xe[3], const cs_real_t xf[3])
	Compute the area of the convex_hull generated by 3 points. This corresponds to the computation of the surface of a triangle. More...
double	cs_math_voltet (const cs_real_t xv[3], const cs_real_t xe[3], const cs_real_t xf[3], const cs_real_t xc[3])
	Compute the volume of the convex_hull generated by 4 points. This is equivalent to the computation of the volume of a tetrahedron. More...
void	cs_math_33_eig_val_vec (const cs_real_t m_in[3][3], const cs_real_t tol_err, cs_real_t eig_val[restrict 3], cs_real_t eig_vec[restrict 3][3])
	Evaluate eigenvalues and eigenvectors of a real symmetric matrix m1[3,3]: m1*m2 = lambda*m2. More...
void	cs_math_fact_lu (cs_lnum_t n_blocks, int b_size, const cs_real_t *a, cs_real_t *a_lu)
	Compute LU factorization of an array of dense matrices of identical size. More...
void	cs_math_fw_and_bw_lu (const cs_real_t a_lu[], int n, cs_real_t x[], const cs_real_t b[])
	Block Jacobi utilities. Compute forward and backward to solve an LU P*P system. More...
void	cs_math_sym_44_factor_ldlt (cs_real_t ldlt[10])
	LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http://mathforcollege.com/nm/mws/gen/04sle/mws_gen_sle_txt_cholesky.pdf. More...
cs_real_t	cs_math_sym_44_partial_solve_ldlt (const cs_real_t ldlt[10], const cs_real_t rhs[4])
	LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http://mathforcollege.com/nm/mws/gen/04sle/mws_gen_sle_txt_cholesky.pdf. More...

Detailed Description

Mathematical base functions.

Function Documentation

cs_f_math_3_normalize()

void cs_f_math_3_normalize	(	const cs_real_t	vin[3],
		cs_real_t	vout[3]
	)

cs_f_math_reduce_sym_prod_33_to_66()

void cs_f_math_reduce_sym_prod_33_to_66	(	const cs_real_t	s[3][3],
		cs_real_t	sout[6][6]
	)

cs_f_math_sym_33_inv_cramer()

void cs_f_math_sym_33_inv_cramer	(	const cs_real_t	s[6],
		cs_real_t	sout[6]
	)

cs_f_math_sym_33_product()

void cs_f_math_sym_33_product	(	const cs_real_t	s1[6],
		const cs_real_t	s2[6],
		cs_real_t	sout[6]
	)

cs_math_33_eig_val_vec()

void cs_math_33_eig_val_vec	(	const cs_real_t	m_in[3][3],
		const cs_real_t	tol_err,
		cs_real_t	eig_val[restrict 3],
		cs_real_t	eig_vec[restrict 3][3]
	)

Evaluate eigenvalues and eigenvectors of a real symmetric matrix m1[3,3]: m1*m2 = lambda*m2.

Use of Jacobi method for symmetric matrices (adapted from the book Numerical Recipes in C, Chapter 11.1)

Parameters

[in]	m_in	matrix of 3x3 real values (initial)
[in]	tol_err	absolute tolerance (sum of off-diagonal elements)
[out]	eig_val	vector of 3 real values (eigenvalues)
[out]	eig_vec	matrix of 3x3 real values (eigenvectors)

cs_math_33_eigen()

void cs_math_33_eigen	(	const cs_real_t	m[3][3],
		cs_real_t *	eig_ratio,
		cs_real_t *	eig_max
	)

Compute max/min eigenvalues ratio and max. eigenvalue of a 3x3 symmetric matrix with non-symmetric storage.

Based on: Oliver K. Smith "eigenvalues of a symmetric 3x3 matrix", Communication of the ACM (April 1961) (Wikipedia article entitled "Eigenvalue algorithm")

Parameters

[in]	m	3x3 matrix
[out]	eig_ratio	max/min
[out]	eig_max	max. eigenvalue

cs_math_3_length_unitv()

void cs_math_3_length_unitv ( const cs_real_t  xa[3], const cs_real_t  xb[3], cs_real_t *  len, cs_real_3_t  unitv  )

inline

Compute the length (Euclidean norm) between two points xa and xb in a Cartesian coordinate system of dimension 3.

Parameters

[in]	xa	coordinate of the first extremity
[in]	xb	coordinate of the second extremity
[out]	len	pointer to the length of the vector va -> vb
[out]	unitv	unitary vector along va -> vb

cs_math_fact_lu()

void cs_math_fact_lu	(	cs_lnum_t	n_blocks,
		int	b_size,
		const cs_real_t *	a,
		cs_real_t *	a_lu
	)

Compute LU factorization of an array of dense matrices of identical size.

Parameters

[in]	n_blocks	number of blocks
[in]	b_size	block size
[in]	a	matrix blocks
[out]	a_lu	LU factorizations of matrix blocks

cs_math_fw_and_bw_lu()

void cs_math_fw_and_bw_lu	(	const cs_real_t	a_lu[],
		int	n,
		cs_real_t	x[],
		const cs_real_t	b[]
	)

Block Jacobi utilities. Compute forward and backward to solve an LU P*P system.

Parameters

[in]	a_lu	matrix LU factorization
[in]	n	matrix size
[out]	x	solution
[out]	b	right hand side

cs_math_surftri()

double cs_math_surftri ( const cs_real_t  xv[3], const cs_real_t  xe[3], const cs_real_t  xf[3]  )

inline

Compute the area of the convex_hull generated by 3 points. This corresponds to the computation of the surface of a triangle.

Parameters

[in]	xv	coordinates of the first vertex
[in]	xe	coordinates of the second vertex
[in]	xf	coordinates of the third vertex

Returns

the surface of a triangle

cs_math_sym_33_eigen()

void cs_math_sym_33_eigen	(	const cs_real_t	m[6],
		cs_real_t	eig_vals[3]
	)

Compute all eigenvalues of a 3x3 symmetric matrix with symmetric storage.

Based on: Oliver K. Smith "eigenvalues of a symmetric 3x3 matrix", Communication of the ACM (April 1961) (Wikipedia article entitled "Eigenvalue algorithm")

Parameters

[in]	m	3x3 symmetric matrix (m11, m22, m33, m12, m23, m13)
[out]	eig_vals	size 3 vector

cs_math_sym_44_factor_ldlt()

void cs_math_sym_44_factor_ldlt

(

cs_real_t

ldlt[10]

)

LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http://mathforcollege.com/nm/mws/gen/04sle/mws_gen_sle_txt_cholesky.pdf.

Parameters

[in,out]

ldlt

pointer to matrix coefficients: (m00, m10, m11, m20, m21, m22, m30, m31, m32, m33) in (f00, l10, f11, l20, l21, f22, l30, l31, l32, f33) out

cs_math_sym_44_partial_solve_ldlt()

cs_real_t cs_math_sym_44_partial_solve_ldlt	(	const cs_real_t	ldlt[10],
		const cs_real_t	rhs[4]
	)

LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http://mathforcollege.com/nm/mws/gen/04sle/mws_gen_sle_txt_cholesky.pdf.

Here we only need to use the last element of the solution vector, so we return that value only and simplify the computation.

Parameters

[in,out]	ldlt	pointer to matrix coefficients: (f00, l10, f11, l20, l21, f22, l30, l31, l32, f33)
[in]	rhs	right-hand side

cs_math_voltet()

double cs_math_voltet	(	const cs_real_t	xv[3],
		const cs_real_t	xe[3],
		const cs_real_t	xf[3],
		const cs_real_t	xc[3]
	)

Compute the volume of the convex_hull generated by 4 points. This is equivalent to the computation of the volume of a tetrahedron.

Parameters

[in]	xv	coordinates of the first vertex
[in]	xe	coordinates of the second vertex
[in]	xf	coordinates of the third vertex
[in]	xc	coordinates of the fourth vertex

Returns

the volume of the tetrahedron.