dense_matrix_impl.h

浏览该文件的文档.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2023 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


#ifndef LIBMESH_DENSE_MATRIX_IMPL_H
#define LIBMESH_DENSE_MATRIX_IMPL_H

// C++ Includes
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath> // for sqrt

// Local Includes
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/int_range.h"
#include "libmesh/libmesh.h"

#ifdef LIBMESH_HAVE_METAPHYSICL
#include "metaphysicl/dualnumber_decl.h"
#endif

namespace libMesh
{



// ------------------------------------------------------------
// Dense Matrix 成员函数

template<typename T>
void DenseMatrix<T>::left_multiply(const DenseMatrixBase<T> & M2)
{
  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 DenseMatrix<T>::left_multiply(const DenseMatrixBase<T2> & M2)
{
  // M3 是 *this 调整大小之前的副本
  DenseMatrix<T> M3(*this);

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

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

template<typename T>
void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T> & A)
{
  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 DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T2> & A)
{
  this->_left_multiply_transpose(A);
}

template<typename T>
template<typename T2>
void DenseMatrix<T>::_left_multiply_transpose(const DenseMatrix<T2> & A)
{
  // 检查是否正在执行 (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 DenseMatrix<T>::right_multiply(const DenseMatrixBase<T> & M3)
{
  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 DenseMatrix<T>::right_multiply(const DenseMatrixBase<T2> & M3)
{
  // M2 是 *this 调整大小之前的副本
  DenseMatrix<T> M2(*this);

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

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

template<typename T>
void DenseMatrix<T>::right_multiply_transpose(const DenseMatrix<T> & B)
{
  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 DenseMatrix<T>::right_multiply_transpose(const DenseMatrix<T2> & B)
{
  this->_right_multiply_transpose(B);
}

template<typename T>
template<typename T2>
void DenseMatrix<T>::_right_multiply_transpose(const DenseMatrix<T2> & B)
{
  // 检查是否正在执行 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>
void DenseMatrix<T>::vector_mult(DenseVector<T> & dest, const DenseVector<T> & arg) const
{
  // 确保输入大小兼容
  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 DenseMatrix<T>::vector_mult(DenseVector<typename CompareTypes<T, T2>::supertype> & dest,
                                 const DenseVector<T2> & arg) const
{
  // 确保输入大小兼容
  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 DenseMatrix<T>::vector_mult_transpose(DenseVector<T> & dest,
                                           const DenseVector<T> & arg) const
{
  // 确保输入大小兼容
  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 DenseMatrix<T>::vector_mult_transpose(DenseVector<typename CompareTypes<T, T2>::supertype> & dest,
                                           const DenseVector<T2> & arg) const
{
  // 确保输入大小兼容
  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 DenseMatrix<T>::vector_mult_add(DenseVector<T> & dest,
                                     const T factor,
                                     const DenseVector<T> & arg) const
{
  // 如果矩阵为空，直接返回
  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 DenseMatrix<T>::vector_mult_add(DenseVector<typename CompareTypes<T, typename CompareTypes<T2, T3>::supertype>::supertype> & dest,
                                     const T2 factor,
                                     const DenseVector<T3> & arg) const
{
  // 如果矩阵为空，直接返回
  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 DenseMatrix<T>::get_principal_submatrix(unsigned int sub_m,
                                             unsigned int sub_n,
                                             DenseMatrix<T> & dest) const
{
  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 DenseMatrix<T>::outer_product(const DenseVector<T> & a,
                                   const DenseVector<T> & b)
{
  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 DenseMatrix<T>::get_principal_submatrix(unsigned int sub_m, DenseMatrix<T> & dest) const
{
  get_principal_submatrix(sub_m, sub_m, dest);
}

template<typename T>
void DenseMatrix<T>::get_transpose(DenseMatrix<T> & dest) const
{
  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>
void DenseMatrix<T>::lu_solve(const DenseVector<T> & b, DenseVector<T> & x)
{
  // 确保矩阵是方阵
  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>
void DenseMatrix<T>::_lu_back_substitute(const DenseVector<T> & b, DenseVector<T> & x) const
{
  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>
void DenseMatrix<T>::_lu_decompose()
{
  // 如果调用了此函数，则矩阵不应该有任何先前的分解。
  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>
void DenseMatrix<T>::svd(DenseVector<Real> & sigma)
{
  // 使用 LAPACK svd 实现
  _svd_lapack(sigma);
}

template<typename T>
void DenseMatrix<T>::svd(DenseVector<Real> & sigma,
                         DenseMatrix<Number> & U,
                         DenseMatrix<Number> & VT)
{
  // 使用 LAPACK svd 实现
  _svd_lapack(sigma, U, VT);
}



template <typename T>
template <typename T2>
void DenseMatrix<T>::cholesky_solve(const DenseVector<T2> &b, DenseVector<T2> &x)
{
  // 检查是否已进行先前的分解
  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 DenseMatrix<T>::_cholesky_decompose()
{
  // 如果调用了此函数，则矩阵不应该有任何先前的分解。
  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 <typename T2>
void DenseMatrix<T>::_cholesky_back_substitute(const DenseVector<T2> &b,
                                               DenseVector<T2> &x) const
{
  // 确保矩阵是方阵
  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);
  }
}








// This routine is commented out since it is not really a memory
// efficient implementation.  Also, you don't *need* the inverse
// for anything, instead just use lu_solve to solve Ax=b.
// template<typename T>
// void DenseMatrix<T>::inverse ()
// {
//   // First LU decompose the matrix
//   // Note that the lu_decompose routine will check to see if the
//   // matrix is square so we don't worry about it.
//   if (!this->_lu_decomposed)
//     this->_lu_decompose();

//   // A unit vector which will be used as a rhs
//   // to pick off a single value each time.
//   DenseVector<T> e;
//   e.resize(this->m());

//   // An empty vector which will be used to hold the solution
//   // to the back substitutions.
//   DenseVector<T> x;
//   x.resize(this->m());

//   // An empty dense matrix to store the resulting inverse
//   // temporarily until we can overwrite A.
//   DenseMatrix<T> inv;
//   inv.resize(this->m(), this->n());

//   // Resize the passed in matrix to hold the inverse
//   inv.resize(this->m(), this->n());

//   for (unsigned int j=0; j<this->n(); ++j)
//     {
//       e.zero();
//       e(j) = 1.;
//       this->_lu_back_substitute(e, x, false);
//       for (unsigned int i=0; i<this->n(); ++i)
// inv(i,j) = x(i);
//     }

//   // Now overwrite all the entries
//   *this = inv;
// }


} // namespace libMesh

#endif // LIBMESH_DENSE_MATRIX_IMPL_H