Objeto Petsc Mat en clase

Question

[Relativamente nuevo en Petsc] Estoy escribiendo un proyecto orientado a objetos y mi idea es tener objetos paralelos cuando el usuario construye el objeto con argumentos MPI. Así que tenga los datos de los miembros Mat y llénelos / ensamble en el constructor. Aunque estoy luchando para hacer esto correctamente, esa es mi primera pregunta (vea el código debajo. Por favor, avíseme si necesita un ejemplo de trabajo mínimo para depurar). Mi segunda pregunta general es ¿cómo harías esto? ¿Cómo interactuaría con este objeto sin romper la encapsulación y sin tener buenas prácticas generales?

Aquí hay un ejemplo de lo que estoy pensando, primero el archivo 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;
   };

Una cosa a tener en cuenta, me gustaría crear un destructor personalizado que destruya los objetos Petsc Mat, pero no estoy seguro de cómo hacerlo. En general, ¿cómo interactuaría con este objeto Mat en una clase? Solo estoy devolviendo el tapete completo, pero no estoy seguro de si esto funciona. Sin embargo, ¡me encantaría escuchar mejores diseños!

Ahora la implementación del constructor paralelo (no estoy seguro de la inicialización de Mat en la lista de miembros).

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 ​​no funciona correctamente. Aquí está la salida. Parece solo iteraciones y no mis valores laplacianos.

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

Aquí está el código secuencial de trabajo

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

En este fragmento de código, m_dvr_rep es el tapete Petsc que estoy tratando de crear en paralelo. Como puede ver, estoy cruzando tres bucles, así que no estoy seguro de cómo mapear esta estructura en Mat, pero puede ver mi intento arriba. Cuando veo Matrix, solo tiene iteraciones, así que creo que estoy haciendo algo tonto.

Answer 1

Si está buscando un ejemplo, eche un vistazo a la MatrixBaseclase aquí: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. Pon el MatDestroy en tu destructor.
2. Su lógica para configurar elementos es completamente incorrecta. La solución más sencilla es:

// cosas

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