libMesh::DenseMatrix< T > 模板类 参考

定义用于有限元类型计算的密集矩阵。 用于在求和成全局矩阵之前存储单元刚度矩阵。所有被覆盖的虚函数都记录在dense_matrix_base.h中。 更多...

#include <dof_map.h>

类 libMesh::DenseMatrix< T > 继承关系图:

类
struct	UseBlasLapack
	用于确定是否使用blas_lapack的辅助结构。 更多...

Public 成员函数
	DenseMatrix (const unsigned int new_m=0, const unsigned int new_n=0)
	构造函数。创建一个维度为 m x n 的密集矩阵。 更多...
template<typename T2 >
	DenseMatrix (unsigned int nrow, unsigned int ncol, std::initializer_list< T2 > init_list)
	构造函数。根据行数、列数和初始化列表创建矩阵。 更多...
	DenseMatrix (DenseMatrix &&)=default
	移动构造函数。移动构造函数使用默认实现，因为该类不管理任何内存。 更多...
	DenseMatrix (const DenseMatrix &)=default
	复制构造函数。复制构造函数使用默认实现，因为该类不管理任何内存。 更多...
DenseMatrix &	operator= (const DenseMatrix &)=default
	复制赋值运算符。复制赋值运算符使用默认实现，因为该类不管理任何内存。 更多...
DenseMatrix &	operator= (DenseMatrix &&)=default
	移动赋值运算符。移动赋值运算符使用默认实现，因为该类不管理任何内存。 更多...
virtual	~DenseMatrix ()=default
	析构函数。析构函数使用默认实现，因为该类不管理任何内存。 更多...
virtual void	zero () overridefinal
	将矩阵的所有元素设置为0，并重置任何可能先前设置的分解标志。 这允许例如计算新的LU分解，同时重用相同的存储空间。 更多...
