Objet Petsc Mat en classe

Question

[Relativement nouveau dans Petsc] J'écris un projet orienté objet et mon idée est d'avoir des objets parallèles lorsque l'utilisateur construit l'objet avec des arguments MPI. Ayez donc les données de membre Mat et remplissez / assemblez-les dans le constructeur. Bien que j'ai du mal à faire cela correctement, c'est ma première question (voir ci-dessous le code. Veuillez me faire savoir si vous avez besoin d'un exemple de travail minimal pour déboguer). Ma deuxième question générale est de savoir comment feriez-vous cela? Comment interagiriez-vous avec cet objet sans casser l'encapsulation et avoir de bonnes pratiques générales?

Voici un exemple de ce à quoi je pense, d'abord le fichier 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;
   };

Une chose à noter, je voudrais créer un destructeur personnalisé qui détruit les objets Petsc Mat mais je ne sais pas comment le faire. En général, comment interagiriez-vous avec cet objet Mat dans une classe? Je retourne juste le tapis entier mais je ne sais pas si cela fonctionne. Cependant, j'aimerais entendre de meilleurs designs!

Maintenant, l'implémentation du constructeur parallèle (je ne suis pas sûr de l'initialisation de Mat dans la liste des membres).

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

Les MatSetValues ​​ne fonctionnent pas correctement. Voici la sortie. Cela ressemble à des itérations et non à mes valeurs laplaciennes.

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

Voici le code séquentiel de travail

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

Dans cet extrait de code, m_dvr_rep est le tapis Petsc que j'essaie de créer en parallèle. Comme vous pouvez le voir, je traverse trois boucles, donc je ne sais pas comment mapper cette structure sur Mat mais vous pouvez voir ma tentative ci-dessus. Quand je regarde la matrice, elle ne contient que des itérations, donc je pense que je fais quelque chose de stupide.

Answer 1

Si vous cherchez un exemple, jetez un œil à la MatrixBaseclasse ici: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

Answer 2

1. Mettez le MatDestroy dans votre destructeur.
2. Votre logique de réglage des éléments est complètement fausse. La solution la plus simple est:

// truc

for ( .... i .... )
   for ( ... j ... )
      for ( .... l ... )
         globalnumber = l * nbas * nbas + i * nbas + j
         if ( globalnumber>=low && globalnumber<high ) 
            MatSetValue