类中的 Petsc Mat 对象

计算科学 宠物 mpi 哎呀
2021-12-25 02:17:05

[Petsc 相对较新] 我正在编写一个面向对象的项目，我的想法是在用户使用 MPI 参数构造对象时拥有并行对象。所以有成员数据 Mat 并在构造函数中填充/组装它。虽然我正在努力正确地做到这一点，但这是我的第一个问题（见下面的代码。如果您需要一个最小的工作示例来调试，请告诉我）。我的第二个一般性问题是您将如何做这件事？在不破坏封装和具有一般良好实践的情况下，您将如何与该对象进行交互？

这是我在想的一个例子，首先是 hpp 文件

   class Toperator {
   public:
     /// Toperator: $T = \Del^2 + \frac{l(l+1)}{2x^2}$
     Toperator(std::shared_ptr<FEMDVR> a_radial_grid, const int &a_lmax_times_2);
  
     Toperator(const PetscMPIInt &a_numprocs, const PetscMPIInt &id,
               std::shared_ptr<FEMDVR> a_radial_grid, const int &a_lmax_times_2);
  
     /// Destructor
     ~Toperator();
  
     /// Needed to destroy the TXX Petsc Matrix
     void destroyTXXPetscMatrix();
  
     /// Needed to destroy the TIXX Petsc Matrix
     void destroyTIXXPetscMatrix();
  
     /// Returns the TXX Petsc Matrix, i.e. return type Mat
     Mat getTXXPetscMat();
  
     /// Returns the TIXX Petsc Matrix, i.e. return type Mat
     Mat getTIXXPetscMat();
  
     /// Getter for the T operator in the DVR representation.
     std::complex<double> getTXX(int index) const;
  
     /// Getter for the inverse of the T operator in the DVR representation.
     /// Used in the poison solution for $\frac{1}{|r_1-r_2|}.
     std::complex<double> getTIXX(int index) const;
  
   private:
     std::unique_ptr<std::complex<double>[]> m_dvr_rep, m_inverse_dvr_rep;
     Mat m_TXX, m_TIXX;
   };

需要注意的一点是，我想创建一个自定义析构函数来销毁 Petsc Mat 对象，但不知道如何做到这一点。一般来说，你将如何在一个类中与这个 Mat 对象交互？我只是退回整个垫子，但不确定这是否有效。不过，我很想听到更好的设计！

现在并行构造函数实现（我不确定成员列表中 Mat 的初始化）。

Toperator::Toperator(const PetscMPIInt &a_numprocs, const PetscMPIInt &id,
                     std::shared_ptr<FEMDVR> a_radial_grid,
                     const int &a_lmax_times_2)
    : m_dvr_rep(std::unique_ptr<std::complex<double>[]>(
          new std::complex<double>[a_radial_grid->getNbas() *
                                   a_radial_grid->getNbas() *
                                   a_lmax_times_2]())),
      m_inverse_dvr_rep(std::unique_ptr<std::complex<double>[]>(
          new std::complex<double>[a_radial_grid->getNbas() *
                                   a_radial_grid->getNbas() *
                                   a_lmax_times_2]())), m_TXX(nullptr), m_TIXX(nullptr) {
  PetscErrorCode ierr;
  PetscInt nbas = a_radial_grid->getNbas();
  ierr = MatCreate(PETSC_COMM_WORLD, &m_TXX);
  ierr = MatSetSizes(m_TXX, PETSC_DECIDE, PETSC_DECIDE, a_lmax_times_2,
                     nbas * nbas);
  ierr = MatSetFromOptions(m_TXX);
  ierr = MatSetUp(m_TXX);

  int start, end;
  MatGetOwnershipRange(m_TXX, &start, &end);

  for (int l = start; l < end; ++l) {
    for (int i = 0; i < nbas; ++i) {
      for (int j = 0; j < nbas; ++j) {
        if (i == j) {
          int index = i * nbas + j;
          PetscScalar tmp_TXX = a_radial_grid->getLaplacian(i * nbas + j) +
                                (std::complex<double>)(l * (l + 1)) /
                                    pow(a_radial_grid->getPoint(j), 2);
          /* MatSetValue(m_TXX, l, index, tmp_TXX, ADD_VALUES); */
          MatSetValues(m_TXX, 1, &l, 1, &index, &tmp_TXX, INSERT_VALUES);
        } else {
          int index = i * nbas + j;
          PetscScalar tmp_TXX = a_radial_grid->getLaplacian(i * nbas + j);
          /* MatSetValue(m_TXX, l, index, tmp_TXX, ADD_VALUES); */
          MatSetValues(m_TXX, 1, &l, 1, &index, &tmp_TXX, INSERT_VALUES);
        }
      }
    }
  }
  MatAssemblyBegin(m_TXX, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(m_TXX, MAT_FINAL_ASSEMBLY);
}

Toperator::~Toperator() {}

void Toperator::destroyTXXPetscMatrix() { MatDestroy(&m_TXX); }

void Toperator::destroyTIXXPetscMatrix() { MatDestroy(&m_TIXX); }

Mat Toperator::getTXXPetscMat() { return m_TXX; }

Mat Toperator::getTIXXPetscMat() { return m_TIXX; }

std::complex<double> Toperator::getTXX(int index) const {
  return m_dvr_rep[index];
}

std::complex<double> Toperator::getTIXX(int index) const {
  return m_inverse_dvr_rep[index];
}

MatSetValues 无法正常工作。这是输出。看起来只是迭代，而不是我的拉普拉斯值。

 local_row 19 local_column 8
Mat Object: 2 MPI processes
  type: mpiaij
row 0: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 1: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 2: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 3: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 4: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 5: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 6: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 7: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 8: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 9: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 10: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 11: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 12: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 13: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 14: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 15: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 16: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 17: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)
row 18: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)

这是工作顺序代码


  int nbas = a_radial_grid->getNbas();
  for (int l = 0; l < a_lmax_times_2; ++l) {
    for (int i = 0; i < nbas; ++i) {
      for (int j = 0; j < nbas; ++j) {
        if (i == j) {
          m_dvr_rep[l * nbas * nbas + i * nbas + j] =
              a_radial_grid->getLaplacian(i * nbas + j) +
              (std::complex<double>)(l * (l + 1)) /
                  pow(a_radial_grid->getPoint(j), 2);
        } else {
          m_dvr_rep[l * nbas * nbas + i * nbas + j] =
              a_radial_grid->getLaplacian(i * nbas + j);
        }
      }
    }
  }

在此代码片段中，m_dvr_rep 是我尝试并行创建的 Petsc Mat。如您所见，我正在跨越三个循环，因此我不确定如何将此结构映射到 Mat 但您可以在上面看到我的尝试。当我查看矩阵时，它只有迭代，所以我认为我在做一些愚蠢的事情。

2个回答
  1. 将 MatDestroy 放入您的析构函数中。
  2. 您设置元素的逻辑是完全错误的。最简单的解决方案是：

// 东西

for ( .... i .... )
   for ( ... j ... )
      for ( .... l ... )
         globalnumber = l * nbas * nbas + i * nbas + j
         if ( globalnumber>=low && globalnumber<high ) 
            MatSetValue
其它你可能感兴趣的问题
上一篇如果存在两个最长边，如何细化四面体？ 下一篇1D 中移动的非碰撞粒子系统