Petsc Mat-Objekt in der Klasse
[Relativ neu bei Petsc] Ich schreibe ein objektorientiertes Projekt und meine Idee ist es, parallele Objekte zu haben, wenn der Benutzer das Objekt mit MPI-Argumenten erstellt. Lassen Sie also die Mitgliedsdaten Mat und füllen / montieren Sie sie im Konstruktor. Obwohl ich Schwierigkeiten habe, dies richtig zu machen, ist dies meine erste Frage (siehe unten unter dem Code. Bitte lassen Sie mich wissen, wenn Sie ein minimales Arbeitsbeispiel zum Debuggen benötigen). Meine zweite allgemeine Frage ist, wie würden Sie das tun? Wie würden Sie mit diesem Objekt kommunizieren, ohne die Kapselung zu unterbrechen und allgemeine bewährte Verfahren anzuwenden?
Hier ist ein Beispiel für das, was ich denke, zuerst die HPP-Datei
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;
};
Eine Sache zu beachten, ich möchte einen benutzerdefinierten Destruktor erstellen, der die Petsc Mat-Objekte zerstört, aber nicht sicher, wie das geht. Wie würden Sie im Allgemeinen mit diesem Mat-Objekt in einer Klasse kommunizieren? Ich gebe nur die ganze Matte zurück, bin mir aber nicht sicher, ob das funktioniert. Allerdings würde ich gerne bessere Designs hören!
Nun die parallele Konstruktorimplementierung (ich bin mir nicht sicher über die Initialisierung von Mat in der Mitgliederliste).
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];
}
Die MatSetValues funktionieren nicht richtig. Hier ist die Ausgabe. Sieht aus wie nur Iterationen und nicht wie meine Laplace-Werte.
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.)
Hier ist der funktionierende sequentielle 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 diesem Code-Snippet ist m_dvr_rep die Petsc-Matte, die ich parallel erstellen möchte. Wie Sie sehen können, gehe ich über drei Schleifen, daher bin ich mir nicht sicher, wie ich diese Struktur Mat zuordnen soll, aber Sie können meinen Versuch oben sehen. Wenn ich die Matrix betrachte, enthält sie nur Iterationen, also denke ich, dass ich etwas Dummes mache.
Antworten
Wenn Sie nach einem Beispiel suchen, schauen Sie sich die MatrixBaseKlasse hier an: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
- Setze den MatDestroy in deinen Destruktor ein.
- Ihre Logik zum Setzen von Elementen ist völlig falsch. Die einfachste Lösung ist:
// Zeug
for ( .... i .... )
for ( ... j ... )
for ( .... l ... )
globalnumber = l * nbas * nbas + i * nbas + j
if ( globalnumber>=low && globalnumber<high )
MatSetValue