Objeto Petsc Mat na aula

Question

[Relativamente novo no Petsc] Estou escrevendo um projeto orientado a objetos e minha ideia é ter objetos paralelos quando o usuário construir o objeto com argumentos MPI. Portanto, tenha os dados do membro Mat e preencha / monte-os no construtor. Embora eu esteja lutando para fazer isso corretamente, essa é minha primeira pergunta (veja o código abaixo. Por favor, deixe-me saber se você precisa de um exemplo mínimo de trabalho para depurar). Minha segunda pergunta geral é como você faria isso? Como você faria a interface com este objeto sem quebrar o encapsulamento e ter boas práticas gerais?

Aqui está um exemplo do que estou pensando, primeiro o arquivo 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;
   };

Uma coisa a ser observada: gostaria de criar um destruidor personalizado que destrua os objetos Petsc Mat, mas não tenho certeza de como fazer isso. Em geral, como você faria a interface com este objeto Mat em uma classe? Estou apenas devolvendo o tapete inteiro, mas não tenho certeza se isso funciona. Porém, eu adoraria ouvir designs melhores!

Agora, a implementação do construtor paralelo (não tenho certeza sobre a inicialização do Mat na lista de membros).

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

O MatSetValues ​​não funciona direito. Aqui está a saída. Parece apenas iterações e não meus 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.)

Aqui está o código sequencial de trabalho

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

Neste trecho de código, m_dvr_rep é o Petsc Mat que estou tentando criar em paralelo. Como você pode ver, estou percorrendo três loops, então não tenho certeza de como mapear essa estrutura para Mat, mas você pode ver minha tentativa acima. Quando vejo a Matrix, ela contém apenas iterações, então acho que estou fazendo algo bobo.

Answer 1

Se você está procurando um exemplo, dê uma olhada na MatrixBaseclasse aqui: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. Coloque o MatDestroy em seu destruidor.
2. Sua lógica para definir os elementos está completamente errada. A solução mais fácil é:

// rechear

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