DenseMatrix	sub_matrix (unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
	获取具有最小行和列索引以及子矩阵大小的子矩阵。 更多...
T	operator() (const unsigned int i, const unsigned int j) const
	获取矩阵的 (i,j) 元素。 更多...
T &	operator() (const unsigned int i, const unsigned int j)
	获取矩阵的 (i,j) 元素的可写引用。 更多...
virtual T	el (const unsigned int i, const unsigned int j) const overridefinal
	获取矩阵的 (i,j) 元素。 更多...
virtual T &	el (const unsigned int i, const unsigned int j) overridefinal
	获取矩阵的 (i,j) 元素的可写引用。 更多...
virtual void	left_multiply (const DenseMatrixBase< T > &M2) overridefinal
	左乘以矩阵 M2。 更多...
template<typename T2 >
void	left_multiply (const DenseMatrixBase< T2 > &M2)
	左乘以不同类型的矩阵 M2。 更多...
virtual void	right_multiply (const DenseMatrixBase< T > &M2) overridefinal
	右乘以矩阵 M2。 更多...
template<typename T2 >
void	right_multiply (const DenseMatrixBase< T2 > &M2)
	右乘以不同类型的矩阵 M2。 更多...
void	vector_mult (DenseVector< T > &dest, const DenseVector< T > &arg) const
	执行矩阵-向量乘法，dest := (*this) * arg。 更多...
template<typename T2 >
void	vector_mult (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
	执行矩阵-向量乘法，dest := (*this) * arg ,支持不同类型的矩阵和向量。 更多...
void	vector_mult_transpose (DenseVector< T > &dest, const DenseVector< T > &arg) const
	执行矩阵-向量乘法，dest := (*this)^T * arg。 更多...
template<typename T2 >
void	vector_mult_transpose (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
	执行矩阵-向量乘法，dest := (*this)^T * arg，支持不同类型的矩阵和向量。 更多...
void	vector_mult_add (DenseVector< T > &dest, const T factor, const DenseVector< T > &arg) const
	执行缩放矩阵-向量乘法，结果存储在 dest 中。 dest += factor * (*this) * arg. 更多...
template<typename T2 , typename T3 >
void	vector_mult_add (DenseVector< typename CompareTypes< T, typename CompareTypes< T2, T3 >::supertype >::supertype > &dest, const T2 factor, const DenseVector< T3 > &arg) const
	执行混合类型的缩放矩阵-向量乘法，结果存储在 dest 中。 dest += factor * (*this) * arg. 更多...
void	get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
	将 sub_m x sub_n 的主子矩阵放入 dest 中。 更多...
void	get_principal_submatrix (unsigned int sub_m, DenseMatrix< T > &dest) const
	将 sub_m x sub_m 的主子矩阵放入 dest 中。 更多...
void	outer_product (const DenseVector< T > &a, const DenseVector< T > &b)
	计算两个向量的外积并将结果存储在 (*this) 中。 更多...
template<typename T2 >
DenseMatrix< T > &	operator= (const DenseMatrix< T2 > &other_matrix)
	从另一种矩阵类型复制矩阵。 更多...
void	swap (DenseMatrix< T > &other_matrix)
	STL 风格的交换方法,交换值和分解方法。 更多...
void	resize (const unsigned int new_m, const unsigned int new_n)
	调整矩阵的大小，并调用 zero() 方法。 更多...
void	scale (const T factor)
	将矩阵的每个元素乘以 factor。 更多...
void	scale_column (const unsigned int col, const T factor)
	将矩阵的第 col 列的每个元素乘以 factor。 更多...
DenseMatrix< T > &	operator*= (const T factor)
	将矩阵的每个元素乘以 factor。 更多...
template<typename T2 , typename T3 >
boostcopy::enable_if_c < ScalarTraits< T2 >::value, void >::type	add (const T2 factor, const DenseMatrix< T3 > &mat)
	将 factor 倍的矩阵 mat 添加到此矩阵。 更多...
bool	operator== (const DenseMatrix< T > &mat) const
	检查矩阵是否与另一个矩阵 mat 完全相等。 更多...
bool	operator!= (const DenseMatrix< T > &mat) const
	检查矩阵是否与另一个矩阵 mat 完全不相等。 更多...
DenseMatrix< T > &	operator+= (const DenseMatrix< T > &mat)
	将矩阵 mat 添加到此矩阵。 更多...
DenseMatrix< T > &	operator-= (const DenseMatrix< T > &mat)
	从此矩阵中减去矩阵 mat。 更多...
auto	min () const -> decltype(libmesh_real(T(0)))
	返回矩阵的最小元素或复数情况下的最小实部。 更多...
auto	max () const -> decltype(libmesh_real(T(0)))
	返回矩阵的最大元素或复数情况下的最大实部。 更多...
auto	l1_norm () const -> decltype(std::abs(T(0)))
	返回矩阵的 l1-范数，即最大列和 更多...
auto	linfty_norm () const -> decltype(std::abs(T(0)))
	返回矩阵的 linfty-范数，即最大行和 更多...
void	left_multiply_transpose (const DenseMatrix< T > &A)
	用矩阵 A 的转置左乘。 更多...
template<typename T2 >
void	left_multiply_transpose (const DenseMatrix< T2 > &A)
	用包含不同数值类型的矩阵 A 的转置左乘。 更多...
void	right_multiply_transpose (const DenseMatrix< T > &A)
	用矩阵 A 的转置右乘。 更多...
template<typename T2 >
void	right_multiply_transpose (const DenseMatrix< T2 > &A)
	用包含不同数值类型的矩阵 A 的转置右乘。 更多...
T	transpose (const unsigned int i, const unsigned int j) const
	返回转置矩阵的 (i, j) 元素。 更多...
void	get_transpose (DenseMatrix< T > &dest) const
	将转置矩阵存储在 dest 中。 更多...
std::vector< T > &	get_values ()
	返回对应于存储向量的引用的底层数据存储矢量。 更多...
const std::vector< T > &	get_values () const
	返回底层数据存储矢量的常量引用。 更多...
void	condense (const unsigned int i, const unsigned int j, const T val, DenseVector< T > &rhs)
	将矩阵的 (i,j) 元素凝聚出来，强制其取值为 val。 更多...
void	lu_solve (const DenseVector< T > &b, DenseVector< T > &x)
	解方程组 给定输入向量 b。 更多...
template<typename T2 >
void	cholesky_solve (const DenseVector< T2 > &b, DenseVector< T2 > &x)
	对于对称正定矩阵（SPD），进行Cholesky分解。 分解为 A = L L^T，比标准LU分解大约快两倍。 因此，如果事先知道矩阵是SPD，则可以使用此方法。如果矩阵不是SPD，则会生成错误。 Cholesky分解的一个优点是它不需要主元交换以确保稳定性。 更多...
void	svd (DenseVector< Real > &sigma)
	计算矩阵的奇异值分解（SVD）。 在退出时，sigma包含所有奇异值（按降序排列）。 更多...
void	svd (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
	计算矩阵的“简化”奇异值分解。 在退出时，sigma包含所有奇异值（按降序排列）， U包含左奇异向量，VT包含右奇异向量的转置。 在简化的SVD中，U有min(m,n)列，VT有min(m,n)行。（在“完整”SVD中，U和VT将是方形的。） 更多...
void	svd_solve (const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond=std::numeric_limits< Real >::epsilon()) const
	以最小二乘意义解方程组 。 可能是非方阵和/或秩亏的。 可以通过更改“rcond”参数来控制哪些奇异值被视为零。 对于满足S(i) <= rcond*S(1)的奇异值S(i)，在解的目的上被视为零。 传递负数的rcond将强制使用“机器精度”值。 更多...
void	evd (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag)
	计算一般矩阵的特征值（实部和虚部）。 更多...
void	evd_left (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VL)
	计算一般矩阵 的特征值（实部和虚部）和左特征向量。 更多...
void	evd_right (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VR)
	计算一般矩阵 的特征值（实部和虚部）和右特征向量。 更多...
void	evd_left_and_right (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VL, DenseMatrix< T > &VR)
	计算一般矩阵的特征值（实部和虚部），以及左右特征向量。 更多...
T	det ()
	返回矩阵的行列式。 更多...
unsigned int	m () const
	返回矩阵的行维度。 更多...
unsigned int	n () const
	返回矩阵的列维度。 更多...
void	print (std::ostream &os=libMesh::out) const
	漂亮地打印矩阵，默认为 libMesh::out。 更多...
void	print_scientific (std::ostream &os, unsigned precision=8) const
	以科学计数法格式打印矩阵。 更多...
template<typename T2 , typename T3 >
boostcopy::enable_if_c < ScalarTraits< T2 >::value, void >::type	add (const T2 factor, const DenseMatrixBase< T3 > &mat)
	将 factor 添加到矩阵的每个元素。 这应该仅在 T += T2 * T3 是有效的 C++，且 T2 是标量的情况下工作。返回类型是 void。 更多...
DenseVector< T >	diagonal () const
	返回矩阵的对角线。 更多...

Public 属性
bool	use_blas_lapack
	计算密矩阵的逆（假设可逆性）， 首先计算LU分解，然后执行多次回代步骤。 遵循在Web上提供的Numerical Recipes in C中的算法。 更多...

Protected 成员函数
void	condense (const unsigned int i, const unsigned int j, const T val, DenseVectorBase< T > &rhs)
	将矩阵的 (i,j) 条目压缩出来，强制它取值为 val。这对于在数值模拟中应用边界条件很有用。 保留矩阵的对称性。 更多...

静态 Protected 成员函数
static void	multiply (DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
	辅助函数 - 执行计算 M1 = M2 * M3 其中: M1 = (m x n) M2 = (m x p) M3 = (p x n) 更多...

Protected 属性
unsigned int	_m
	行维度。 更多...
unsigned int	_n
	列维度。 更多...

Private 类型
enum	DecompositionType { LU =0, CHOLESKY =1, LU_BLAS_LAPACK, NONE }
	_cholesky_back_substitute 的分解方案会改变矩阵A的条目。因此，在调用A.lu_solve()后再次调用A.cholesky_solve()是错误的， 因为结果可能不符合任何预期结果。这个typedef跟踪在此矩阵上执行了哪种分解。 更多...
enum	_BLAS_Multiply_Flag { LEFT_MULTIPLY = 0, RIGHT_MULTIPLY, LEFT_MULTIPLY_TRANSPOSE, RIGHT_MULTIPLY_TRANSPOSE }
	用于确定_multiply_blas函数行为的枚举。 更多...
typedef PetscBLASInt	pivot_index_t
	用于存储枢轴索引的数组。可能被当前活动的任何因式分解使用， 类的客户端不应出于任何原因依赖它。 更多...
typedef int	pivot_index_t

Private 成员函数
void	_lu_decompose ()
	形成矩阵的LU分解。此函数是私有的，因为它仅在lu_solve(...)函数的实现部分中被调用。 更多...
void	_lu_back_substitute (const DenseVector< T > &b, DenseVector< T > &x) const
	通过回代步骤解方程Ax=b。此函数是私有的，因为它仅在lu_solve(...)函数的实现部分中被调用。 更多...
void	_cholesky_decompose ()
	将对称正定矩阵分解为两个下三角矩阵的乘积，即A = LL^T。 更多...
template<typename T2 >
void	_cholesky_back_substitute (const DenseVector< T2 > &b, DenseVector< T2 > &x) const
	根据矩阵A的Cholesky分解解方程Ax=b，得到未知值x和rhs b。 更多...
void	_multiply_blas (const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
	_multiply_blas函数使用BLAS gemm函数计算A <- op(A) * op(B)。 用于right_multiply()、left_multiply()、right_multiply_transpose()和 left_multiply_transpose()例程。 [ 在dense_matrix_blas_lapack.C中实现 ] 更多...
void	_lu_decompose_lapack ()
	使用Lapack例程“getrf”计算矩阵的LU分解。此例程只能由lu_solve()函数的“use_blas_lapack”分支使用。 调用此函数后，矩阵将被其因式分解版本替换，并且DecompositionType将设置为LU_BLAS_LAPACK。 [ 在dense_matrix_blas_lapack.C中实现 ] 更多...
void	_svd_lapack (DenseVector< Real > &sigma)
	使用Lapack例程“getsvd”计算矩阵的奇异值分解。 [ 在dense_matrix_blas_lapack.C中实现 ] 更多...
void	_svd_lapack (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
	使用Lapack例程“getsvd”计算矩阵的“缩减”奇异值分解。 [ 在dense_matrix_blas_lapack.C中实现 ] 更多...
void	_svd_solve_lapack (const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond) const
	由svd_solve(rhs)调用。 更多...
void	_svd_helper (char JOBU, char JOBVT, std::vector< Real > &sigma_val, std::vector< Number > &U_val, std::vector< Number > &VT_val)
	实际执行SVD的辅助函数。 [ 在dense_matrix_blas_lapack.C中实现 ] 更多...
void	_evd_lapack (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=nullptr, DenseMatrix< T > *VR=nullptr)
	使用Lapack例程“DGEEV”计算矩阵的特征值。如果VR和/或VL不为nullptr， 则此函数还计算并返回右和/或左特征向量的矩阵。 [ 在dense_matrix_blas_lapack.C中实现 ] 更多...
void	_lu_back_substitute_lapack (const DenseVector< T > &b, DenseVector< T > &x)
	与_lu_decompose_lapack()相对应的伴随函数。不要直接使用，通过公共lu_solve()接口调用。 由于我们只是调用LAPACK例程，该函数在逻辑上是const的，因为它不修改矩阵， 但由于它比const_cast更少的麻烦，所以直接声明函数为非const。 [ 在dense_matrix_blas_lapack.C中实现 ] 更多...
void	_matvec_blas (T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
	使用BLAS GEMV函数（通过PETSc）计算 更多...
template<typename T2 >
void	_left_multiply_transpose (const DenseMatrix< T2 > &A)
	由矩阵 A 的转置左乘，该矩阵可能包含不同的数值类型。 更多...
template<typename T2 >
void	_right_multiply_transpose (const DenseMatrix< T2 > &A)
	由矩阵 A 的转置右乘，该矩阵可能包含不同的数值类型。 更多...

Private 属性
std::vector< T >	_val
	实际的数据值，存储为1D数组。 更多...
DecompositionType	_decomposition_type
	此标志跟踪在矩阵上执行的分解类型。 更多...
std::vector< pivot_index_t >	_pivots

详细描述

template<typename T>
class libMesh::DenseMatrix< T >

定义用于有限元类型计算的密集矩阵。 用于在求和成全局矩阵之前存储单元刚度矩阵。所有被覆盖的虚函数都记录在dense_matrix_base.h中。

作者

Benjamin S. Kirk

日期

2002 用于有限元组装和数值计算的密集矩阵对象。

在文件 dof_map.h 第 65 行定义.

成员类型定义说明

template<typename T>

typedef PetscBLASInt libMesh::DenseMatrix< T >::pivot_index_t

private

用于存储枢轴索引的数组。可能被当前活动的任何因式分解使用， 类的客户端不应出于任何原因依赖它。

在文件 dense_matrix.h 第 849 行定义.

template<typename T>

typedef int libMesh::DenseMatrix< T >::pivot_index_t

private

在文件 dense_matrix.h 第 851 行定义.

成员枚举类型说明

template<typename T>

enum libMesh::DenseMatrix::_BLAS_Multiply_Flag

private

用于确定_multiply_blas函数行为的枚举。

枚举值
LEFT_MULTIPLY
RIGHT_MULTIPLY
LEFT_MULTIPLY_TRANSPOSE
RIGHT_MULTIPLY_TRANSPOSE

在文件 dense_matrix.h 第 756 行定义.

                           {
    LEFT_MULTIPLY = 0,
    RIGHT_MULTIPLY,
    LEFT_MULTIPLY_TRANSPOSE,
    RIGHT_MULTIPLY_TRANSPOSE
  };

template<typename T>

enum libMesh::DenseMatrix::DecompositionType

private

_cholesky_back_substitute 的分解方案会改变矩阵A的条目。因此，在调用A.lu_solve()后再次调用A.cholesky_solve()是错误的， 因为结果可能不符合任何预期结果。这个typedef跟踪在此矩阵上执行了哪种分解。

枚举值
LU
CHOLESKY
LU_BLAS_LAPACK
NONE

在文件 dense_matrix.h 第 746 行定义.

{LU=0, CHOLESKY=1, LU_BLAS_LAPACK, NONE};

构造及析构函数说明

template<typename T >

libMesh::DenseMatrix< T >::DenseMatrix ( const unsigned int  new_m = 0, const unsigned int  new_n = 0  )

inline

构造函数。创建一个维度为 m x n 的密集矩阵。

参数

new_m	矩阵的行数。
new_n	矩阵的列数。

在文件 dense_matrix.h 第 954 行定义.

参考 libMesh::DenseMatrix< T >::resize().

                                                      :
  DenseMatrixBase<T>(new_m,new_n),
  use_blas_lapack(DenseMatrix<T>::UseBlasLapack::value),
  _val(),
  _decomposition_type(NONE)
{
  this->resize(new_m,new_n);
}

template<typename T >

template<typename T2 >

libMesh::DenseMatrix< T >::DenseMatrix	(	unsigned int	nrow,
		unsigned int	ncol,
		std::initializer_list< T2 >	init_list
	)

构造函数。根据行数、列数和初始化列表创建矩阵。

参数

new_m	矩阵的行数。
new_n	矩阵的列数。
init_list	包含行主元素值的初始化列表，其长度必须为 nrow * ncol。

在文件 dense_matrix.h 第 966 行定义.

                                                               :
  DenseMatrixBase<T>(nrow, ncol),
  use_blas_lapack(DenseMatrix<T>::UseBlasLapack::value),
  _val(init_list.begin(), init_list.end()),
  _decomposition_type(NONE)
{
 // 确保传递的数据量与矩阵大小一致。
  libmesh_assert_equal_to(nrow * ncol, init_list.size());
}

template<typename T>

libMesh::DenseMatrix< T >::DenseMatrix ( DenseMatrix< T > &&  )

default

移动构造函数。移动构造函数使用默认实现，因为该类不管理任何内存。

template<typename T>

libMesh::DenseMatrix< T >::DenseMatrix ( const DenseMatrix< T > &  )

default

复制构造函数。复制构造函数使用默认实现，因为该类不管理任何内存。

template<typename T>

virtual libMesh::DenseMatrix< T >::~DenseMatrix ( )

virtualdefault

析构函数。析构函数使用默认实现，因为该类不管理任何内存。

成员函数说明

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::_cholesky_back_substitute ( const DenseVector< T2 > &  b, DenseVector< T2 > &  x  ) const

private

根据矩阵A的Cholesky分解解方程Ax=b，得到未知值x和rhs b。

使用 Cholesky 分解的回代步骤，求解 Ly=b 和 L^T x = y。

注解

当A是实值，b和x是复值时，可以使用此方法。

模板参数

T2	解方程时使用的数据类型。

参数

b	右侧向量。
x	存储解的向量。

模板参数

T2	输入向量 b 和解向量 x 的数据类型。

参数

b	输入向量 b。
x	存储解的目标向量 x。

在文件 dense_matrix_impl.h 第 885 行定义.

参考 libMesh::DenseVector< T >::resize().

{
  // 确保矩阵是方阵
  libmesh_assert_equal_to(this->m(), this->n());

  // 行和列的数量的简写。
  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();

  // 方便引用 *this
  const DenseMatrix<T> &A = *this;

  // 计算 Ax = b 的解。
  x.resize(n_rows);

  // 解 Ly=b
  for (unsigned int i = 0; i < n_cols; ++i)
  {
    T2 temp = b(i);

    for (unsigned int k = 0; k < i; ++k)
      temp -= A(i, k) * x(k);

    x(i) = temp / A(i, i);
  }

  // 解 L^T x = y
  for (unsigned int i = 0; i < n_cols; ++i)
  {
    const unsigned int ib = (n_cols - 1) - i;

    for (unsigned int k = (ib + 1); k < n_cols; ++k)
      x(ib) -= A(k, ib) * x(k);

    x(ib) /= A(ib, ib);
  }
}

template<typename T >

void libMesh::DenseMatrix< T >::_cholesky_decompose ( )

private

将对称正定矩阵分解为两个下三角矩阵的乘积，即A = LL^T。

使用 Cholesky 分解的实际分解过程。

注解

如果矩阵不是SPD，则此程序会生成错误。

该算法基于 Numerical Recipes in C 书中的 Cholesky 分解。

在文件 dense_matrix_impl.h 第 835 行定义.

参考 std::sqrt().

{
  // 如果调用了此函数，则矩阵不应该有任何先前的分解。
  libmesh_assert_equal_to(this->_decomposition_type, NONE);

  // 为了确保是方阵，检查行和列的数量是否相同
  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();

  // 仅为了确保是方阵
  libmesh_assert_equal_to(n_rows, n_cols);

  // 方便引用 *this
  DenseMatrix<T> &A = *this;

  for (unsigned int i = 0; i < n_rows; ++i)
  {
    for (unsigned int j = i; j < n_cols; ++j)
    {
      for (unsigned int k = 0; k < i; ++k)
        A(i, j) -= A(i, k) * A(j, k);

      if (i == j)
      {
#ifndef LIBMESH_USE_COMPLEX_NUMBERS
        libmesh_error_msg_if(A(i, j) <= 0.0,
                             "Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
#endif

        A(i, i) = std::sqrt(A(i, j));
      }
      else
        A(j, i) = A(i, j) / A(i, i);
    }
  }

  // 设置 CHOLESKY 分解标志
  this->_decomposition_type = CHOLESKY;
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_evd_lapack ( DenseVector< T > &  lambda_real, DenseVector< T > &  lambda_imag, DenseMatrix< T > *  VL = nullptr, DenseMatrix< T > *  VR = nullptr  )

private

使用Lapack例程“DGEEV”计算矩阵的特征值。如果VR和/或VL不为nullptr， 则此函数还计算并返回右和/或左特征向量的矩阵。 [ 在dense_matrix_blas_lapack.C中实现 ]

参数

lambda_real	存储实部特征值的向量。
lambda_imag	存储虚部特征值的向量。
VL	存储左特征向量的矩阵。
VR	存储右特征向量的矩阵。

在文件 dense_matrix_blas_lapack.C 第 695 行定义.

参考 libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libMesh::pPS(), libMesh::DenseVector< T >::resize() , 以及 libMesh::DenseMatrix< T >::resize().

{
  // This algorithm only works for square matrices, so verify this up front.
  libmesh_error_msg_if(this->m() != this->n(), "Can only compute eigen-decompositions for square matrices.");

  // If the user requests left or right eigenvectors, we have to make
  // sure and pass the transpose of this matrix to Lapack, otherwise
  // it will compute the inverse transpose of what we are
  // after... since we know the matrix is square, we can just swap
  // entries in place.  If the user does not request eigenvectors, we
  // can skip this extra step, since the eigenvalues for A and A^T are
  // the same.
  if (VL || VR)
    {
      for (auto i : make_range(this->_m))
        for (auto j : make_range(i))
          std::swap((*this)(i,j), (*this)(j,i));
    }

  // The calling sequence for dgeev is:
  // DGEEV (character  JOBVL,
  //        character  JOBVR,
  //        integer N,
  //        double precision, dimension( lda, * )  A,
  //        integer  LDA,
  //        double precision, dimension( * )  WR,
  //        double precision, dimension( * )  WI,
  //        double precision, dimension( ldvl, * )  VL,
  //        integer  LDVL,
  //        double precision, dimension( ldvr, * )  VR,
  //        integer  LDVR,
  //        double precision, dimension( * )  WORK,
  //        integer  LWORK,
  //        integer  INFO)

  // JOBVL (input)
  //       = 'N': left eigenvectors of A are not computed;
  //       = 'V': left eigenvectors of A are computed.
  char JOBVL = VL ? 'V' : 'N';

  // JOBVR (input)
  //       = 'N': right eigenvectors of A are not computed;
  //       = 'V': right eigenvectors of A are computed.
  char JOBVR = VR ? 'V' : 'N';

  // N (input)
  //   The number of rows/cols of the matrix A.  N >= 0.
  PetscBLASInt N = this->m();

  // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  //   On entry, the N-by-N matrix A.
  //   On exit, A has been overwritten.

  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,N).
  PetscBLASInt LDA = N;

  // WR (output) double precision array, dimension (N)
  // WI (output) double precision array, dimension (N)
  //    WR and WI contain the real and imaginary parts,
  //    respectively, of the computed eigenvalues.  Complex
  //    conjugate pairs of eigenvalues appear consecutively
  //    with the eigenvalue having the positive imaginary part
  //    first.
  lambda_real.resize(N);
  lambda_imag.resize(N);

  // VL (output) double precision array, dimension (LDVL,N)
  //    If JOBVL = 'V', the left eigenvectors u(j) are stored one
  //    after another in the columns of VL, in the same order
  //    as their eigenvalues.
  //    If JOBVL = 'N', VL is not referenced.
  //    If the j-th eigenvalue is real, then u(j) = VL(:,j),
  //    the j-th column of VL.
  //    If the j-th and (j+1)-st eigenvalues form a complex
  //    conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
  //    u(j+1) = VL(:,j) - i*VL(:,j+1).
  // Will be set below if needed.

  // LDVL (input)
  //      The leading dimension of the array VL.  LDVL >= 1; if
  //      JOBVL = 'V', LDVL >= N.
  PetscBLASInt LDVL = VL ? N : 1;

  // VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
  //    If JOBVR = 'V', the right eigenvectors v(j) are stored one
  //    after another in the columns of VR, in the same order
  //    as their eigenvalues.
  //    If JOBVR = 'N', VR is not referenced.
  //    If the j-th eigenvalue is real, then v(j) = VR(:,j),
  //    the j-th column of VR.
  //    If the j-th and (j+1)-st eigenvalues form a complex
  //    conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
  //    v(j+1) = VR(:,j) - i*VR(:,j+1).
  // Will be set below if needed.

  // LDVR (input)
  //      The leading dimension of the array VR.  LDVR >= 1; if
  //      JOBVR = 'V', LDVR >= N.
  PetscBLASInt LDVR = VR ? N : 1;

  // WORK (workspace/output) double precision array, dimension (MAX(1,LWORK))
  //      On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  //
  // LWORK (input)
  //       The dimension of the array WORK.  LWORK >= max(1,3*N), and
  //       if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
  //       performance, LWORK must generally be larger.
  //
  //       If LWORK = -1, then a workspace query is assumed; the routine
  //       only calculates the optimal size of the WORK array, returns
  //       this value as the first entry of the WORK array, and no error
  //       message related to LWORK is issued by XERBLA.
  PetscBLASInt LWORK = (VR || VL) ? 4*N : 3*N;
  std::vector<T> WORK(LWORK);

  // INFO (output)
  //      = 0:  successful exit
  //      < 0:  if INFO = -i, the i-th argument had an illegal value.
  //      > 0:  if INFO = i, the QR algorithm failed to compute all the
  //            eigenvalues, and no eigenvectors or condition numbers
  //            have been computed; elements 1:ILO-1 and i+1:N of WR
  //            and WI contain eigenvalues which have converged.
  PetscBLASInt INFO = 0;

  // Get references to raw data
  std::vector<T> & lambda_real_val = lambda_real.get_values();
  std::vector<T> & lambda_imag_val = lambda_imag.get_values();

  // Set up eigenvector storage if necessary.
  T * VR_ptr = nullptr;
  if (VR)
    {
      VR->resize(N, N);
      VR_ptr = VR->get_values().data();
    }

  T * VL_ptr = nullptr;
  if (VL)
    {
      VL->resize(N, N);
      VL_ptr = VL->get_values().data();
    }

  // Ready to call the Lapack routine through PETSc's interface
  LAPACKgeev_(&JOBVL,
              &JOBVR,
              &N,
              pPS(_val.data()),
              &LDA,
              pPS(lambda_real_val.data()),
              pPS(lambda_imag_val.data()),
              pPS(VL_ptr),
              &LDVL,
              pPS(VR_ptr),
              &LDVR,
              pPS(WORK.data()),
              &LWORK,
              &INFO);

  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack eigenvalue calculation!");

  // If the user requested either right or left eigenvectors, LAPACK
  // has now computed the transpose of the desired matrix, i.e. V^T
  // instead of V.  We could leave this up to user code to handle, but
  // rather than risking getting very unexpected results, we'll just
  // transpose it in place before handing it back.
  if (VR)
    {
      for (auto i : make_range(N))
        for (auto j : make_range(i))
          std::swap((*VR)(i,j), (*VR)(j,i));
    }

  if (VL)
    {
      for (auto i : make_range(N))
        for (auto j : make_range(i))
          std::swap((*VL)(i,j), (*VL)(j,i));
    }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::_left_multiply_transpose ( const DenseMatrix< T2 > &  A )

private

由矩阵 A 的转置左乘，该矩阵可能包含不同的数值类型。

执行左乘转置矩阵的实际操作。

模板参数

T2	矩阵元素的新数据类型。

参数

A

另一个矩阵。

如果正在执行 (A^T)*A，将执行一种优化方法；否则，将执行通用的左乘转置操作。

参数

A

要左乘的转置矩阵。

在文件 dense_matrix_impl.h 第 125 行定义.

参考 libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n() , 以及 libMesh::DenseMatrix< T >::transpose().

{
  // 检查是否正在执行 (A^T)*A
  if (this == &A)
  {
    // 复制 *this 的副本
    DenseMatrix<T> B(*this);

    // 简单但低效的方式
    // 返回 this->left_multiply_transpose(B);

    // 更有效，但代码更多的方式
    // 如果 A 是 mxn，则结果将是大小为 n x n 的方阵。
    const unsigned int n_rows = A.m();
    const unsigned int n_cols = A.n();

    // 调整 *this 的大小并将所有条目清零。
    this->resize(n_cols, n_cols);

    // 计算下三角部分
    for (unsigned int i = 0; i < n_cols; ++i)
      for (unsigned int j = 0; j <= i; ++j)
        for (unsigned int k = 0; k < n_rows; ++k) // 内积在 n_rows 上
          (*this)(i, j) += B(k, i) * B(k, j);

    // 将下三角部分复制到上三角部分
    for (unsigned int i = 0; i < n_cols; ++i)
      for (unsigned int j = i + 1; j < n_cols; ++j)
        (*this)(i, j) = (*this)(j, i);
  }
  else
  {
    DenseMatrix<T> B(*this);

    this->resize(A.n(), B.n());

    libmesh_assert_equal_to(A.m(), B.m());
    libmesh_assert_equal_to(this->m(), A.n());
    libmesh_assert_equal_to(this->n(), B.n());

    const unsigned int m_s = A.n();
    const unsigned int p_s = A.m();
    const unsigned int n_s = this->n();

    // 以这种方式执行是因为
    // 很大的机会（至少对于约束矩阵而言）
    // A.transpose(i,k) = 0。
    for (unsigned int i = 0; i < m_s; i++)
      for (unsigned int k = 0; k < p_s; k++)
        if (A.transpose(i, k) != 0.)
          for (unsigned int j = 0; j < n_s; j++)
            (*this)(i, j) += A.transpose(i, k) * B(k, j);
  }
}

template<typename T >

void libMesh::DenseMatrix< T >::_lu_back_substitute ( const DenseVector< T > &  b, DenseVector< T > &  x  ) const

private

通过回代步骤解方程Ax=b。此函数是私有的，因为它仅在lu_solve(...)函数的实现部分中被调用。

使用 LU 分解的反向代入解线性系统。

参数

b	输入向量 b。
x	存储结果的目标向量 x。

在文件 dense_matrix_impl.h 第 656 行定义.

参考 libMesh::DenseVector< T >::resize() , 以及 libMesh::DenseVector< T >::size().

{
  const unsigned int n_cols = this->n();

  libmesh_assert_equal_to(this->m(), n_cols);
  libmesh_assert_equal_to(this->m(), b.size());

  x.resize(n_cols);

  // 方便引用 *this
  const DenseMatrix<T> & A = *this;

  // 临时向量存储。使用这个向量而不是修改 RHS。
  DenseVector<T> z = b;

  // 下三角 "从上到下" 解法步骤，考虑了主元素
  for (unsigned int i = 0; i < n_cols; ++i)
  {
    // 交换
    if (_pivots[i] != static_cast<pivot_index_t>(i))
      std::swap(z(i), z(_pivots[i]));

    x(i) = z(i);

    for (unsigned int j = 0; j < i; ++j)
      x(i) -= A(i, j) * x(j);

    x(i) /= A(i, i);
  }

  // 上三角 "从下到上" 解法步骤
  const unsigned int last_row = n_cols - 1;

  for (int i = last_row; i >= 0; --i)
  {
    for (int j = i + 1; j < static_cast<int>(n_cols); ++j)
      x(i) -= A(i, j) * x(j);
  }
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_lu_back_substitute_lapack ( const DenseVector< T > &  b, DenseVector< T > &  x  )

private

与_lu_decompose_lapack()相对应的伴随函数。不要直接使用，通过公共lu_solve()接口调用。 由于我们只是调用LAPACK例程，该函数在逻辑上是const的，因为它不修改矩阵， 但由于它比const_cast更少的麻烦，所以直接声明函数为非const。 [ 在dense_matrix_blas_lapack.C中实现 ]

参数

b	右侧向量。
x	存储解的向量。

在文件 dense_matrix_blas_lapack.C 第 901 行定义.

参考 libMesh::DenseVector< T >::get_values() , 以及 libMesh::pPS().

{
  // The calling sequence for getrs is:
  // dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)

  // trans (input)
  //       'n' for no transpose, 't' for transpose
  char TRANS[] = "t";

  // N (input)
  //   The order of the matrix A.  N >= 0.
  PetscBLASInt N = this->m();


  // NRHS (input)
  //      The number of right hand sides, i.e., the number of columns
  //      of the matrix B.  NRHS >= 0.
  PetscBLASInt NRHS = 1;

  // A (input) double precision array, dimension (LDA,N)
  //   The factors L and U from the factorization A = P*L*U
  //   as computed by dgetrf.

  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,N).
  PetscBLASInt LDA = N;

  // ipiv (input) int array, dimension (N)
  //      The pivot indices from DGETRF; for 1<=i<=N, row i of the
  //      matrix was interchanged with row IPIV(i).
  // Here, we pass _pivots.data() which was computed in _lu_decompose_lapack

  // B (input/output) double precision array, dimension (LDB,NRHS)
  //   On entry, the right hand side matrix B.
  //   On exit, the solution matrix X.
  // Here, we pass a copy of the rhs vector's data array in x, so that the
  // passed right-hand side b is unmodified.  I don't see a way around this
  // copy if we want to maintain an unmodified rhs in LibMesh.
  x = b;
  std::vector<T> & x_vec = x.get_values();

  // We can avoid the copy if we don't care about overwriting the RHS: just
  // pass b to the Lapack routine and then swap with x before exiting
  // std::vector<T> & x_vec = b.get_values();

  // LDB (input)
  //     The leading dimension of the array B.  LDB >= max(1,N).
  PetscBLASInt LDB = N;

  // INFO (output)
  //      = 0:  successful exit
  //      < 0:  if INFO = -i, the i-th argument had an illegal value
  PetscBLASInt INFO = 0;

  // Finally, ready to call the Lapack getrs function
  LAPACKgetrs_(TRANS, &N, &NRHS, pPS(_val.data()), &LDA,
               _pivots.data(), pPS(x_vec.data()), &LDB, &INFO);

  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack LU solve!");

  // Don't do this if you already made a copy of b above
  // Swap b and x.  The solution will then be in x, and whatever was originally
  // in x, maybe garbage, maybe nothing, will be in b.
  // FIXME: Rewrite the LU and Cholesky solves to just take one input, and overwrite
  // the input.  This *should* make user code simpler, as they don't have to create
  // an extra vector just to pass it in to the solve function!
  // b.swap(x);
}

template<typename T >

void libMesh::DenseMatrix< T >::_lu_decompose ( )

private

形成矩阵的LU分解。此函数是私有的，因为它仅在lu_solve(...)函数的实现部分中被调用。

使用 LU 分解的实际分解过程。

在文件 dense_matrix_impl.h 第 700 行定义.

参考 std::abs(), libMesh::DenseMatrix< T >::resize() , 以及 libMesh::zero.

{
  // 如果调用了此函数，则矩阵不应该有任何先前的分解。
  libmesh_assert_equal_to(this->_decomposition_type, NONE);

  // 获取矩阵大小并确保其是方阵
  const unsigned int n_rows = this->m();

  // 方便引用 *this
  DenseMatrix<T> & A = *this;

  _pivots.resize(n_rows);

  for (unsigned int i = 0; i < n_rows; ++i)
  {
    // 通过向下搜索第 i 列找到主元素所在的行
    _pivots[i] = i;

    // std::abs(complex) 必须返回一个实数!
    auto the_max = std::abs(A(i, i));
    for (unsigned int j = i + 1; j < n_rows; ++j)
    {
      auto candidate_max = std::abs(A(j, i));
      if (the_max < candidate_max)
      {
        the_max = candidate_max;
        _pivots[i] = j;
      }
    }

    // 如果最大值在不同的行中找到，则交换行
    // 这里交换整个行，在高斯消元中，只交换 A(i,j) 和 A(p(i),j) 的子行，对于 j>i
    if (_pivots[i] != static_cast<pivot_index_t>(i))
    {
      for (unsigned int j = 0; j < n_rows; ++j)
        std::swap(A(i, j), A(_pivots[i], j));
    }

    // 如果找到的最大绝对值条目为零，则矩阵是奇异的
    libmesh_error_msg_if(A(i, i) == libMesh::zero, "Matrix A is singular!");

    // 将第 i 行的上三角条目缩放为对角线条目
    // 注意：不要缩放对角线条目本身！
    const T diag_inv = 1. / A(i, i);
    for (unsigned int j = i + 1; j < n_rows; ++j)
      A(i, j) *= diag_inv;

    // 通过减去（对角线已缩放的）第 i 行的上三角部分，更新剩余子矩阵 A[i+1:m][i+1:m]
    // 通过每行的第 i 列条目进行缩放。在行操作方面，这是：
    // 对于每个 r > i
    //   SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
    //
    // 如果我们也缩放了第 i 列，就像在高斯消元中一样，这将“零”第 i 列的条目。
    for (unsigned int row = i + 1; row < n_rows; ++row)
      for (unsigned int col = i + 1; col < n_rows; ++col)
        A(row, col) -= A(row, i) * A(i, col);
  } // end i 循环

  // 设置 LU 分解标志
  this->_decomposition_type = LU;
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_lu_decompose_lapack ( )

private

使用Lapack例程“getrf”计算矩阵的LU分解。此例程只能由lu_solve()函数的“use_blas_lapack”分支使用。 调用此函数后，矩阵将被其因式分解版本替换，并且DecompositionType将设置为LU_BLAS_LAPACK。 [ 在dense_matrix_blas_lapack.C中实现 ]

在文件 dense_matrix_blas_lapack.C 第 210 行定义.

参考 libMesh::pPS().

{
  // If this function was called, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert_equal_to (this->_decomposition_type, NONE);

  // The calling sequence for dgetrf is:
  // dgetrf(M, N, A, lda, ipiv, info)

  // M (input)
  //   The number of rows of the matrix A.  M >= 0.
  // In C/C++, pass the number of *cols* of A
  PetscBLASInt M = this->n();

  // N (input)
  //   The number of columns of the matrix A.  N >= 0.
  // In C/C++, pass the number of *rows* of A
  PetscBLASInt N = this->m();

  // A (input/output) double precision array, dimension (LDA,N)
  //   On entry, the M-by-N matrix to be factored.
  //   On exit, the factors L and U from the factorization
  //   A = P*L*U; the unit diagonal elements of L are not stored.

  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,M).
  PetscBLASInt LDA = M;

  // ipiv (output) integer array, dimension (min(m,n))
  //      The pivot indices; for 1 <= i <= min(m,n), row i of the
  //      matrix was interchanged with row IPIV(i).
  // Here, we pass _pivots.data(), a private class member used to store pivots
  this->_pivots.resize( std::min(M,N) );

  // info (output)
  //      = 0:  successful exit
  //      < 0:  if INFO = -i, the i-th argument had an illegal value
  //      > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
  //            has been completed, but the factor U is exactly
  //            singular, and division by zero will occur if it is used
  //            to solve a system of equations.
  PetscBLASInt INFO = 0;

  // Ready to call the actual factorization routine through PETSc's interface
  LAPACKgetrf_(&M, &N, pPS(this->_val.data()), &LDA, _pivots.data(),
               &INFO);

  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack LU factorization!");

  // Set the flag for LU decomposition
  this->_decomposition_type = LU_BLAS_LAPACK;
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_matvec_blas ( T  alpha, T  beta, DenseVector< T > &  dest, const DenseVector< T > &  arg, bool  trans = false  ) const

private

使用BLAS GEMV函数（通过PETSc）计算

dest := alpha*A*arg + beta*dest

其中alpha和beta是标量，A是此矩阵，arg和dest是适当大小的输入向量。 如果trans为true，则计算转置matvec。默认情况下，trans==false。

[ 在dense_matrix_blas_lapack.C 中实现 ]

参数

alpha	标量alpha。
beta	标量beta。
dest	存储结果的向量。
arg	输入向量。
trans	如果为true，则计算转置matvec。

在文件 dense_matrix_blas_lapack.C 第 990 行定义.

参考 libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libMesh::pPS() , 以及 libMesh::DenseVector< T >::size().

{
  // Ensure that dest and arg sizes are compatible
  if (!trans)
    {
      // dest  ~ A     * arg
      // (mx1)   (mxn) * (nx1)
      libmesh_error_msg_if((dest.size() != this->m()) || (arg.size() != this->n()),
                           "Improper input argument sizes!");
    }

  else // trans == true
    {
      // Ensure that dest and arg are proper size
      // dest  ~ A^T   * arg
      // (nx1)   (nxm) * (mx1)
      libmesh_error_msg_if((dest.size() != this->n()) || (arg.size() != this->m()),
                           "Improper input argument sizes!");
    }

  // Calling sequence for dgemv:
  //
  // dgemv(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)

  // TRANS (input)
  //       't' for transpose, 'n' for non-transpose multiply
  // We store everything in row-major order, so pass the transpose flag for
  // non-transposed matvecs and the 'n' flag for transposed matvecs
  char TRANS[] = "t";
  if (trans)
    TRANS[0] = 'n';

  // M (input)
  //   On entry, M specifies the number of rows of the matrix A.
  // In C/C++, pass the number of *cols* of A
  PetscBLASInt M = this->n();

  // N (input)
  //   On entry, N specifies the number of columns of the matrix A.
  // In C/C++, pass the number of *rows* of A
  PetscBLASInt N = this->m();

  // ALPHA (input)
  // The scalar constant passed to this function

  // A (input) double precision array of DIMENSION ( LDA, n ).
  //   Before entry, the leading m by n part of the array A must
  //   contain the matrix of coefficients.
  // The matrix, *this.  Note that _matvec_blas is called from
  // a const function, vector_mult(), and so we have made this function const
  // as well.  Since BLAS knows nothing about const, we have to cast it away
  // now.
  DenseMatrix<T> & a_ref = const_cast<DenseMatrix<T> &> ( *this );
  std::vector<T> & a = a_ref.get_values();

  // LDA (input)
  //     On entry, LDA specifies the first dimension of A as declared
  //     in the calling (sub) program. LDA must be at least
  //     max( 1, m ).
  PetscBLASInt LDA = M;

  // X (input) double precision array of DIMENSION at least
  //   ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  //   and at least
  //   ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  //   Before entry, the incremented array X must contain the
  //   vector x.
  // Here, we must cast away the const-ness of "arg" since BLAS knows
  // nothing about const
  DenseVector<T> & x_ref = const_cast<DenseVector<T> &> ( arg );
  std::vector<T> & x = x_ref.get_values();

  // INCX (input)
  //      On entry, INCX specifies the increment for the elements of
  //      X. INCX must not be zero.
  PetscBLASInt INCX = 1;

  // BETA (input)
  //      On entry, BETA specifies the scalar beta. When BETA is
  //      supplied as zero then Y need not be set on input.
  // The second scalar constant passed to this function

  // Y (input) double precision array of DIMENSION at least
  //   ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  //   and at least
  //   ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  //   Before entry with BETA non-zero, the incremented array Y
  //   must contain the vector y. On exit, Y is overwritten by the
  //   updated vector y.
  // The input vector "dest"
  std::vector<T> & y = dest.get_values();

  // INCY (input)
  //      On entry, INCY specifies the increment for the elements of
  //      Y. INCY must not be zero.
  PetscBLASInt INCY = 1;

  // Finally, ready to call the BLAS function
  BLASgemv_(TRANS, &M, &N, pPS(&alpha), pPS(a.data()), &LDA,
            pPS(x.data()), &INCX, pPS(&beta), pPS(y.data()), &INCY);
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_multiply_blas ( const DenseMatrixBase< T > &  other, _BLAS_Multiply_Flag  flag  )

private

_multiply_blas函数使用BLAS gemm函数计算A <- op(A) * op(B)。 用于right_multiply()、left_multiply()、right_multiply_transpose()和 left_multiply_transpose()例程。 [ 在dense_matrix_blas_lapack.C中实现 ]

参数

other	另一个矩阵。
flag	用于指定操作类型的标志。

在文件 dense_matrix_blas_lapack.C 第 37 行定义.

参考 libMesh::DenseMatrix< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n() , 以及 libMesh::pPS().

{
  int result_size = 0;

  // For each case, determine the size of the final result make sure
  // that the inner dimensions match
  switch (flag)
    {
    case LEFT_MULTIPLY:
      {
        result_size = other.m() * this->n();
        if (other.n() == this->m())
          break;
      }
      libmesh_fallthrough();
    case RIGHT_MULTIPLY:
      {
        result_size = other.n() * this->m();
        if (other.m() == this->n())
          break;
      }
      libmesh_fallthrough();
    case LEFT_MULTIPLY_TRANSPOSE:
      {
        result_size = other.n() * this->n();
        if (other.m() == this->m())
          break;
      }
      libmesh_fallthrough();
    case RIGHT_MULTIPLY_TRANSPOSE:
      {
        result_size = other.m() * this->m();
        if (other.n() == this->n())
          break;
      }
      libmesh_fallthrough();
    default:
      libmesh_error_msg("Unknown flag selected or matrices are incompatible for multiplication.");
    }

  // For this to work, the passed arg. must actually be a DenseMatrix<T>
  const DenseMatrix<T> * const_that = cast_ptr<const DenseMatrix<T> *>(&other);

  // Also, although 'that' is logically const in this BLAS routine,
  // the PETSc BLAS interface does not specify that any of the inputs are
  // const.  To use it, I must cast away const-ness.
  DenseMatrix<T> * that = const_cast<DenseMatrix<T> *> (const_that);

  // Initialize A, B pointers for LEFT_MULTIPLY* cases
  DenseMatrix<T> * A = this;
  DenseMatrix<T> * B = that;

  // For RIGHT_MULTIPLY* cases, swap the meaning of A and B.
  // Here is a full table of combinations we can pass to BLASgemm, and what the answer is when finished:
  // pass A B   -> (Fortran) -> A^T B^T -> (C++) -> (A^T B^T)^T -> (identity) -> B A   "lt multiply"
  // pass B A   -> (Fortran) -> B^T A^T -> (C++) -> (B^T A^T)^T -> (identity) -> A B   "rt multiply"
  // pass A B^T -> (Fortran) -> A^T B   -> (C++) -> (A^T B)^T   -> (identity) -> B^T A "lt multiply t"
  // pass B^T A -> (Fortran) -> B A^T   -> (C++) -> (B A^T)^T   -> (identity) -> A B^T "rt multiply t"
  if (flag==RIGHT_MULTIPLY || flag==RIGHT_MULTIPLY_TRANSPOSE)
    std::swap(A,B);

  // transa, transb values to pass to blas
  char
    transa[] = "n",
    transb[] = "n";

  // Integer values to pass to BLAS:
  //
  // M
  // In Fortran, the number of rows of op(A),
  // In the BLAS documentation, typically known as 'M'.
  //
  // In C/C++, we set:
  // M = n_cols(A) if (transa='n')
  //     n_rows(A) if (transa='t')
  PetscBLASInt M = static_cast<PetscBLASInt>( A->n() );

  // N
  // In Fortran, the number of cols of op(B), and also the number of cols of C.
  // In the BLAS documentation, typically known as 'N'.
  //
  // In C/C++, we set:
  // N = n_rows(B) if (transb='n')
  //     n_cols(B) if (transb='t')
  PetscBLASInt N = static_cast<PetscBLASInt>( B->m() );

  // K
  // In Fortran, the number of cols of op(A), and also
  // the number of rows of op(B). In the BLAS documentation,
  // typically known as 'K'.
  //
  // In C/C++, we set:
  // K = n_rows(A) if (transa='n')
  //     n_cols(A) if (transa='t')
  PetscBLASInt K = static_cast<PetscBLASInt>( A->m() );

  // LDA (leading dimension of A). In our cases,
  // LDA is always the number of columns of A.
  PetscBLASInt LDA = static_cast<PetscBLASInt>( A->n() );

  // LDB (leading dimension of B).  In our cases,
  // LDB is always the number of columns of B.
  PetscBLASInt LDB = static_cast<PetscBLASInt>( B->n() );

  if (flag == LEFT_MULTIPLY_TRANSPOSE)
    {
      transb[0] = 't';
      N = static_cast<PetscBLASInt>( B->n() );
    }

  else if (flag == RIGHT_MULTIPLY_TRANSPOSE)
    {
      transa[0] = 't';
      std::swap(M,K);
    }

  // LDC (leading dimension of C).  LDC is the
  // number of columns in the solution matrix.
  PetscBLASInt LDC = M;

  // Scalar values to pass to BLAS
  //
  // scalar multiplying the whole product AB
  T alpha = 1.;

  // scalar multiplying C, which is the original matrix.
  T beta  = 0.;

  // Storage for the result
  std::vector<T> result (result_size);

  // Finally ready to call the BLAS
  BLASgemm_(transa, transb, &M, &N, &K, pPS(&alpha),
            pPS(A->get_values().data()), &LDA,
            pPS(B->get_values().data()), &LDB, pPS(&beta),
            pPS(result.data()), &LDC);

  // Update the relevant dimension for this matrix.
  switch (flag)
    {
    case LEFT_MULTIPLY:            { this->_m = other.m(); break; }
    case RIGHT_MULTIPLY:           { this->_n = other.n(); break; }
    case LEFT_MULTIPLY_TRANSPOSE:  { this->_m = other.n(); break; }
    case RIGHT_MULTIPLY_TRANSPOSE: { this->_n = other.m(); break; }
    default:
      libmesh_error_msg("Unknown flag selected.");
    }

  // Swap my data vector with the result
  this->_val.swap(result);
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::_right_multiply_transpose ( const DenseMatrix< T2 > &  B )

private

由矩阵 A 的转置右乘，该矩阵可能包含不同的数值类型。

执行右乘转置矩阵的实际操作。

模板参数

T2	矩阵元素的新数据类型。

参数

A

另一个矩阵。

如果正在执行 B*(B^T)，将执行一种优化方法；否则，将执行通用的右乘转置操作。

参数

B

要右乘的转置矩阵。

在文件 dense_matrix_impl.h 第 261 行定义.

参考 libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n() , 以及 libMesh::DenseMatrix< T >::transpose().

{
  // 检查是否正在执行 B*(B^T)
  if (this == &B)
  {
    // 复制 *this 的副本
    DenseMatrix<T> A(*this);

    // 简单但低效的方式
    // 返回 this->right_multiply_transpose(A);

    // 更有效，但代码更多的方式
    // 如果 B 是 mxn，则结果将是大小为 m x m 的方阵。
    const unsigned int n_rows = B.m();
    const unsigned int n_cols = B.n();

    // 调整 *this 的大小并将所有条目清零。
    this->resize(n_rows, n_rows);

    // 计算下三角部分
    for (unsigned int i = 0; i < n_rows; ++i)
      for (unsigned int j = 0; j <= i; ++j)
        for (unsigned int k = 0; k < n_cols; ++k) // 内积在 n_cols 上
          (*this)(i, j) += A(i, k) * A(j, k);

    // 将下三角部分复制到上三角部分
    for (unsigned int i = 0; i < n_rows; ++i)
      for (unsigned int j = i + 1; j < n_rows; ++j)
        (*this)(i, j) = (*this)(j, i);
  }
  else
  {
    DenseMatrix<T> A(*this);

    this->resize(A.m(), B.m());

    libmesh_assert_equal_to(A.n(), B.n());
    libmesh_assert_equal_to(this->m(), A.m());
    libmesh_assert_equal_to(this->n(), B.m());

    const unsigned int m_s = A.m();
    const unsigned int p_s = A.n();
    const unsigned int n_s = this->n();

    // 以这种方式执行是因为
    // 很多情况下（至少对于约束矩阵而言）B.transpose(k,j) = 0
    for (unsigned int j = 0; j < n_s; j++)
      for (unsigned int k = 0; k < p_s; k++)
        if (B.transpose(k, j) != 0.)
          for (unsigned int i = 0; i < m_s; i++)
            (*this)(i, j) += A(i, k) * B.transpose(k, j);
  }
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_helper ( char  JOBU, char  JOBVT, std::vector< Real > &  sigma_val, std::vector< Number > &  U_val, std::vector< Number > &  VT_val  )

private

实际执行SVD的辅助函数。 [ 在dense_matrix_blas_lapack.C中实现 ]

参数

JOBU	用于控制计算左奇异向量的选项。
JOBVT	用于控制计算右奇异向量的选项。
sigma_val	存储奇异值的向量。
U_val	存储左奇异向量的向量。
VT_val	存储右奇异向量的向量。

在文件 dense_matrix_blas_lapack.C 第 384 行定义.

参考 libMesh::pPS().

{

  // M (input)
  //   The number of rows of the matrix A.  M >= 0.
  // In C/C++, pass the number of *cols* of A
  PetscBLASInt M = this->n();

  // N (input)
  //   The number of columns of the matrix A.  N >= 0.
  // In C/C++, pass the number of *rows* of A
  PetscBLASInt N = this->m();

  PetscBLASInt min_MN = (M < N) ? M : N;
  PetscBLASInt max_MN = (M > N) ? M : N;

  // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  //   On entry, the M-by-N matrix A.
  //   On exit,
  //   if JOBU  = 'O', A is overwritten with the first min(m,n)
  //                   columns of U (the left singular vectors,
  //                   stored columnwise);
  //   if JOBVT = 'O', A is overwritten with the first min(m,n)
  //                   rows of V**T (the right singular vectors,
  //                   stored rowwise);
  //   if JOBU != 'O' and JOBVT != 'O', the contents of A are destroyed.

  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,M).
  PetscBLASInt LDA = M;

  // S (output) DOUBLE PRECISION array, dimension (min(M,N))
  //   The singular values of A, sorted so that S(i) >= S(i+1).
  sigma_val.resize( min_MN );

  // LDU (input)
  //     The leading dimension of the array U.  LDU >= 1; if
  //     JOBU = 'S' or 'A', LDU >= M.
  PetscBLASInt LDU = M;

  // U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
  //   (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
  //   If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
  //   if JOBU = 'S', U contains the first min(m,n) columns of U
  //   (the left singular vectors, stored columnwise);
  //   if JOBU = 'N' or 'O', U is not referenced.
  if (JOBU == 'S')
    U_val.resize( LDU*min_MN );
  else
    U_val.resize( LDU*M );

  // LDVT (input)
  //      The leading dimension of the array VT.  LDVT >= 1; if
  //      JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
  PetscBLASInt LDVT = N;
  if (JOBVT == 'S')
    LDVT = min_MN;

  // VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
  //    If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
  //    V**T;
  //    if JOBVT = 'S', VT contains the first min(m,n) rows of
  //    V**T (the right singular vectors, stored rowwise);
  //    if JOBVT = 'N' or 'O', VT is not referenced.
  VT_val.resize( LDVT*N );

  // LWORK (input)
  //       The dimension of the array WORK.
  //       LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
  //       For good performance, LWORK should generally be larger.
  //
  //       If LWORK = -1, then a workspace query is assumed; the routine
  //       only calculates the optimal size of the WORK array, returns
  //       this value as the first entry of the WORK array, and no error
  //       message related to LWORK is issued by XERBLA.
  PetscBLASInt larger = (3*min_MN+max_MN > 5*min_MN) ? 3*min_MN+max_MN : 5*min_MN;
  PetscBLASInt LWORK  = (larger > 1) ? larger : 1;


  // WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  //      On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
  //      if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
  //      superdiagonal elements of an upper bidiagonal matrix B
  //      whose diagonal is in S (not necessarily sorted). B
  //      satisfies A = U * B * VT, so it has the same singular values
  //      as A, and singular vectors related by U and VT.
  std::vector<Number> WORK( LWORK );

  // INFO (output)
  //      = 0:  successful exit.
  //      < 0:  if INFO = -i, the i-th argument had an illegal value.
  //      > 0:  if DBDSQR did not converge, INFO specifies how many
  //            superdiagonals of an intermediate bidiagonal form B
  //            did not converge to zero. See the description of WORK
  //            above for details.
  PetscBLASInt INFO = 0;

  // Ready to call the actual factorization routine through PETSc's interface.
#ifdef LIBMESH_USE_REAL_NUMBERS
  // Note that the call to LAPACKgesvd_ may modify _val
  LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, pPS(_val.data()), &LDA,
               pPS(sigma_val.data()), pPS(U_val.data()), &LDU,
               pPS(VT_val.data()), &LDVT, pPS(WORK.data()), &LWORK,
               &INFO);
#else
  // When we have LIBMESH_USE_COMPLEX_NUMBERS then we must pass an array of Complex
  // numbers to LAPACKgesvd_, but _val may contain Reals so we copy to Number below to
  // handle both the real-valued and complex-valued cases.
  std::vector<Number> val_copy(_val.size());
  for (auto i : make_range(_val.size()))
    val_copy[i] = _val[i];

  std::vector<Real> RWORK(5 * min_MN);
  LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, val_copy.data(), &LDA, sigma_val.data(), U_val.data(),
               &LDU, VT_val.data(), &LDVT, WORK.data(), &LWORK, RWORK.data(), &INFO);
#endif

  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack SVD calculation!");
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_lapack ( DenseVector< Real > &  sigma )

private

使用Lapack例程“getsvd”计算矩阵的奇异值分解。 [ 在dense_matrix_blas_lapack.C中实现 ]

参数

sigma

存储奇异值的向量。

在文件 dense_matrix_blas_lapack.C 第 277 行定义.

参考 libMesh::DenseVector< T >::resize() , 以及 libMesh::DenseVector< T >::size().

{
  // The calling sequence for dgetrf is:
  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )

  // JOBU (input)
  //      Specifies options for computing all or part of the matrix U:
  //      = 'A':  all M columns of U are returned in array U:
  //      = 'S':  the first min(m,n) columns of U (the left singular
  //              vectors) are returned in the array U;
  //      = 'O':  the first min(m,n) columns of U (the left singular
  //              vectors) are overwritten on the array A;
  //      = 'N':  no columns of U (no left singular vectors) are
  //              computed.
  char JOBU = 'N';

  // JOBVT (input)
  //       Specifies options for computing all or part of the matrix
  //       V**T:
  //       = 'A':  all N rows of V**T are returned in the array VT;
  //       = 'S':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are returned in the array VT;
  //       = 'O':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are overwritten on the array A;
  //       = 'N':  no rows of V**T (no right singular vectors) are
  //               computed.
  char JOBVT = 'N';

  std::vector<Real> sigma_val;
  std::vector<Number> U_val;
  std::vector<Number> VT_val;

  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);

  // Copy the singular values into sigma, ignore U_val and VT_val
  sigma.resize(cast_int<unsigned int>(sigma_val.size()));
  for (auto i : make_range(sigma.size()))
    sigma(i) = sigma_val[i];
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_lapack ( DenseVector< Real > &  sigma, DenseMatrix< Number > &  U, DenseMatrix< Number > &  VT  )

private

使用Lapack例程“getsvd”计算矩阵的“缩减”奇异值分解。 [ 在dense_matrix_blas_lapack.C中实现 ]

参数

sigma	存储奇异值的向量。
U	存储左奇异向量的矩阵。
VT	存储右奇异向量的矩阵的转置。

在文件 dense_matrix_blas_lapack.C 第 318 行定义.

参考 libMesh::DenseMatrix< T >::get_values(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize() , 以及 libMesh::DenseVector< T >::size().

{
  // The calling sequence for dgetrf is:
  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )

  // JOBU (input)
  //      Specifies options for computing all or part of the matrix U:
  //      = 'A':  all M columns of U are returned in array U:
  //      = 'S':  the first min(m,n) columns of U (the left singular
  //              vectors) are returned in the array U;
  //      = 'O':  the first min(m,n) columns of U (the left singular
  //              vectors) are overwritten on the array A;
  //      = 'N':  no columns of U (no left singular vectors) are
  //              computed.
  char JOBU = 'S';

  // JOBVT (input)
  //       Specifies options for computing all or part of the matrix
  //       V**T:
  //       = 'A':  all N rows of V**T are returned in the array VT;
  //       = 'S':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are returned in the array VT;
  //       = 'O':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are overwritten on the array A;
  //       = 'N':  no rows of V**T (no right singular vectors) are
  //               computed.
  char JOBVT = 'S';

  // Note: Lapack is going to compute the singular values of A^T.  If
  // A=U * S * V^T, then A^T = V * S * U^T, which means that the
  // values returned in the "U_val" array actually correspond to the
  // entries of the V matrix, and the values returned in the VT_val
  // array actually correspond to the entries of U^T.  Therefore, we
  // pass VT in the place of U and U in the place of VT below!
  std::vector<Real> sigma_val;
  int M = this->n();
  int N = this->m();
  int min_MN = (M < N) ? M : N;

  // Size user-provided storage appropriately. Inside svd_helper:
  // U_val is sized to (M x min_MN)
  // VT_val is sized to (min_MN x N)
  // So, we set up U to have the shape of "VT_val^T", and VT to
  // have the shape of "U_val^T".
  //
  // Finally, since the results are stored in column-major order by
  // Lapack, but we actually want the transpose of what Lapack
  // returns, this means (conveniently) that we don't even have to
  // copy anything after the call to _svd_helper, it should already be
  // in the correct order!
  U.resize(N, min_MN);
  VT.resize(min_MN, M);

  _svd_helper(JOBU, JOBVT, sigma_val, VT.get_values(), U.get_values());

  // Copy the singular values into sigma.
  sigma.resize(cast_int<unsigned int>(sigma_val.size()));
  for (auto i : make_range(sigma.size()))
    sigma(i) = sigma_val[i];
}

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_solve_lapack ( const DenseVector< T > &  rhs, DenseVector< T > &  x, Real  rcond  ) const

private

由svd_solve(rhs)调用。

参数

rhs	右侧向量。
x	存储解的向量。
rcond	控制哪些奇异值被视为零的参数。

在文件 dense_matrix_blas_lapack.C 第 529 行定义.

参考 libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::pPS(), libMesh::DenseVector< T >::resize() , 以及 libMesh::DenseVector< T >::size().

{
  // Since BLAS is expecting column-major storage, we first need to
  // make a transposed copy of *this, then pass it to the gelss
  // routine instead of the original.  This extra copy is kind of a
  // bummer, it might be better if we could use the full SVD to
  // compute the least-squares solution instead...  Note that it isn't
  // completely terrible either, since A_trans gets overwritten by
  // Lapack, and we usually would end up making a copy of A outside
  // the function call anyway.
  DenseMatrix<T> A_trans;
  this->get_transpose(A_trans);

  // M
  // The number of rows of the input matrix. M >= 0.
  // This is actually the number of *columns* of A_trans.
  PetscBLASInt M = A_trans.n();

  // N
  // The number of columns of the matrix A. N >= 0.
  // This is actually the number of *rows* of A_trans.
  PetscBLASInt N = A_trans.m();

  // We'll use the min and max of (M,N) several times below.
  PetscBLASInt max_MN = std::max(M,N);
  PetscBLASInt min_MN = std::min(M,N);

  // NRHS
  // The number of right hand sides, i.e., the number of columns
  // of the matrices B and X. NRHS >= 0.
  // This could later be generalized to solve for multiple right-hand
  // sides...
  PetscBLASInt NRHS = 1;

  // A is double precision array, dimension (LDA,N)
  // On entry, the M-by-N matrix A.
  // On exit, the first min(m,n) rows of A are overwritten with
  // its right singular vectors, stored rowwise.
  //
  // The data vector that will be passed to Lapack.
  std::vector<T> & A_trans_vals = A_trans.get_values();

  // LDA
  // The leading dimension of the array A.  LDA >= max(1,M).
  PetscBLASInt LDA = M;

  // B is double precision array, dimension (LDB,NRHS)
  // On entry, the M-by-NRHS right hand side matrix B.
  // On exit, B is overwritten by the N-by-NRHS solution
  // matrix X.  If m >= n and RANK = n, the residual
  // sum-of-squares for the solution in the i-th column is given
  // by the sum of squares of elements n+1:m in that column.
  //
  // Since we don't want the user's rhs vector to be overwritten by
  // the solution, we copy the rhs values into the solution vector "x"
  // now.  x needs to be long enough to hold both the (Nx1) solution
  // vector or the (Mx1) rhs, so size it to the max of those.
  x.resize(max_MN);
  for (auto i : make_range(rhs.size()))
    x(i) = rhs(i);

  // Make the syntax below simpler by grabbing a reference to this array.
  std::vector<T> & B = x.get_values();

  // LDB
  // The leading dimension of the array B. LDB >= max(1,max(M,N)).
  PetscBLASInt LDB = x.size();

  // S is double precision array, dimension (min(M,N))
  // The singular values of A in decreasing order.
  // The condition number of A in the 2-norm = S(1)/S(min(m,n)).
  std::vector<T> S(min_MN);

  // RCOND
  // Used to determine the effective rank of A.  Singular values
  // S(i) <= RCOND*S(1) are treated as zero.  If RCOND < 0, machine
  // precision is used instead.
  PetscScalar RCOND = rcond;

  // RANK
  // The effective rank of A, i.e., the number of singular values
  // which are greater than RCOND*S(1).
  PetscBLASInt RANK = 0;

  // LWORK
  // The dimension of the array WORK. LWORK >= 1, and also:
  // LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
  // For good performance, LWORK should generally be larger.
  //
  // If LWORK = -1, then a workspace query is assumed; the routine
  // only calculates the optimal size of the WORK array, returns
  // this value as the first entry of the WORK array, and no error
  // message related to LWORK is issued by XERBLA.
  //
  // The factor of 3/2 is arbitrary and is used to satisfy the "should
  // generally be larger" clause.
  PetscBLASInt LWORK = (3*min_MN + std::max(2*min_MN, std::max(max_MN, NRHS))) * 3/2;

  // WORK is double precision array, dimension (MAX(1,LWORK))
  // On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  std::vector<T> WORK(LWORK);

  // INFO
  // = 0:  successful exit
  // < 0:  if INFO = -i, the i-th argument had an illegal value.
  // > 0:  the algorithm for computing the SVD failed to converge;
  //       if INFO = i, i off-diagonal elements of an intermediate
  //       bidiagonal form did not converge to zero.
  PetscBLASInt INFO = 0;

  // LAPACKgelss_(const PetscBLASInt *, // M
  //              const PetscBLASInt *, // N
  //              const PetscBLASInt *, // NRHS
  //              PetscScalar *,        // A
  //              const PetscBLASInt *, // LDA
  //              PetscScalar *,        // B
  //              const PetscBLASInt *, // LDB
  //              PetscReal *,          // S(out) = singular values of A in increasing order
  //              const PetscReal *,    // RCOND = tolerance for singular values
  //              PetscBLASInt *,       // RANK(out) = number of "non-zero" singular values
  //              PetscScalar *,        // WORK
  //              const PetscBLASInt *, // LWORK
  //              PetscBLASInt *);      // INFO
  LAPACKgelss_(&M, &N, &NRHS, pPS(A_trans_vals.data()), &LDA,
               pPS(B.data()), &LDB, pPS(S.data()), &RCOND, &RANK,
               pPS(WORK.data()), &LWORK, &INFO);

  // Check for errors in the Lapack call
  libmesh_error_msg_if(INFO < 0, "Error, argument " << -INFO << " to LAPACKgelss_ had an illegal value.");
  libmesh_error_msg_if(INFO > 0, "The algorithm for computing the SVD failed to converge!");

  // Debugging: print singular values and information about condition number:
  // libMesh::err << "RCOND=" << RCOND << std::endl;
  // libMesh::err << "Singular values: " << std::endl;
  // for (auto i : index_range(S))
  //   libMesh::err << S[i] << std::endl;
  // libMesh::err << "The condition number of A is approximately: " << S[0]/S.back() << std::endl;

  // Lapack has already written the solution into B, but it will be
  // the wrong size for non-square problems, so we need to resize it
  // correctly.  The size of the solution vector should be the number
  // of columns of the original A matrix.  Unfortunately, resizing a
  // DenseVector currently also zeros it out (unlike a std::vector) so
  // we'll resize the underlying storage directly (the size is not
  // stored independently elsewhere).
  x.get_values().resize(this->n());
}

template<typename T >

template<typename T2 , typename T3 >

boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type libMesh::DenseMatrixBase< T >::add ( const T2  factor, const DenseMatrixBase< T3 > &  mat  )

inlineinherited

将 factor 添加到矩阵的每个元素。 这应该仅在 T += T2 * T3 是有效的 C++，且 T2 是标量的情况下工作。返回类型是 void。

参数

factor	要添加的标量因子。
mat	乘法的右操作数矩阵。

在文件 dense_matrix_base.h 第 217 行定义.

参考 libMesh::DenseMatrixBase< T >::el(), libMesh::DenseMatrixBase< T >::m() , 以及 libMesh::DenseMatrixBase< T >::n().

{
  libmesh_assert_equal_to (this->m(), mat.m());
  libmesh_assert_equal_to (this->n(), mat.n());

  for (auto j : make_range(this->n()))
    for (auto i : make_range(this->m()))
      this->el(i,j) += factor*mat.el(i,j);
}

template<typename T >

template<typename T2 , typename T3 >

boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type libMesh::DenseMatrix< T >::add ( const T2  factor, const DenseMatrix< T3 > &  mat  )

inline

将 factor 倍的矩阵 mat 添加到此矩阵。

参数

factor	缩放因子。
mat	要添加的矩阵。

在文件 dense_matrix.h 第 1135 行定义.

{
  libmesh_assert_equal_to (this->m(), mat.m());
  libmesh_assert_equal_to (this->n(), mat.n());

  for (auto i : make_range(this->m()))
    for (auto j : make_range(this->n()))
      (*this)(i,j) += factor * mat(i,j);
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::cholesky_solve	(	const DenseVector< T2 > &	b,
		DenseVector< T2 > &	x
	)

对于对称正定矩阵（SPD），进行Cholesky分解。 分解为 A = L L^T，比标准LU分解大约快两倍。 因此，如果事先知道矩阵是SPD，则可以使用此方法。如果矩阵不是SPD，则会生成错误。 Cholesky分解的一个优点是它不需要主元交换以确保稳定性。

使用 Cholesky 分解解线性系统。

注意：一旦调用cholesky_solve()，就不能通过对operator(i,j)的调用修改矩阵的条目， 并在没有先调用zero()或resize()的情况下再次调用cholesky_solve()，否则代码将跳过计算矩阵的分解， 直接进行回代步骤。这是出于效率考虑而故意这样做的，以便相同的分解可用于多个右侧。但这也可能导致一些问题，所以要小心！

注解

此方法也可用于A为实值，而x和b为复值的情况。

模板参数

T2	解向量x和右侧向量b的元素类型。

参数

b	输入向量b。
x	解向量x。

Cholesky 分解首先通过 cholesky_decompose 对矩阵进行分解， 然后使用 cholesky_back_substitute 找到解 x。

模板参数

T2	输入向量 b 的数据类型。

参数

b	输入向量 b。
x	存储结果的目标向量 x。

在文件 dense_matrix_impl.h 第 804 行定义.

{
  // 检查是否已进行先前的分解
  switch (this->_decomposition_type)
  {
    case NONE:
    {
      this->_cholesky_decompose();
      break;
    }

    case CHOLESKY:
    {
      // 已经分解，只需调用 back_substitute。
      break;
    }

    default:
      libmesh_error_msg("Error! This matrix already has a different decomposition...");
  }

  // 执行回代
  this->_cholesky_back_substitute(b, x);
}

template<typename T >

void libMesh::DenseMatrixBase< T >::condense ( const unsigned int  iv, const unsigned int  jv, const T  val, DenseVectorBase< T > &  rhs  )

protectedinherited

将矩阵的 (i,j) 条目压缩出来，强制它取值为 val。这对于在数值模拟中应用边界条件很有用。 保留矩阵的对称性。

将矩阵列约减为已知值。

参数

i	要压缩的行索引。
j	要压缩的列索引。
val	压缩出的值。
rhs	与压缩操作相关联的右侧向量。

模板参数

T

矩阵元素类型。

参数

iv	列索引。
jv	列索引。
val	已知值。
rhs	包含右侧向量的 DenseVectorBase 对象。

在文件 dense_matrix_base_impl.h 第 94 行定义.

参考 libMesh::DenseVectorBase< T >::el() , 以及 libMesh::DenseVectorBase< T >::size().

参考自 libMesh::DenseMatrix< T >::condense().

{
    libmesh_assert_equal_to(this->_m, rhs.size());
    libmesh_assert_equal_to(iv, jv);

    // 将已知值移入 RHS，并将列置零
    for (auto i : make_range(this->m()))
    {
        rhs.el(i) -= this->el(i, jv) * val;
        this->el(i, jv) = 0.;
    }

    // 将行置零
    for (auto j : make_range(this->n()))
        this->el(iv, j) = 0.;

    this->el(iv, jv) = 1.;
    rhs.el(iv) = val;
}

template<typename T>

void libMesh::DenseMatrix< T >::condense ( const unsigned int  i, const unsigned int  j, const T  val, DenseVector< T > &  rhs  )

inline

将矩阵的 (i,j) 元素凝聚出来，强制其取值为 val。

这在数值模拟中用于应用边界条件。保留矩阵的对称性。

参数

i	矩阵行索引。
j	矩阵列索引。
val	设定的值。
rhs	右侧向量。

在文件 dense_matrix.h 第 507 行定义.

参考 libMesh::DenseMatrixBase< T >::condense().

  { DenseMatrixBase<T>::condense (i, j, val, rhs); }

template<typename T>

T libMesh::DenseMatrix< T >::det

(

)

返回矩阵的行列式。

注解

通过计算LU分解，然后取对角线项的乘积实现。因此，这是一个 修改矩阵条目的非const方法。

返回

矩阵的行列式。

template<typename T >

DenseVector< T > libMesh::DenseMatrixBase< T >::diagonal ( ) const

inherited

返回矩阵的对角线。

返回稠密矩阵的对角线。

返回

矩阵对角线的 DenseVector。

模板参数

T

矩阵元素类型。

返回

返回包含对角线元素的 DenseVector 对象。

在文件 dense_matrix_base_impl.h 第 43 行定义.

参考自 libMesh::DiagonalMatrix< T >::add_matrix().

{
    DenseVector<T> ret(_m);
    for (decltype(_m) i = 0; i < _m; ++i)
        ret(i) = el(i, i);
    return ret;
}

template<typename T>

virtual T libMesh::DenseMatrix< T >::el ( const unsigned int  i, const unsigned int  j  ) const

inlinefinaloverridevirtual

获取矩阵的 (i,j) 元素。

参数

i	元素的行索引。
j	元素的列索引。

返回

矩阵的 (i,j) 元素。

实现了 libMesh::DenseMatrixBase< T >.

在文件 dense_matrix.h 第 160 行定义.

  { return (*this)(i,j); }

template<typename T>

virtual T& libMesh::DenseMatrix< T >::el ( const unsigned int  i, const unsigned int  j  )

inlinefinaloverridevirtual

获取矩阵的 (i,j) 元素的可写引用。

参数

i	元素的行索引。
j	元素的列索引。

返回

矩阵的 (i,j) 元素的可写引用。

实现了 libMesh::DenseMatrixBase< T >.

在文件 dense_matrix.h 第 171 行定义.

  { return (*this)(i,j); }

template<typename T>

void libMesh::DenseMatrix< T >::evd	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag
	)

计算一般矩阵的特征值（实部和虚部）。

警告：此函数会覆盖*this的内容！

实现需要由PETSc包装的LAPACKgeev_函数。

参数

lambda_real	存储实部特征值的向量。
lambda_imag	存储虚部特征值的向量。

template<typename T>

void libMesh::DenseMatrix< T >::evd_left	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag,
		DenseMatrix< T > &	VL
	)

计算一般矩阵 的特征值（实部和虚部）和左特征向量。

警告：此函数会覆盖*this的内容！

矩阵 的左特征向量 满足： ，其中 表示 的共轭转置。

如果第j和(j+1)个特征值形成复共轭对，则VL的第j和(j+1)列“共享”它们的实值存储，方式如下： u_j = VL(:,j) + i*VL(:,j+1) 和 u_{j+1} = VL(:,j) - i*VL(:,j+1)。

实现需要由PETSc提供的LAPACKgeev_例程。

参数

lambda_real	存储实部特征值的向量。
lambda_imag	存储虚部特征值的向量。
VL	存储左特征向量的矩阵。

template<typename T>

void libMesh::DenseMatrix< T >::evd_left_and_right	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag,
		DenseMatrix< T > &	VL,
		DenseMatrix< T > &	VR
	)

计算一般矩阵的特征值（实部和虚部），以及左右特征向量。

警告：此函数会覆盖*this的内容！

有关更多信息，请参见evd_left() 和 evd_right() 函数的文档。 实现需要由PETSc提供的LAPACKgeev_例程。

参数

lambda_real	存储实部特征值的向量。
lambda_imag	存储虚部特征值的向量。
VL	存储左特征向量的矩阵。
VR	存储右特征向量的矩阵。

template<typename T>

void libMesh::DenseMatrix< T >::evd_right	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag,
		DenseMatrix< T > &	VR
	)

计算一般矩阵 的特征值（实部和虚部）和右特征向量。

警告：此函数会覆盖*this的内容！

矩阵 的右特征向量 满足： ，其中 是其相应的特征值。

注解

如果第j和(j+1)个特征值形成复共轭对，则VR的第j和(j+1)列“共享”它们的实值存储，方式如下： v_j = VR(:,j) + i*VR(:,j+1) 和 v_{j+1} = VR(:,j) - i*VR(:,j+1)。

实现需要由PETSc提供的LAPACKgeev_例程。

参数

lambda_real	存储实部特征值的向量。
lambda_imag	存储虚部特征值的向量。
VR	存储右特征向量的矩阵。

template<typename T >

void libMesh::DenseMatrix< T >::get_principal_submatrix	(	unsigned int	sub_m,
		unsigned int	sub_n,
		DenseMatrix< T > &	dest
	)		const

将 sub_m x sub_n 的主子矩阵放入 dest 中。

获取主对角线子矩阵。

参数

sub_m	子矩阵的行数。
sub_n	子矩阵的列数。
dest	存储结果的矩阵。
sub_m	子矩阵的行数。
sub_n	子矩阵的列数。
dest	存储结果的目标矩阵。

在文件 dense_matrix_impl.h 第 529 行定义.

参考 libMesh::DenseMatrix< T >::resize().

{
  libmesh_assert((sub_m <= this->m()) && (sub_n <= this->n()));

  dest.resize(sub_m, sub_n);
  for (unsigned int i = 0; i < sub_m; i++)
    for (unsigned int j = 0; j < sub_n; j++)
      dest(i, j) = (*this)(i, j);
}

template<typename T >

void libMesh::DenseMatrix< T >::get_principal_submatrix	(	unsigned int	sub_m,
		DenseMatrix< T > &	dest
	)		const

将 sub_m x sub_m 的主子矩阵放入 dest 中。

获取主对角线子矩阵。

参数

sub_m	子矩阵的行数和列数。
dest	存储结果的矩阵。
sub_m	子矩阵的行数和列数。
dest	存储结果的目标矩阵。

在文件 dense_matrix_impl.h 第 567 行定义.

{
  get_principal_submatrix(sub_m, sub_m, dest);
}

template<typename T >

void libMesh::DenseMatrix< T >::get_transpose

(

DenseMatrix< T > &

dest

)

const

将转置矩阵存储在 dest 中。

获取矩阵的转置。

参数

dest	存储结果的矩阵。
dest	存储结果的目标矩阵。

在文件 dense_matrix_impl.h 第 578 行定义.

参考 libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n() , 以及 libMesh::DenseMatrix< T >::resize().

{
  dest.resize(this->n(), this->m());

  for (auto i : make_range(dest.m()))
    for (auto j : make_range(dest.n()))
      dest(i, j) = (*this)(j, i);
}

template<typename T>

std::vector<T>& libMesh::DenseMatrix< T >::get_values ( )

inline

返回对应于存储向量的引用的底层数据存储矢量。

警告：在使用时应谨慎（即不应更改向量的大小等），但用于与期望简单数组的低级BLAS例程进行交互很有用。

返回

对底层数据存储矢量的引用。

在文件 dense_matrix.h 第 488 行定义.

参考 libMesh::DenseMatrix< T >::_val.

参考自 libMesh::DenseMatrix< T >::_evd_lapack(), libMesh::DenseMatrix< T >::_matvec_blas(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_svd_lapack(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::EpetraMatrix< T >::add_matrix() , 以及 libMesh::PetscMatrix< T >::add_matrix().

{ return _val; }

template<typename T>

const std::vector<T>& libMesh::DenseMatrix< T >::get_values ( ) const

inline

返回底层数据存储矢量的常量引用。

返回

对底层数据存储矢量的常量引用。

在文件 dense_matrix.h 第 495 行定义.

参考 libMesh::DenseMatrix< T >::_val.

{ return _val; }

template<typename T >

auto libMesh::DenseMatrix< T >::l1_norm ( ) const -> decltype(std::abs(T(0)))

inline

返回矩阵的 l1-范数，即最大列和

这是与向量的 l1-范数兼容的自然矩阵范数，即 。

返回

矩阵的 l1-范数。

在文件 dense_matrix.h 第 1242 行定义.

参考 std::abs().

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  auto columnsum = std::abs(T(0));
  for (unsigned int i=0; i!=this->_m; i++)
    {
      columnsum += std::abs((*this)(i,0));
    }
  auto my_max = columnsum;
  for (unsigned int j=1; j!=this->_n; j++)
    {
      columnsum = 0.;
      for (unsigned int i=0; i!=this->_m; i++)
        {
          columnsum += std::abs((*this)(i,j));
        }
      my_max = (my_max > columnsum? my_max : columnsum);
    }
  return my_max;
}

template<typename T >

void libMesh::DenseMatrix< T >::left_multiply ( const DenseMatrixBase< T > &  M2 )

finaloverridevirtual

左乘以矩阵 M2。

左乘另一个矩阵（M2 * (*this)）。

参数

M2	另一个矩阵。
M2	与之左乘的矩阵。

实现了 libMesh::DenseMatrixBase< T >.

在文件 dense_matrix_impl.h 第 51 行定义.

参考 libMesh::DenseMatrixBase< T >::m() , 以及 libMesh::DenseMatrixBase< T >::n().

{
  if (this->use_blas_lapack)
    this->_multiply_blas(M2, LEFT_MULTIPLY);
  else
  {
    // M3 是 *this 调整大小之前的副本
    DenseMatrix<T> M3(*this);

    // 调整 *this 的大小以适应结果
    this->resize(M2.m(), M3.n());

    // 在基类中调用乘法函数
    this->multiply(*this, M2, M3);
  }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::left_multiply

(

const DenseMatrixBase< T2 > &

M2

)

左乘以不同类型的矩阵 M2。

左乘另一个矩阵（M2 * (*this)）。

模板参数

T2	矩阵 M2 的元素类型。

参数

M2	另一个矩阵。
M2	与之左乘的矩阵。

在文件 dense_matrix_impl.h 第 76 行定义.

参考 libMesh::DenseMatrixBase< T >::m() , 以及 libMesh::DenseMatrixBase< T >::n().

{
  // M3 是 *this 调整大小之前的副本
  DenseMatrix<T> M3(*this);

  // 调整 *this 的大小以适应结果
  this->resize(M2.m(), M3.n());

  // 在基类中调用乘法函数
  this->multiply(*this, M2, M3);
}

template<typename T >

void libMesh::DenseMatrix< T >::left_multiply_transpose

(

const DenseMatrix< T > &

A

)

用矩阵 A 的转置左乘。

左乘转置矩阵。

参数

A

要与左乘的矩阵。

如果启用了 BLAS/LAPACK，它将使用优化的 BLAS 例程；否则，它将执行 _left_multiply_transpose 函数。

参数

A

要左乘的转置矩阵。

在文件 dense_matrix_impl.h 第 96 行定义.

参考自 libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual() , 以及 libMesh::DofMap::heterogeneously_constrain_element_matrix_and_vector().

{
  if (this->use_blas_lapack)
    this->_multiply_blas(A, LEFT_MULTIPLY_TRANSPOSE);
  else
    this->_left_multiply_transpose(A);
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::left_multiply_transpose

(

const DenseMatrix< T2 > &

A

)

用包含不同数值类型的矩阵 A 的转置左乘。

左乘转置矩阵。

参数

T2	矩阵元素的其他数值类型。
A	要与左乘的矩阵。
A	要左乘的转置矩阵。

在文件 dense_matrix_impl.h 第 111 行定义.

{
  this->_left_multiply_transpose(A);
}

template<typename T >

auto libMesh::DenseMatrix< T >::linfty_norm ( ) const -> decltype(std::abs(T(0)))

inline

返回矩阵的 linfty-范数，即最大行和

这是与向量的 linfty-范数兼容的自然矩阵范数，即 。

返回

矩阵的 linfty-范数。

在文件 dense_matrix.h 第 1269 行定义.

参考 std::abs().

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  auto rowsum = std::abs(T(0));
  for (unsigned int j=0; j!=this->_n; j++)
    {
      rowsum += std::abs((*this)(0,j));
    }
  auto my_max = rowsum;
  for (unsigned int i=1; i!=this->_m; i++)
    {
      rowsum = 0.;
      for (unsigned int j=0; j!=this->_n; j++)
        {
          rowsum += std::abs((*this)(i,j));
        }
      my_max = (my_max > rowsum? my_max : rowsum);
    }
  return my_max;
}

template<typename T >

void libMesh::DenseMatrix< T >::lu_solve	(	const DenseVector< T > &	b,
		DenseVector< T > &	x
	)

解方程组 给定输入向量 b。

使用 LU 分解解线性系统。

默认情况下执行部分主元消去以保持算法对舍入误差的影响稳定。

注意：一旦调用 lu_solve()，就不能通过对 operator(i, j) 的调用修改矩阵的条目， 也不能在没有先调用 zero() 或 resize() 的情况下再次调用 lu_solve()， 否则代码将跳过计算矩阵的分解而直接进入回代步骤。 这是出于效率考虑而故意这样做的，以便相同的LU分解可以用于多个右侧。 但这也可能导致一些问题，所以要小心！

参数

b	输入向量 b。
x	解向量 x。
b	输入向量 b。
x	存储结果的目标向量 x。 使用 LU 分解解线性系统。

在进行 LU 求解之前，确保矩阵是方阵。通常情况下，可以计算非方阵的 LU 分解，但是：

• 过定定系统（m>n）只有在足够的方程线性相关时才有解。
• 欠定定系统（m<n）通常有无穷多个解。

不希望在这里处理这两种不明确的情况

参数

b	输入向量 b。
x	存储结果的目标向量 x。

在文件 dense_matrix_impl.h 第 607 行定义.

{
  // 确保矩阵是方阵
  libmesh_assert_equal_to(this->m(), this->n());

  switch (this->_decomposition_type)
  {
    case NONE:
    {
      if (this->use_blas_lapack)
        this->_lu_decompose_lapack();
      else
        this->_lu_decompose();
      break;
    }

    case LU_BLAS_LAPACK:
    {
      // 已经分解，只需调用 back_substitute。
      if (this->use_blas_lapack)
        break;
    }
    libmesh_fallthrough();

    case LU:
    {
      // 已经分解，只需调用 back_substitute。
      if (!(this->use_blas_lapack))
        break;
    }
    libmesh_fallthrough();

    default:
      libmesh_error_msg("Error! This matrix already has a different decomposition...");
  }

  if (this->use_blas_lapack)
    this->_lu_back_substitute_lapack(b, x);
  else
    this->_lu_back_substitute(b, x);
}

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::m ( ) const

inlineinherited

返回矩阵的行维度。

返回

矩阵的行维度。

在文件 dense_matrix_base.h 第 115 行定义.

参考 libMesh::DenseMatrixBase< T >::_m.

参考自 libMesh::DenseMatrix< T >::_left_multiply_transpose(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_right_multiply_transpose(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::DenseMatrixBase< T >::add(), libMesh::SparseMatrix< T >::add_block_matrix(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::DiagonalMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::PetscMatrix< T >::add_matrix(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::constrain_element_vector(), libMesh::DofMap::extract_local_vector(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), libMesh::DofMap::heterogeneously_constrain_element_matrix_and_vector(), libMesh::DofMap::heterogeneously_constrain_element_residual(), libMesh::DofMap::heterogeneously_constrain_element_vector(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::DenseMatrix< T >::right_multiply() , 以及 MetaPhysicL::RawType< libMesh::DenseMatrix< T > >::value().

{ return _m; }

template<typename T >

auto libMesh::DenseMatrix< T >::max ( ) const -> decltype(libmesh_real(T(0)))

inline

返回矩阵的最大元素或复数情况下的最大实部。

返回

矩阵的最大元素或复数情况下的最大实部。

在文件 dense_matrix.h 第 1221 行定义.

参考 libMesh::libmesh_real().

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  auto my_max = libmesh_real((*this)(0,0));

  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          auto current = libmesh_real((*this)(i,j));
          my_max = (my_max > current? my_max : current);
        }
    }
  return my_max;
}

template<typename T >

auto libMesh::DenseMatrix< T >::min ( ) const -> decltype(libmesh_real(T(0)))

inline

返回矩阵的最小元素或复数情况下的最小实部。

返回

矩阵的最小元素或复数情况下的最小实部。

在文件 dense_matrix.h 第 1200 行定义.

参考 libMesh::libmesh_real().

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  auto my_min = libmesh_real((*this)(0,0));

  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          auto current = libmesh_real((*this)(i,j));
          my_min = (my_min < current? my_min : current);
        }
    }
  return my_min;
}

template<typename T >

void libMesh::DenseMatrixBase< T >::multiply ( DenseMatrixBase< T > &  M1, const DenseMatrixBase< T > &  M2, const DenseMatrixBase< T > &  M3  )

staticprotectedinherited

辅助函数 - 执行计算 M1 = M2 * M3 其中: M1 = (m x n) M2 = (m x p) M3 = (p x n)

将两个矩阵相乘。

参数

M1	乘法的结果矩阵。
M2	乘法的左操作数矩阵。
M3	乘法的右操作数矩阵。

模板参数

T

矩阵元素类型。

参数

M1	存储结果的矩阵。
M2	第一个矩阵。
M3	第二个矩阵。

在文件 dense_matrix_base_impl.h 第 61 行定义.

参考 libMesh::DenseMatrixBase< T >::el(), libMesh::DenseMatrixBase< T >::m() , 以及 libMesh::DenseMatrixBase< T >::n().

{
    // 断言确保我们传递了正确维度的矩阵
    libmesh_assert_equal_to(M1.m(), M2.m());
    libmesh_assert_equal_to(M1.n(), M3.n());
    libmesh_assert_equal_to(M2.n(), M3.m());

    const unsigned int m_s = M2.m();
    const unsigned int p_s = M2.n();
    const unsigned int n_s = M1.n();

    // 通过这种方式执行是因为有很大的机会（至少对于约束矩阵来说）
    // 在右乘时 M3(k,j) = 0.
    for (unsigned int k = 0; k < p_s; k++)
        for (unsigned int j = 0; j < n_s; j++)
            if (M3.el(k, j) != 0.)
                for (unsigned int i = 0; i < m_s; i++)
                    M1.el(i, j) += M2.el(i, k) * M3.el(k, j);
}

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::n ( ) const

inlineinherited

返回矩阵的列维度。

返回

矩阵的列维度。

在文件 dense_matrix_base.h 第 122 行定义.

参考 libMesh::DenseMatrixBase< T >::_n.

参考自 libMesh::DenseMatrix< T >::_left_multiply_transpose(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_right_multiply_transpose(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::DenseMatrixBase< T >::add(), libMesh::SparseMatrix< T >::add_block_matrix(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::DiagonalMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::PetscMatrix< T >::add_matrix(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::constrain_element_residual(), libMesh::DofMap::constrain_element_vector(), libMesh::DofMap::extract_local_vector(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), libMesh::DofMap::heterogeneously_constrain_element_matrix_and_vector(), libMesh::DofMap::heterogeneously_constrain_element_residual(), libMesh::DofMap::heterogeneously_constrain_element_vector(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::DenseMatrix< T >::right_multiply() , 以及 MetaPhysicL::RawType< libMesh::DenseMatrix< T > >::value().

{ return _n; }

template<typename T >

bool libMesh::DenseMatrix< T >::operator!= ( const DenseMatrix< T > &  mat ) const

inline

检查矩阵是否与另一个矩阵 mat 完全不相等。

参数

mat

要比较的矩阵。

返回

如果 mat 与当前矩阵完全不相等，则为 true，否则为 false。

在文件 dense_matrix.h 第 1163 行定义.

{
  for (auto i : index_range(_val))
    if (_val[i] != mat._val[i])
      return true;

  return false;
}

template<typename T >

T libMesh::DenseMatrix< T >::operator() ( const unsigned int  i, const unsigned int  j  ) const

inline

获取矩阵的 (i,j) 元素。

参数

i	元素的行索引。
j	元素的列索引。

返回

矩阵的 (i,j) 元素。

在文件 dense_matrix.h 第 1070 行定义.

{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);


  //  return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}

template<typename T >

T & libMesh::DenseMatrix< T >::operator() ( const unsigned int  i, const unsigned int  j  )

inline

获取矩阵的 (i,j) 元素的可写引用。

参数

i	元素的行索引。
j	元素的列索引。

返回

矩阵的 (i,j) 元素的可写引用。

在文件 dense_matrix.h 第 1086 行定义.

{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);

  //return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}

template<typename T >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator*= ( const T  factor )

inline

将矩阵的每个元素乘以 factor。

返回

对 *this 的引用。

在文件 dense_matrix.h 第 1122 行定义.

{
  this->scale(factor);
  return *this;
}

template<typename T >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator+= ( const DenseMatrix< T > &  mat )

inline

将矩阵 mat 添加到此矩阵。

返回

对 *this 的引用。

在文件 dense_matrix.h 第 1176 行定义.

{
  for (auto i : index_range(_val))
    _val[i] += mat._val[i];

  return *this;
}

template<typename T >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator-= ( const DenseMatrix< T > &  mat )

inline

从此矩阵中减去矩阵 mat。

返回

对 *this 的引用。

在文件 dense_matrix.h 第 1188 行定义.

{
  for (auto i : index_range(_val))
    _val[i] -= mat._val[i];

  return *this;
}

template<typename T>

DenseMatrix& libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T > &  )

default

复制赋值运算符。复制赋值运算符使用默认实现，因为该类不管理任何内存。

template<typename T>

DenseMatrix& libMesh::DenseMatrix< T >::operator= ( DenseMatrix< T > &&  )

default

移动赋值运算符。移动赋值运算符使用默认实现，因为该类不管理任何内存。

template<typename T >

template<typename T2 >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T2 > &  other_matrix )

inline

从另一种矩阵类型复制矩阵。

将类型为 T2 的矩阵复制到当前矩阵。这对于将实矩阵复制到复矩阵以进行后续操作很有用。

参数

other_matrix

要复制的矩阵。

返回

对 *this 的引用。

在文件 dense_matrix.h 第 997 行定义.

{
  unsigned int mat_m = mat.m(), mat_n = mat.n();
  this->resize(mat_m, mat_n);
  for (unsigned int i=0; i<mat_m; i++)
    for (unsigned int j=0; j<mat_n; j++)
      (*this)(i,j) = mat(i,j);

  return *this;
}

template<typename T >

bool libMesh::DenseMatrix< T >::operator== ( const DenseMatrix< T > &  mat ) const

inline

检查矩阵是否与另一个矩阵 mat 完全相等。

参数

mat

要比较的矩阵。

返回

如果 mat 与当前矩阵完全相等，则为 true，否则为 false。

在文件 dense_matrix.h 第 1150 行定义.

{
  for (auto i : index_range(_val))
    if (_val[i] != mat._val[i])
      return false;

  return true;
}

template<typename T >

void libMesh::DenseMatrix< T >::outer_product	(	const DenseVector< T > &	a,
		const DenseVector< T > &	b
	)

计算两个向量的外积并将结果存储在 (*this) 中。

外积操作。

对于两个实值向量 和 的外积是

对于两个复值向量 和 的外积是

其中 表示向量的共轭转置， 表示复共轭。

参数

a	与行相关的向量。
b	与列相关的向量。
a	输入向量 a。
b	输入向量 b。

在文件 dense_matrix_impl.h 第 548 行定义.

参考 libMesh::libmesh_conj() , 以及 libMesh::DenseVector< T >::size().

{
  const unsigned int m = a.size();
  const unsigned int n = b.size();

  this->resize(m, n);
  for (unsigned int i = 0; i < m; ++i)
    for (unsigned int j = 0; j < n; ++j)
      (*this)(i, j) = a(i) * libmesh_conj(b(j));
}

template<typename T >

void libMesh::DenseMatrixBase< T >::print ( std::ostream &  os = libMesh::out ) const

inherited

漂亮地打印矩阵，默认为 libMesh::out。

以普通格式打印矩阵。

参数

os	输出流，默认为 libMesh::out。

模板参数

T

矩阵元素类型。

参数

os

输出流。

在文件 dense_matrix_base_impl.h 第 155 行定义.

{
    for (auto i : make_range(this->m()))
    {
        for (auto j : make_range(this->n()))
            os << std::setw(8)
               << this->el(i, j) << " ";

        os << std::endl;
    }

    return;
}

template<typename T >

void libMesh::DenseMatrixBase< T >::print_scientific ( std::ostream &  os, unsigned  precision = 8  ) const

inherited

以科学计数法格式打印矩阵。

参数

os	输出流。
precision	打印的精度。

模板参数

T

矩阵元素类型。

参数

os	输出流。
precision	打印的精度。

在文件 dense_matrix_base_impl.h 第 126 行定义.

{
    // 保存初始格式标志
    std::ios_base::fmtflags os_flags = os.flags();

    // 打印矩阵元素
    for (auto i : make_range(this->m()))
    {
        for (auto j : make_range(this->n()))
            os << std::setw(15)
               << std::scientific
               << std::setprecision(precision)
               << this->el(i, j) << " ";

        os << std::endl;
    }

    // 恢复原始格式标志
    os.flags(os_flags);
}

template<typename T >

void libMesh::DenseMatrix< T >::resize ( const unsigned int  new_m, const unsigned int  new_n  )

inline

调整矩阵的大小，并调用 zero() 方法。

不会释放内存，但可能会分配更多内存。注意：当矩阵被 zero() 时，任何分解（LU、Cholesky 等）也将被清除， 强制在下一次调用例如 lu_solve() 时重新计算。

参数

new_m	新的行数。
new_n	新的列数。

在文件 dense_matrix.h 第 1012 行定义.

参考 libMesh::zero.

参考自 libMesh::DenseMatrix< T >::_evd_lapack(), libMesh::DenseMatrix< T >::_lu_decompose(), libMesh::DenseMatrix< T >::_svd_lapack(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DenseMatrix< T >::DenseMatrix(), libMesh::DenseMatrix< T >::get_principal_submatrix() , 以及 libMesh::DenseMatrix< T >::get_transpose().

{
  _val.resize(new_m*new_n);

  this->_m = new_m;
  this->_n = new_n;

  // zero and set decomposition_type to NONE
  this->zero();
}

template<typename T >

void libMesh::DenseMatrix< T >::right_multiply ( const DenseMatrixBase< T > &  M3 )

finaloverridevirtual

右乘以矩阵 M2。

右乘另一个矩阵（(*this) * M3）。

参数

M2

另一个矩阵。

如果启用了 BLAS/LAPACK，它将使用优化的 BLAS 例程；否则，它将执行基类中的乘法函数。

参数

M3

与之右乘的矩阵。

实现了 libMesh::DenseMatrixBase< T >.

在文件 dense_matrix_impl.h 第 188 行定义.

参考 libMesh::DenseMatrixBase< T >::m() , 以及 libMesh::DenseMatrixBase< T >::n().

参考自 libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual() , 以及 libMesh::DofMap::heterogeneously_constrain_element_matrix_and_vector().

{
  if (this->use_blas_lapack)
    this->_multiply_blas(M3, RIGHT_MULTIPLY);
  else
  {
    // M2 是 *this 调整大小之前的副本
    DenseMatrix<T> M2(*this);

    // 调整 *this 的大小以适应结果
    this->resize(M2.m(), M3.n());

    this->multiply(*this, M2, M3);
  }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::right_multiply

(

const DenseMatrixBase< T2 > &

M3

)

右乘以不同类型的矩阵 M2。

右乘另一个矩阵（(*this) * M3）。

模板参数

T2	矩阵 M2 的元素类型。

参数

M2

另一个矩阵。

如果启用了 BLAS/LAPACK，它将使用优化的 BLAS 例程；否则，它将执行基类中的乘法函数。

参数

M3

与之右乘的矩阵。

在文件 dense_matrix_impl.h 第 213 行定义.

参考 libMesh::DenseMatrixBase< T >::m() , 以及 libMesh::DenseMatrixBase< T >::n().

{
  // M2 是 *this 调整大小之前的副本
  DenseMatrix<T> M2(*this);

  // 调整 *this 的大小以适应结果
  this->resize(M2.m(), M3.n());

  this->multiply(*this, M2, M3);
}

template<typename T >

void libMesh::DenseMatrix< T >::right_multiply_transpose

(

const DenseMatrix< T > &

B

)

用矩阵 A 的转置右乘。

右乘转置矩阵。

参数

A

要与右乘的矩阵。

如果启用了 BLAS/LAPACK，它将使用优化的 BLAS 例程；否则，它将执行 _right_multiply_transpose 函数。

参数

B

要右乘的转置矩阵。

在文件 dense_matrix_impl.h 第 232 行定义.

{
  if (this->use_blas_lapack)
    this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
  else
    this->_right_multiply_transpose(B);
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::right_multiply_transpose

(

const DenseMatrix< T2 > &

B

)

用包含不同数值类型的矩阵 A 的转置右乘。

右乘转置矩阵。

参数

T2	矩阵元素的其他数值类型。
A	要与右乘的矩阵。
B	要右乘的转置矩阵。

在文件 dense_matrix_impl.h 第 247 行定义.

{
  this->_right_multiply_transpose(B);
}

template<typename T >

void libMesh::DenseMatrix< T >::scale ( const T  factor )

inline

将矩阵的每个元素乘以 factor。

参数

factor

乘法因子。

在文件 dense_matrix.h 第 1103 行定义.

{
  for (auto & v : _val)
    v *= factor;
}

template<typename T >

void libMesh::DenseMatrix< T >::scale_column ( const unsigned int  col, const T  factor  )

inline

将矩阵的第 col 列的每个元素乘以 factor。

参数

col	列索引。
factor	乘法因子。

在文件 dense_matrix.h 第 1112 行定义.

{
  for (auto i : make_range(this->m()))
    (*this)(i, col) *= factor;
}

template<typename T >

DenseMatrix< T > libMesh::DenseMatrix< T >::sub_matrix ( unsigned int  row_id, unsigned int  row_size, unsigned int  col_id, unsigned int  col_size  ) const

inline

获取具有最小行和列索引以及子矩阵大小的子矩阵。

参数

row_id	子矩阵的起始行索引。
row_size	子矩阵的行数。
col_id	子矩阵的起始列索引。
col_size	子矩阵的列数。

返回

具有指定大小和位置的子矩阵。

在文件 dense_matrix.h 第 1039 行定义.

{
  libmesh_assert_less (row_id + row_size - 1, this->_m);
  libmesh_assert_less (col_id + col_size - 1, this->_n);

  DenseMatrix<T> sub;
  sub._m = row_size;
  sub._n = col_size;
  sub._val.resize(row_size * col_size);

  unsigned int end_col = this->_n - col_size - col_id;
  unsigned int p = row_id * this->_n;
  unsigned int q = 0;
  for (unsigned int i=0; i<row_size; i++)
  {
    // 跳到开头的列
    p += col_id;
    for (unsigned int j=0; j<col_size; j++)
      sub._val[q++] = _val[p++];
    // 跳过剩余的列
    p += end_col;
  }

  return sub;
}

template<typename T >

void libMesh::DenseMatrix< T >::svd

(

DenseVector< Real > &

sigma

)

计算矩阵的奇异值分解（SVD）。 在退出时，sigma包含所有奇异值（按降序排列）。

使用 LAPACK svd 实现的奇异值分解。

实现使用PETSc对BLAS/LAPACK的接口。如果不可用，此函数会引发错误。

参数

sigma	存储奇异值的向量。
sigma	存储奇异值的目标向量。

在文件 dense_matrix_impl.h 第 768 行定义.

{
  // 使用 LAPACK svd 实现
  _svd_lapack(sigma);
}

template<typename T >

void libMesh::DenseMatrix< T >::svd	(	DenseVector< Real > &	sigma,
		DenseMatrix< Number > &	U,
		DenseMatrix< Number > &	VT
	)

计算矩阵的“简化”奇异值分解。 在退出时，sigma包含所有奇异值（按降序排列）， U包含左奇异向量，VT包含右奇异向量的转置。 在简化的SVD中，U有min(m,n)列，VT有min(m,n)行。（在“完整”SVD中，U和VT将是方形的。）

使用 LAPACK svd 实现的奇异值分解。

实现使用PETSc对BLAS/LAPACK的接口。如果不可用，此函数会引发错误。

参数

sigma	存储奇异值的向量。
U	存储左奇异向量的矩阵。
VT	存储右奇异向量的转置的矩阵。
sigma	存储奇异值的目标向量。
U	存储左奇异向量的目标矩阵。
VT	存储右奇异向量的目标矩阵。

在文件 dense_matrix_impl.h 第 782 行定义.

{
  // 使用 LAPACK svd 实现
  _svd_lapack(sigma, U, VT);
}

template<typename T>

void libMesh::DenseMatrix< T >::svd_solve	(	const DenseVector< T > &	rhs,
		DenseVector< T > &	x,
		Real	rcond = std::numeric_limits< Real >::epsilon()
	)		const

以最小二乘意义解方程组 。 可能是非方阵和/或秩亏的。 可以通过更改“rcond”参数来控制哪些奇异值被视为零。 对于满足S(i) <= rcond*S(1)的奇异值S(i)，在解的目的上被视为零。 传递负数的rcond将强制使用“机器精度”值。

由于各种实现细节，此函数被标记为const，我们不需要修改A的内容来计算SVD（内部会进行复制）。

需要PETSc >= 3.1，因为这是提供LAPACKgelss_包装器的第一个版本。

参数

rhs	右侧向量rhs。
x	解向量x。
rcond	控制奇异值被视为零的阈值。

template<typename T >

void libMesh::DenseMatrix< T >::swap ( DenseMatrix< T > &  other_matrix )

inline

STL 风格的交换方法,交换值和分解方法。

参数

other_matrix

要交换值的矩阵。

在文件 dense_matrix.h 第 982 行定义.

{
  std::swap(this->_m, other_matrix._m);
  std::swap(this->_n, other_matrix._n);
  _val.swap(other_matrix._val);
  DecompositionType _temp = _decomposition_type;
  _decomposition_type = other_matrix._decomposition_type;
  other_matrix._decomposition_type = _temp;
}

template<typename T >

T libMesh::DenseMatrix< T >::transpose ( const unsigned int  i, const unsigned int  j  ) const

inline

返回转置矩阵的 (i, j) 元素。

参数

i	行索引。
j	列索引。

返回

转置矩阵的 (i, j) 元素。

在文件 dense_matrix.h 第 1296 行定义.

参考自 libMesh::DenseMatrix< T >::_left_multiply_transpose() , 以及 libMesh::DenseMatrix< T >::_right_multiply_transpose().

{
  // Implement in terms of operator()
  return (*this)(j,i);
}

template<typename T >

void libMesh::DenseMatrix< T >::vector_mult	(	DenseVector< T > &	dest,
		const DenseVector< T > &	arg
	)		const

执行矩阵-向量乘法，dest := (*this) * arg。

矩阵与向量的乘法。

参数

dest	存储结果的向量。
arg	输入向量。

确保输入大小兼容。调整并清除 dest。 注意：DenseVector::resize() 也会将向量清零。

参数

dest	存储结果的目标向量。
arg	与之相乘的输入向量。

在文件 dense_matrix_impl.h 第 326 行定义.

参考 libMesh::DenseVector< T >::resize() , 以及 libMesh::DenseVector< T >::size().

参考自 libMesh::DofMap::heterogeneously_constrain_element_matrix_and_vector() , 以及 libMesh::DofMap::heterogeneously_constrain_element_vector().

{
  // 确保输入大小兼容
  libmesh_assert_equal_to(this->n(), arg.size());

  // 调整并清除 dest
  // 注意：DenseVector::resize() 也会将向量清零
  dest.resize(this->m());

  // 如果矩阵为空，直接返回
  if (this->m() == 0 || this->n() == 0)
    return;

  if (this->use_blas_lapack)
    this->_matvec_blas(1., 0., dest, arg);
  else
  {
    const unsigned int n_rows = this->m();
    const unsigned int n_cols = this->n();

    for (unsigned int i = 0; i < n_rows; i++)
      for (unsigned int j = 0; j < n_cols; j++)
        dest(i) += (*this)(i, j) * arg(j);
  }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::vector_mult	(	DenseVector< typename CompareTypes< T, T2 >::supertype > &	dest,
		const DenseVector< T2 > &	arg
	)		const

执行矩阵-向量乘法，dest := (*this) * arg ,支持不同类型的矩阵和向量。

重载的 vector_mult 函数，支持不同数据类型的向量。

模板参数

T2	输入向量 arg 的元素类型。

参数

dest	存储结果的向量。
arg	输入向量。
dest	存储结果的目标向量。
arg	与之相乘的输入向量。

在文件 dense_matrix_impl.h 第 360 行定义.

参考 libMesh::DenseVector< T >::size().

{
  // 确保输入大小兼容
  libmesh_assert_equal_to(this->n(), arg.size());

  // 调整并清除 dest
  // 注意：DenseVector::resize() 也会将向量清零
  dest.resize(this->m());

  // 如果矩阵为空，直接返回
  if (this->m() == 0 || this->n() == 0)
    return;

  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();

  for (unsigned int i = 0; i < n_rows; i++)
    for (unsigned int j = 0; j < n_cols; j++)
      dest(i) += (*this)(i, j) * arg(j);
}

template<typename T >

void libMesh::DenseMatrix< T >::vector_mult_add	(	DenseVector< T > &	dest,
		const T	factor,
		const DenseVector< T > &	arg
	)		const

执行缩放矩阵-向量乘法，结果存储在 dest 中。 dest += factor * (*this) * arg.

矩阵与向量相加。

参数

dest	结果向量，存储缩放后的矩阵乘以向量的和。
factor	缩放因子。
arg	乘法的右侧向量。

如果矩阵为空，直接返回。根据是否启用 BLAS/LAPACK 执行不同的操作。

参数

dest	存储结果的目标向量。
factor	与乘积相加的因子。
arg	与之相乘的输入向量。

在文件 dense_matrix_impl.h 第 473 行定义.

参考 libMesh::DenseVector< T >::add(), libMesh::DenseVector< T >::resize() , 以及 libMesh::DenseVector< T >::size().

参考自 libMesh::DofMap::build_constraint_matrix_and_vector().

{
  // 如果矩阵为空，直接返回
  if (this->m() == 0)
  {
    dest.resize(0);
    return;
  }

  if (this->use_blas_lapack)
    this->_matvec_blas(factor, 1., dest, arg);
  else
  {
    DenseVector<T> temp(arg.size());
    this->vector_mult(temp, arg);
    dest.add(factor, temp);
  }
}

template<typename T >

template<typename T2 , typename T3 >

void libMesh::DenseMatrix< T >::vector_mult_add	(	DenseVector< typename CompareTypes< T, typename CompareTypes< T2, T3 >::supertype >::supertype > &	dest,
		const T2	factor,
		const DenseVector< T3 > &	arg
	)		const

执行混合类型的缩放矩阵-向量乘法，结果存储在 dest 中。 dest += factor * (*this) * arg.

重载的 vector_mult_add 函数，支持不同数据类型的向量。

参数

dest	结果向量，存储缩放后的矩阵乘以向量的和。
factor	缩放因子。
arg	乘法的右侧向量。

如果矩阵为空，直接返回。根据是否启用 BLAS/LAPACK 执行不同的操作。

参数

dest	存储结果的目标向量。
factor	与乘积相加的因子。
arg	与之相乘的输入向量。

在文件 dense_matrix_impl.h 第 505 行定义.

参考 libMesh::DenseVector< T >::size().

{
  // 如果矩阵为空，直接返回
  if (this->m() == 0)
  {
    dest.resize(0);
    return;
  }

  DenseVector<typename CompareTypes<T, T3>::supertype> temp(arg.size());
  this->vector_mult(temp, arg);
  dest.add(factor, temp);
}

template<typename T >

void libMesh::DenseMatrix< T >::vector_mult_transpose	(	DenseVector< T > &	dest,
		const DenseVector< T > &	arg
	)		const

执行矩阵-向量乘法，dest := (*this)^T * arg。

矩阵转置与向量的乘法。

参数

dest	存储结果的向量。
arg	输入向量。

确保输入大小兼容。调整并清除 dest。 注意：DenseVector::resize() 也会将向量清零。

参数

dest	存储结果的目标向量。
arg	与之相乘的输入向量。

在文件 dense_matrix_impl.h 第 392 行定义.

参考 libMesh::DenseVector< T >::resize() , 以及 libMesh::DenseVector< T >::size().

参考自 libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::constrain_element_residual(), libMesh::DofMap::constrain_element_vector(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), libMesh::DofMap::heterogeneously_constrain_element_matrix_and_vector(), libMesh::DofMap::heterogeneously_constrain_element_residual() , 以及 libMesh::DofMap::heterogeneously_constrain_element_vector().

{
  // 确保输入大小兼容
  libmesh_assert_equal_to(this->m(), arg.size());

  // 调整并清除 dest
  // 注意：DenseVector::resize() 也会将向量清零
  dest.resize(this->n());

  // 如果矩阵为空，直接返回
  if (this->m() == 0)
    return;

  if (this->use_blas_lapack)
  {
    this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
  }
  else
  {
    const unsigned int n_rows = this->m();
    const unsigned int n_cols = this->n();

    // WORKS
    // for (unsigned int j=0; j<n_cols; j++)
    //   for (unsigned int i=0; i<n_rows; i++)
    //     dest(j) += (*this)(i,j)*arg(i);

    // ALSO WORKS, (i,j) just swapped
    for (unsigned int i = 0; i < n_cols; i++)
      for (unsigned int j = 0; j < n_rows; j++)
        dest(i) += (*this)(j, i) * arg(j);
  }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::vector_mult_transpose	(	DenseVector< typename CompareTypes< T, T2 >::supertype > &	dest,
		const DenseVector< T2 > &	arg
	)		const

执行矩阵-向量乘法，dest := (*this)^T * arg，支持不同类型的矩阵和向量。

重载的 vector_mult_transpose 函数，支持不同数据类型的向量。

模板参数

T2	输入向量 arg 的元素类型。

参数

dest	存储结果的向量。
arg	输入向量。
dest	存储结果的目标向量。
arg	与之相乘的输入向量。

在文件 dense_matrix_impl.h 第 435 行定义.

参考 libMesh::DenseVector< T >::size().

{
  // 确保输入大小兼容
  libmesh_assert_equal_to(this->m(), arg.size());

  // 调整并清除 dest
  // 注意：DenseVector::resize() 也会将向量清零
  dest.resize(this->n());

  // 如果矩阵为空，直接返回
  if (this->m() == 0)
    return;

  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();

  // WORKS
  // for (unsigned int j=0; j<n_cols; j++)
  //   for (unsigned int i=0; i<n_rows; i++)
  //     dest(j) += (*this)(i,j)*arg(i);

  // ALSO WORKS, (i,j) just swapped
  for (unsigned int i = 0; i < n_cols; i++)
    for (unsigned int j = 0; j < n_rows; j++)
      dest(i) += (*this)(j, i) * arg(j);
}

template<typename T >

void libMesh::DenseMatrix< T >::zero ( )

inlinefinaloverridevirtual

将矩阵的所有元素设置为0，并重置任何可能先前设置的分解标志。 这允许例如计算新的LU分解，同时重用相同的存储空间。

实现了 libMesh::DenseMatrixBase< T >.

在文件 dense_matrix.h 第 1028 行定义.

{
  _decomposition_type = NONE;

  std::fill (_val.begin(), _val.end(), static_cast<T>(0));
}

类成员变量说明

template<typename T>

DecompositionType libMesh::DenseMatrix< T >::_decomposition_type

private

此标志跟踪在矩阵上执行的分解类型。

在文件 dense_matrix.h 第 751 行定义.

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::_m

protectedinherited

行维度。

在文件 dense_matrix_base.h 第 204 行定义.

参考自 libMesh::DenseMatrixBase< T >::m().

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::_n

protectedinherited

列维度。

在文件 dense_matrix_base.h 第 209 行定义.

参考自 libMesh::DenseMatrixBase< T >::n().

template<typename T>

std::vector<pivot_index_t> libMesh::DenseMatrix< T >::_pivots

private

在文件 dense_matrix.h 第 853 行定义.

template<typename T>

std::vector<T> libMesh::DenseMatrix< T >::_val

private

实际的数据值，存储为1D数组。

在文件 dense_matrix.h 第 709 行定义.

参考自 libMesh::DenseMatrix< T >::get_values().

template<typename T>

bool libMesh::DenseMatrix< T >::use_blas_lapack

计算密矩阵的逆（假设可逆性）， 首先计算LU分解，然后执行多次回代步骤。 遵循在Web上提供的Numerical Recipes in C中的算法。

该程序被注释掉，因为它实际上不是一种内存或计算上高效的实现。 而且，通常不需要实际的逆矩阵，可以使用lu_solve()之类的东西代替。 在运行时选择性地启用/禁用BLAS支持的选项。 这主要用于测试目的

在文件 dense_matrix.h 第 694 行定义.

---

该类的文档由以下文件生成:

• /home/lwz/libmesh/include/base/dof_map.h
• /home/lwz/libmesh/include/numerics/dense_matrix.h
• /home/lwz/libmesh/src/numerics/dense_matrix.C
• /home/lwz/libmesh/src/numerics/dense_matrix_blas_lapack.C
• /home/lwz/libmesh/include/numerics/dense_matrix_impl.h