#include <vector>
#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

Classes
class	Matrix
	Matrix used in continuum membrane model. More...

Functions
Matrix	mat_calloc (const int &nrow, const int &ncol)
	This function allocates a new Matrix object with the specified number of rows and columns, and initializes all elements to zero.

std::ostream &	operator<< (std::ostream &stream, const Matrix &matrix)
	Bit push all element of a matrix in the format of:

Matrix &	operator+= (Matrix &m1, const Matrix &m2)
	This implementation takes two gsl_matrix pointers as arguments and adds the elements of the second matrix to the first matrix. It first checks that the matrices have the same dimensions, and throws an exception if they do not. Then it iterates over the elements of the matrices and adds the corresponding elements. Finally, it returns a pointer to the modified first matrix. Note that this implementation modifies the first matrix in place, rather than creating a new matrix. If you want to create a new matrix instead, you could modify the implementation accordingly.

Matrix	operator+ (const Matrix &m1, const Matrix &m2)
	This implementation checks that the matrices have the same dimensions and allocates memory for the result of the addition. It then adds the matrices element-wise and stores the result in a newly allocated matrix result, which is then returned.

Matrix &	operator-= (Matrix &m1, const Matrix &m2)
	This implementation takes two gsl_matrix pointers as arguments and subtract the elements of the second matrix from the first matrix. It first checks that the matrices have the same dimensions, and throws an exception if they do not. Then it iterates over the elements of the matrices and adds the corresponding elements. Finally, it returns a pointer to the modified first matrix. Note that this implementation modifies the first matrix in place, rather than creating a new matrix. If you want to create a new matrix instead, you could modify the implementation accordingly.

Matrix	operator- (const Matrix &m1, const Matrix &m2)
	This implementation checks that the matrices have the same dimensions and allocates memory for the result of the subtraction. It then adds the matrices element-wise and stores the result in a newly allocated matrix result, which is then returned.

Matrix &	operator*= (Matrix &mat, const double &scalar)
	This implementation uses the gsl_matrix_scale function to scale each element in the matrix by the scalar value. The function returns the modified matrix, which allows for method chaining when using the *= operator.

Matrix &	operator*= (Matrix m1, const Matrix &m2)
	This implementation checks that the matrices have compatible dimensions for multiplication and uses the BLAS library to perform the matrix multiplication. The result is stored in a temporary matrix, which is then copied back into m1. Finally, the temporary memory is freed before returning m1.

Matrix &	operator/= (Matrix &mat, const double &scalar)
	This implementation uses the gsl_matrix_scale function to scale each element in the matrix by 1/the scalar value. The function returns the modified matrix, which allows for method chaining when using the *= operator.

Matrix	operator/ (Matrix matrix, const double &scalar)
	This operator takes a double scalar and a matrix as its operands. It returns a new gsl_matrix* that is the result of scaling the input matrix by 1/scalar.

double	dot_col (const Matrix &m1, const Matrix &m2)
	This implementation takes two column vectors as arguments and returns their dot product as double.

double	dot_row (const Matrix &m1, const Matrix &m2)
	This implementation takes two row vectors as arguments and returns their dot product as double.

Matrix	cross_col (const Matrix &m1, const Matrix &m2)
	This implementation takes two column vectors as arguments and returns their cross product as a new column vector.

Matrix	cross_row (const Matrix &m1, const Matrix &m2)
	This implementation takes two row vectors as arguments and returns their cross product as a new row vector.

void	negative (const Matrix &m1, Matrix &m_neg)
	Computes the negative of a given matrix.

void	get_unit_vector (Matrix &m)
	Get the unit vector in the direction of m which is assumed to be a (3, 1) Matrix and overwrites the result to m.

void	get_unit_vector (const Matrix &m1, Matrix &m_unit)
	Get the unit vector in the direction of m1 which is a (3, 1) Matrix.

void	cross (const Matrix &m1, const Matrix &m2, Matrix &temp)
	Compute the cross product of two 3D vectors represented as 3x1 matrices.

void	addition (const Matrix &m1, const Matrix &m2, Matrix &tmp)
	Compute the element-wise addition of two matrices and store the result in a third matrix.

