2
$\begingroup$

[Relatively new to Petsc] I am writing an object oriented project and my idea is to have parallel objects when the user constructs the object with MPI arguments. So have member data Mat and fill/assembly it in the constructor. Though I am struggling doing this correctly so that's my first question (see below the code. Please let me know if you need a minimal working example to debug). My second general question is how would you go about doing this? How would you interface with this object without breaking encapsulation and having general good practices?

Here is an example of what I am thinking, first the hpp file

   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;
   };

One thing to note, I would like to create a custom destructor that destroys the Petsc Mat objects but not sure how to do that. In general, how would you interface with this Mat object in a class? I am just returning the whole Mat but not sure if this works. Though, I would love to hear better designs!

Now the parallel constructor implementation (I am not sure about the initialization of Mat in the member list).

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];
}

The MatSetValues does't work right. Here is the output. Looks like just iterations and not my laplacian values.

 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.)

Here is the working sequential code


  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);
        }
      }
    }
  }

In this code snippet, m_dvr_rep is the Petsc Mat I am trying to create in Parallel. As you can see, I am striding across three loops so I am not sure how to map this structure on to Mat but you can see my attempt above. When I view the Matrix it only has iterations in it so I think I am doing something silly.

Improve this question
$\endgroup$
0
2
$\begingroup$

If you're looking for an example, take a look at the MatrixBase class here: https://github.com/dealii/dealii/blob/master/include/deal.II/lac/petsc_matrix_base.h https://github.com/dealii/dealii/blob/master/source/lac/petsc_matrix_base.cc

Improve this answer
$\endgroup$
7
  • $\begingroup$ Take a look at the tests in the tests/petsc/ directory of that project. These are going to be the most to-the-point examples. $\endgroup$ – Wolfgang Bangerth Sep 18 '20 at 15:37
  • $\begingroup$ @zstreet Thanks for the kind words -- pretty much two decades of my academic life have gone into this :-) $\endgroup$ – Wolfgang Bangerth Sep 18 '20 at 20:51
  • $\begingroup$ It shows! And this will help accelerate my growth as well! $\endgroup$ – zstreet Sep 19 '20 at 13:05
  • $\begingroup$ @zstreet Well, thanks again the for the additional kind words :-) $\endgroup$ – Wolfgang Bangerth Sep 21 '20 at 0:50
  • $\begingroup$ How would you map my serial three loop construction of the T matrix to a parallel PETSc version? The serial version is correct so just check that out. The PETSc version is my flawed attempt... $\endgroup$ – zstreet Sep 21 '20 at 17:07
0
$\begingroup$
  1. Put the MatDestroy in your destructor.
  2. Your logic for setting elements is completely wrong. The easiest solution is:

// stuff

for ( .... i .... )
   for ( ... j ... )
      for ( .... l ... )
         globalnumber = l * nbas * nbas + i * nbas + j
         if ( globalnumber>=low && globalnumber<high ) 
            MatSetValue
Improve this answer
$\endgroup$
7
  • $\begingroup$ I initially tried putting the call to MatDestroy in the destructor but it was giving errors. Can you give me an example? I'll try again... I am not following your code snippet, sorry. Can you expand with the actual MatSetValues full call? Also I am trying to use MatSetValues so each row and not MatSetValue. Note, my first time ever using these! $\endgroup$ – zstreet Sep 17 '20 at 16:04
  • $\begingroup$ also, do you have any recommendations for design? For example, would you create a getter that returns Mat like I am showing? Thank you for your help! $\endgroup$ – zstreet Sep 17 '20 at 16:10
  • $\begingroup$ The MatSetValue call is not the point. Your treatment of indices and ownership is. $\endgroup$ – Victor Eijkhout Sep 19 '20 at 9:19
  • $\begingroup$ Yes I know I am not doing it correctly. Though your snippet isn't clear to me, sorry. My example is over 3 loops l, i, and j. Should I assume this is inside the first l loop? How do you stride appropriately for three loops? Is your low high variables my start end that I get from the Ownership call? $\endgroup$ – zstreet Sep 19 '20 at 13:04
  • $\begingroup$ What are your three loops? The i,j coordinates in your domain and what? $\endgroup$ – Victor Eijkhout Sep 20 '20 at 14:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.