Oggetto Petsc Mat in classe
[Relativamente nuovo per Petsc] Sto scrivendo un progetto orientato agli oggetti e la mia idea è di avere oggetti paralleli quando l'utente costruisce l'oggetto con argomenti MPI. Quindi avere Mat dati membro e riempirlo / assemblarlo nel costruttore. Anche se ho difficoltà a farlo correttamente, questa è la mia prima domanda (vedi sotto il codice. Fammi sapere se hai bisogno di un esempio minimo di lavoro per eseguire il debug). La mia seconda domanda generale è come faresti a farlo? Come vorresti interfacciarti con questo oggetto senza rompere l'incapsulamento e avere buone pratiche generali?
Ecco un esempio di quello che sto pensando, prima il file 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 da notare, vorrei creare un distruttore personalizzato che distrugge gli oggetti Petsc Mat ma non sono sicuro di come farlo. In generale, come ti interfacciaeresti con questo oggetto Mat in una classe? Sto solo restituendo l'intero tappetino ma non sono sicuro che funzioni. Tuttavia, mi piacerebbe sentire progetti migliori!
Ora l'implementazione del costruttore parallelo (non sono sicuro dell'inizializzazione di Mat nell'elenco dei membri).
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 non funziona correttamente. Ecco l'output. Sembra solo iterazioni e non i miei valori laplaciani.
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.)
Ecco il codice sequenziale funzionante
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 questo frammento di codice, m_dvr_rep è il Petsc Mat che sto cercando di creare in Parallel. Come puoi vedere, sto attraversando tre cicli quindi non sono sicuro di come mappare questa struttura su Mat, ma puoi vedere il mio tentativo sopra. Quando vedo Matrix, contiene solo iterazioni, quindi penso di fare qualcosa di sciocco.
Risposte
Se stai cercando un esempio, dai un'occhiata alla MatrixBaseclasse qui: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
- Metti il MatDestroy nel tuo distruttore.
- La tua logica per impostare gli elementi è completamente sbagliata. La soluzione più semplice è:
// cose
for ( .... i .... )
for ( ... j ... )
for ( .... l ... )
globalnumber = l * nbas * nbas + i * nbas + j
if ( globalnumber>=low && globalnumber<high )
MatSetValue