void	subtraction (const Matrix &m1, const Matrix &m2, Matrix &tmp)
	Compute the element-wise subtraction of two matrices and store the result in a third matrix.

void	multiplication (const Matrix &m1, const Matrix &m2, Matrix &tmp)
	Multiply a matrix by a constant and store the result in a second matrix.

void	const_multiplication (const Matrix &m1, const double num, Matrix &tmp)
	Multiply a matrix by a constant and store the result in a second matrix.

void	const_division (const Matrix &m1, const double num, Matrix &tmp)
	Divide a matrix by a constant and store the result in a second matrix.

void	colvec_matrix_multiplication (const Matrix &v1, const Matrix &m1, Matrix &tmp)
	Multiply a row vector by a matrix and store the result in a row vector.

Matrix	kron (const Matrix &v1, const Matrix &v2)
	the Kronecker product between two column matrices and returns the result in an output matrix. The Kronecker product is defined as follows:

Matrix	kron (const Matrix &v1, const Matrix &v2, Matrix &m_result)
	the Kronecker product between two column matrices and returns the result in an output matrix. The Kronecker product is defined as follows:

void	assign_rowVec_to_colVec (const Matrix &srcRowVec, Matrix &destColVec)

void	a_cross_b_plus_c_cross_d (const Matrix &a, const Matrix &b, const Matrix &c, const Matrix &d, Matrix &tmp_f, Matrix &tmp_l, Matrix &v_result)


Matrix	operator* (const Matrix &matrix, const double &scalar)
	This operator takes a double scalar and a matrix as its operands. It returns a new gsl_matrix* that is the result of scaling the input matrix by the scalar.

Matrix	operator* (const double &scalar, const Matrix &matrix)


Matrix	operator* (const Matrix &m1, const Matrix &m2)
	This implementation checks that the matrices have compatible dimensions for multiplication and uses the BLAS library to perform the matrix multiplication. The result is stored in a newly allocated matrix result, which is then returned. The temporary memory is not freed here, and it is the responsibility of the caller to free the memory once the result is no longer needed.

Matrix	dot (const Matrix &m1, const Matrix &m2)

Detailed Description

This file serves as an interface for implementing any linear algebra package for the continuum membrane model.

The code here defines a Matrix class that uses the GNU Scientific Library (GSL) to implement various linear algebra operations. The class has four constructors:

Warning: this interface between Linalg package (currently GSL linalg) mutes all pre-operation checks to enhance the computational efficiency as much as possible. However, this may also results in "core-dumped" or "double-free" error if the class functions are not handled correctly. Please make sure to read thru the docstring, especially the prerequistes and warnings carefully before you reference to the functions!

Todo:: benchmark - In this case, if my matrices are small (< 12 * 3), but I need to repeat these matrices calculation for a large number of times, is gsl_blas_dgemm still more eifficent?

For small matrices like the ones you've described, it's possible that an element-wise implementation could be faster than using gsl_blas_dgemm() due to the overhead associated with setting up the function call and memory management. However, this depends on the context in which the function is being used and the specific hardware being used.

If the calculation of these small matrices is being repeated a large number of times, then the overhead of setting up the function call and memory management might be amortized over multiple calculations, making gsl_blas_dgemm() more efficient.

Ultimately, the best approach would be to benchmark both implementations (element-wise and using gsl_blas_dgemm() on your specific system and compare their performance for your use case.

-ChatGPT

Date: 2023-03-20 (Created); 2024-01-29 (Edited)

Copyright: Copyright (c) 2023

Function Documentation

◆ a_cross_b_plus_c_cross_d()

void a_cross_b_plus_c_cross_d	(	const Matrix &	a,
		const Matrix &	b,
		const Matrix &	c,
		const Matrix &	d,
		Matrix &	tmp_f,
		Matrix &	tmp_l,
		Matrix &	v_result
	)

Computes the sum of two cross products and stores the result in an output matrix.

Parameters

[in]	a	First input column matrix.
[in]	b	Second input column matrix.
[in]	c	Third input column matrix.
[in]	d	Fourth input column matrix.
[in,out]	tmp_f	Temporary matrix for storing the first cross product.
[in,out]	tmp_l	Temporary matrix for storing the second cross product.
[out]	v_result	Output matrix to store the final result.

Exceptions

std::invalid_argument if the input matrices have incompatible dimensions.

◆ addition()

void addition	(	const Matrix &	m1,
		const Matrix &	m2,
		Matrix &	tmp
	)

Compute the element-wise addition of two matrices and store the result in a third matrix.

This function computes the element-wise addition of two input matrices of equal dimensions and stores the result in a third matrix. All input matrices are assumed to be of type Matrix, which internally contains a gsl_matrix pointer for efficiency.

Parameters

m1	The first input matrix.
m2	The second input matrix.
tmp	The output matrix to store the resulting sum.

Precondition: All input matrices must have the same dimensions.

Postcondition: The output matrix 'tmp' contains the element-wise sum of 'm1' and 'm2'.

◆ assign_rowVec_to_colVec()

void assign_rowVec_to_colVec	(	const Matrix &	srcRowVec,
		Matrix &	destColVec
	)

Copies the values in a row vector to a column vector.

Parameters

[in]	srcRowVec	Row vector containing the values to be copied.
[in,out]	destColVec	Column vector to receive the copied values.

Exceptions

std::invalid_argument if the input matrices have incompatible dimensions.

◆ colvec_matrix_multiplication()

void colvec_matrix_multiplication	(	const Matrix &	v1,
		const Matrix &	m1,
		Matrix &	tmp
	)

Multiply a row vector by a matrix and store the result in a row vector.

This function multiplies a row vector with a matrix and stores the resulting product in another row vector. All input matrices are assumed to be of type Matrix, which internally contains a gsl_matrix pointer for efficiency.

Parameters

v1	The input row vector to be multiplied.
m1	The input matrix.
tmp	The output row vector to store the resulting product.

Precondition: The number of columns of 'v1' must match the number of rows of 'm1', and the dimensions of 'tmp' must match the number of columns of 'm1'.

Postcondition: The output row vector 'tmp' contains the product of the input row vector 'v1' and the input matrix 'm1'.

Todo:: benchmark - In this case, if my matrices are small (< 12 * 3), but I need to repeat these matrices calculation for a large number of times, is gsl_blas_dgemm still more eifficent?

For small matrices like the ones you've described, it's possible that an element-wise implementation could be faster than using gsl_blas_dgemm() due to the overhead associated with setting up the function call and memory management. However, this depends on the context in which the function is being used and the specific hardware being used.

If the calculation of these small matrices is being repeated a large number of times, then the overhead of setting up the function call and memory management might be amortized over multiple calculations, making gsl_blas_dgemm() more efficient.

Ultimately, the best approach would be to benchmark both implementations (element-wise and using gsl_blas_dgemm() on your specific system and compare their performance for your use case.

◆ const_division()

void const_division	(	const Matrix &	m1,
		const double	num,
		Matrix &	tmp
	)

Divide a matrix by a constant and store the result in a second matrix.

This function divides the elements of an input matrix by a scalar constant and stores the result in a second output matrix. All input matrices are assumed to be of type Matrix.

Parameters

m1	The input matrix to be divided.
num	The scalar constant to divide by.
tmp	The output matrix to store the resulting quotient.

Precondition: The dimensions of 'm1' and 'tmp' must match.

Postcondition: The output matrix 'tmp' contains the element-wise division of 'm1' and 'num'.

◆ const_multiplication()

void const_multiplication	(	const Matrix &	m1,
		const double	num,
		Matrix &	tmp
	)

Multiply a matrix by a constant and store the result in a second matrix.

This function multiplies the elements of an input matrix with a scalar constant and stores the result in a second output matrix. All input matrices are assumed to be of type Matrix

Parameters

m1	The input matrix to be multiplied.
num	The scalar constant to multiply with.
tmp	The output matrix to store the resulting product.

Precondition: The dimensions of 'm1' and 'tmp' must match.

Postcondition: The output matrix 'tmp' contains the element-wise multiplication of 'm1' and 'num'.

◆ cross()

void cross	(	const Matrix &	m1,
		const Matrix &	m2,
		Matrix &	temp
	)

Compute the cross product of two 3D vectors represented as 3x1 matrices.

This function computes the cross product of two input vectors, which are assumed to be 3x1 matrices. The result is stored in another 3x1 matrix passed in as an argument.

Parameters

m1	The first input vector as a 3x1 matrix.
m2	The second input vector as a 3x1 matrix.
temp	The output matrix to store the resulting cross product.

Precondition: All input matrices must have 3 rows and 1 column.

Postcondition: The output matrix 'temp' contains the cross product of 'm1' and 'm2'.

◆ cross_col()

Matrix cross_col	(	const Matrix &	m1,
		const Matrix &	m2
	)

This implementation takes two column vectors as arguments and returns their cross product as a new column vector.

Note: Note that this implementation assumes that the input matrices are already column vectors, so no additional transposition is necessary.

◆ cross_row()

Matrix cross_row	(	const Matrix &	m1,
		const Matrix &	m2
	)

This implementation takes two row vectors as arguments and returns their cross product as a new row vector.

Note: Note that this implementation assumes that the input matrices are already row vectors, so no additional transposition is necessary.

◆ dot()

Matrix dot	(	const Matrix &	m1,
		const Matrix &	m2
	)

◆ dot_col()

double dot_col	(	const Matrix &	m1,
		const Matrix &	m2
	)

This implementation takes two column vectors as arguments and returns their dot product as double.

Note: Note that this implementation assumes that the input matrices are already column vectors, so no additional transposition is necessary.

◆ dot_row()

double dot_row	(	const Matrix &	m1,
		const Matrix &	m2
	)

This implementation takes two row vectors as arguments and returns their dot product as double.

Note: Note that this implementation converts the row vectors into column vectors before computing the dot product.

◆ get_unit_vector() [1/2]

void get_unit_vector	(	const Matrix &	m1,
		Matrix &	m_unit
	)

Get the unit vector in the direction of m1 which is a (3, 1) Matrix.

This function computes the unit vector in the direction of a given input vector, which is assumed to be a 3x1 matrix. The result is stored in another 3x1 matrix.

Warning: m1 and m_unit cannot be the same Matrix instances! Use get_unit_vector(Matrix &m) instead.

Parameters

m1	The input vector to compute the unit vector of.
m_unit	The output matrix to store the resulting unit vector.

Precondition: The input matrix 'm1' must be a 3x1 matrix.

Postcondition: The output matrix 'm_unit' contains the unit vector in the direction of 'm1'.

◆ get_unit_vector() [2/2]

void get_unit_vector ( Matrix & m )

Get the unit vector in the direction of m which is assumed to be a (3, 1) Matrix and overwrites the result to m.

Parameters

m1	The input vector to compute the unit vector of and the output matrix to store the resulting unit vector.

Precondition: The input matrix 'm' must be a 3x1 matrix.

Postcondition: The output matrix 'm' contains the unit vector in the direction of 'm'.

◆ kron() [1/2]

Matrix kron	(	const Matrix &	v1,
		const Matrix &	v2
	)

the Kronecker product between two column matrices and returns the result in an output matrix. The Kronecker product is defined as follows:

kron(v1, v2) = [v1(1)*v2(1)  v1(1)*v2(2)  ...  v1(1)*v2(n2)
                v1(2)*v2(1)  v1(2)*v2(2)  ...  v1(2)*v2(n2)
                ...          ...          ...  ...
                v1(n1)*v2(1)  v1(n1)*v2(2)  ...  v1(n1)*v2(n2)]

where v1 and v2 are column matrices of size n1 and n2 respectively, and the output matrix has size n1*n2 x n1*n2.

Parameters

[in]	v1	First input column matrix.
[in]	v2	Second input column matrix.

Deprecated:: Might cause memory leak!

Returns: Output matrix to store the result.

Exceptions

std::invalid_argument if the output matrix has incompatible dimensions.

◆ kron() [2/2]

Matrix kron	(	const Matrix &	v1,
		const Matrix &	v2,
		Matrix &	m_result
	)

the Kronecker product between two column matrices and returns the result in an output matrix. The Kronecker product is defined as follows:

kron(v1, v2) = [v1(1)*v2(1)  v1(1)*v2(2)  ...  v1(1)*v2(n2)
                v1(2)*v2(1)  v1(2)*v2(2)  ...  v1(2)*v2(n2)
                ...          ...          ...  ...
                v1(n1)*v2(1)  v1(n1)*v2(2)  ...  v1(n1)*v2(n2)]

where v1 and v2 are column matrices of size n1 and n2 respectively, and the output matrix has size n1*n2 x n1*n2.

Parameters

[in]	v1	First input column matrix.
[in]	v2	Second input column matrix.
[out]	m_result	Output matrix to store the result.

◆ mat_calloc()

Matrix mat_calloc	(	const int &	nrow,
		const int &	ncol
	)

This function allocates a new Matrix object with the specified number of rows and columns, and initializes all elements to zero.

Parameters

nrow
ncol

Returns: Matrix

◆ multiplication()

void multiplication	(	const Matrix &	m1,
		const Matrix &	m2,
		Matrix &	tmp
	)

Multiply a matrix by a constant and store the result in a second matrix.

This function multiplies the elements of an input matrix with a scalar constant and stores the result in a second output matrix. All input matrices are assumed to be of type Matrix

Parameters

m1	The input matrix to be multiplied.
num	The scalar constant to multiply with.
tmp	The output matrix to store the resulting product.

Precondition: The dimensions of 'm1' and 'tmp' must match.

Postcondition: The output matrix 'tmp' contains the element-wise multiplication of 'm1' and 'num'.

◆ negative()

void negative	(	const Matrix &	m1,
		Matrix &	m_neg
	)

Computes the negative of a given matrix.

This function computes the element-wise negative of a given matrix and stores the result in another matrix. Element-wise negative is defined as multiplying each element of the matrix with -1. The input matrix 'm1' is not modified by this function.

Parameters

m1	The matrix to compute the negative of.
m_neg	The matrix to store the result in.

Precondition: Both matrices 'm1' and 'm_neg' must have been allocated memory and initialized properly.; Both matrices 'm1' and 'm_neg' must have the same dimensions.

Postcondition: The matrix 'm_neg' contains the element-wise negative of 'm1'.

◆ operator*() [1/3]

Matrix operator*	(	const double &	scalar,
		const Matrix &	matrix
	)

◆ operator*() [2/3]

Matrix operator*	(	const Matrix &	m1,
		const Matrix &	m2
	)

This implementation checks that the matrices have compatible dimensions for multiplication and uses the BLAS library to perform the matrix multiplication. The result is stored in a newly allocated matrix result, which is then returned. The temporary memory is not freed here, and it is the responsibility of the caller to free the memory once the result is no longer needed.

Parameters

m1
m2

Returns: Matrix

◆ operator*() [3/3]

Matrix operator*	(	const Matrix &	matrix,
		const double &	scalar
	)

This operator takes a double scalar and a matrix as its operands. It returns a new gsl_matrix* that is the result of scaling the input matrix by the scalar.

Parameters

scalar
matrix

Returns: Matrix

◆ operator*=() [1/2]

Matrix & operator*=	(	Matrix &	mat,
		const double &	scalar
	)

This implementation uses the gsl_matrix_scale function to scale each element in the matrix by the scalar value. The function returns the modified matrix, which allows for method chaining when using the *= operator.

Parameters

mat
scalar

Returns: Matrix

◆ operator*=() [2/2]

Matrix & operator*=	(	Matrix	m1,
		const Matrix &	m2
	)

This implementation checks that the matrices have compatible dimensions for multiplication and uses the BLAS library to perform the matrix multiplication. The result is stored in a temporary matrix, which is then copied back into m1. Finally, the temporary memory is freed before returning m1.

Parameters

m1
m2

Returns: Matrix

◆ operator+()

Matrix operator+	(	const Matrix &	m1,
		const Matrix &	m2
	)

This implementation checks that the matrices have the same dimensions and allocates memory for the result of the addition. It then adds the matrices element-wise and stores the result in a newly allocated matrix result, which is then returned.

◆ operator+=()

Matrix & operator+=	(	Matrix &	m1,
		const Matrix &	m2
	)

This implementation takes two gsl_matrix pointers as arguments and adds the elements of the second matrix to the first matrix. It first checks that the matrices have the same dimensions, and throws an exception if they do not. Then it iterates over the elements of the matrices and adds the corresponding elements. Finally, it returns a pointer to the modified first matrix. Note that this implementation modifies the first matrix in place, rather than creating a new matrix. If you want to create a new matrix instead, you could modify the implementation accordingly.

◆ operator-()

Matrix operator-	(	const Matrix &	m1,
		const Matrix &	m2
	)

This implementation checks that the matrices have the same dimensions and allocates memory for the result of the subtraction. It then adds the matrices element-wise and stores the result in a newly allocated matrix result, which is then returned.

◆ operator-=()

Matrix & operator-=	(	Matrix &	m1,
		const Matrix &	m2
	)

This implementation takes two gsl_matrix pointers as arguments and subtract the elements of the second matrix from the first matrix. It first checks that the matrices have the same dimensions, and throws an exception if they do not. Then it iterates over the elements of the matrices and adds the corresponding elements. Finally, it returns a pointer to the modified first matrix. Note that this implementation modifies the first matrix in place, rather than creating a new matrix. If you want to create a new matrix instead, you could modify the implementation accordingly.

◆ operator/()

Matrix operator/	(	Matrix	matrix,
		const double &	scalar
	)

This operator takes a double scalar and a matrix as its operands. It returns a new gsl_matrix* that is the result of scaling the input matrix by 1/scalar.

Parameters

scalar
matrix

Returns: Matrix

◆ operator/=()

Matrix & operator/=	(	Matrix &	mat,
		const double &	scalar
	)

This implementation uses the gsl_matrix_scale function to scale each element in the matrix by 1/the scalar value. The function returns the modified matrix, which allows for method chaining when using the *= operator.

Parameters

mat
scalar

Returns: Matrix

◆ operator<<()

std::ostream & operator<<	(	std::ostream &	stream,
		const Matrix &	matrix
	)

Bit push all element of a matrix in the format of:

a11, a12, a13, ..., a21, a22, a23, ..., ...,

Parameters

stream	ostream
matrix	matrix to be output as string

Returns: std::ostream&

◆ subtraction()

void subtraction	(	const Matrix &	m1,
		const Matrix &	m2,
		Matrix &	tmp
	)

Compute the element-wise subtraction of two matrices and store the result in a third matrix.

This function computes the element-wise subtraction of two input matrices of equal dimensions and stores the result in a third matrix. All input matrices are assumed to be of type Matrix, which internally contains a gsl_matrix pointer for efficiency.

Parameters

m1	The first input matrix.
m2	The second input matrix.
tmp	The output matrix to store the resulting difference.

Precondition: All input matrices must have the same dimensions.

Postcondition: The output matrix 'tmp' contains the element-wise difference of 'm1' and 'm2'.

Classes

Functions

Detailed Description

Function Documentation

◆ a_cross_b_plus_c_cross_d()

◆ addition()

◆ assign_rowVec_to_colVec()

◆ colvec_matrix_multiplication()

◆ const_division()

◆ const_multiplication()

◆ cross()

◆ cross_col()

◆ cross_row()

◆ dot()

◆ dot_col()

◆ dot_row()

◆ get_unit_vector() [1/2]

◆ get_unit_vector() [2/2]

◆ kron() [1/2]

◆ kron() [2/2]

◆ mat_calloc()

◆ multiplication()

◆ negative()

◆ operator*() [1/3]

◆ operator*() [2/3]

◆ operator*() [3/3]

◆ operator*=() [1/2]

◆ operator*=() [2/2]

◆ operator+()

◆ operator+=()

◆ operator-()

◆ operator-=()

◆ operator/()

◆ operator/=()

◆ operator<<()

◆ subtraction()