From b3e1c70b95d63962ea749998aa3a7806a0d5dc28 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sun, 14 Aug 2016 20:49:07 -0400 Subject: [PATCH] Allow PETScWrappers::Precondition* re-init This PR allows re-using a preconditioner object by calling initialize() more than once. Update documentation while we are here. --- include/deal.II/lac/petsc_precondition.h | 28 ++--- source/lac/petsc_precondition.cc | 39 ++++++- tests/petsc/reinit_preconditioner_01.cc | 107 ++++++++++++++++++ .../reinit_preconditioner_01.mpirun=3.output | 65 +++++++++++ tests/petsc/reinit_preconditioner_02.cc | 94 +++++++++++++++ tests/petsc/reinit_preconditioner_02.output | 31 +++++ 6 files changed, 345 insertions(+), 19 deletions(-) create mode 100644 tests/petsc/reinit_preconditioner_01.cc create mode 100644 tests/petsc/reinit_preconditioner_01.mpirun=3.output create mode 100644 tests/petsc/reinit_preconditioner_02.cc create mode 100644 tests/petsc/reinit_preconditioner_02.output diff --git a/include/deal.II/lac/petsc_precondition.h b/include/deal.II/lac/petsc_precondition.h index 957f539957..c872a7344e 100644 --- a/include/deal.II/lac/petsc_precondition.h +++ b/include/deal.II/lac/petsc_precondition.h @@ -65,6 +65,12 @@ namespace PETScWrappers */ virtual ~PreconditionerBase (); + /** + * Destroys the preconditioner, leaving an object like just after having + * called the constructor. + */ + void clear (); + /** * Apply the preconditioner once to the given src vector. */ @@ -256,9 +262,7 @@ namespace PETScWrappers * A class that implements the interface to use the PETSc SOR * preconditioner. * - * See the comment in the base class - * @ref PreconditionerBase - * for when this preconditioner may or may not work. + * @note Only works in serial with a PETScWrappers::SparseMatrix. * * @ingroup PETScWrappers * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 @@ -318,9 +322,7 @@ namespace PETScWrappers * A class that implements the interface to use the PETSc SSOR * preconditioner. * - * See the comment in the base class - * @ref PreconditionerBase - * for when this preconditioner may or may not work. + * @note Only works in serial with a PETScWrappers::SparseMatrix. * * @ingroup PETScWrappers * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 @@ -378,7 +380,7 @@ namespace PETScWrappers /** * A class that implements the interface to use the PETSc Eisenstat - * preconditioner. + * preconditioner, which implements SSOR on each processor. * * See the comment in the base class * @ref PreconditionerBase @@ -442,9 +444,7 @@ namespace PETScWrappers * A class that implements the interface to use the PETSc Incomplete * Cholesky preconditioner. * - * See the comment in the base class - * @ref PreconditionerBase - * for when this preconditioner may or may not work. + * @note Only works in serial with a PETScWrappers::SparseMatrix. * * @ingroup PETScWrappers * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 @@ -504,9 +504,7 @@ namespace PETScWrappers * A class that implements the interface to use the PETSc ILU * preconditioner. * - * See the comment in the base class - * @ref PreconditionerBase - * for when this preconditioner may or may not work. + * @note Only works in serial with a PETScWrappers::SparseMatrix. * * @ingroup PETScWrappers * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 @@ -567,9 +565,7 @@ namespace PETScWrappers * The LU decomposition is only implemented for single processor machines. * It should provide a convenient interface to another direct solver. * - * See the comment in the base class - * @ref PreconditionerBase - * for when this preconditioner may or may not work. + * @note Only works in serial with a PETScWrappers::SparseMatrix. * * @ingroup PETScWrappers * @author Oliver Kayser-Herold, 2004 diff --git a/source/lac/petsc_precondition.cc b/source/lac/petsc_precondition.cc index a16898df31..3f9d104c14 100644 --- a/source/lac/petsc_precondition.cc +++ b/source/lac/petsc_precondition.cc @@ -35,9 +35,16 @@ namespace PETScWrappers pc(NULL), matrix(NULL) {} - PreconditionerBase::~PreconditionerBase () { + clear(); + } + + void + PreconditionerBase::clear () + { + matrix = NULL; + if (pc!=NULL) { #if DEAL_II_PETSC_VERSION_LT(3,2,0) @@ -45,6 +52,7 @@ namespace PETScWrappers #else int ierr = PCDestroy(&pc); #endif + pc = NULL; AssertThrow (ierr == 0, ExcPETScError(ierr)); } } @@ -128,6 +136,8 @@ namespace PETScWrappers void PreconditionJacobi::initialize() { + AssertThrow (pc != NULL, StandardExceptions::ExcInvalidState ()); + int ierr; ierr = PCSetType (pc, const_cast(PCJACOBI)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -140,6 +150,8 @@ namespace PETScWrappers PreconditionJacobi::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -191,6 +203,8 @@ namespace PETScWrappers PreconditionBlockJacobi::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -226,6 +240,8 @@ namespace PETScWrappers PreconditionSOR::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -271,6 +287,8 @@ namespace PETScWrappers PreconditionSSOR::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -320,6 +338,8 @@ namespace PETScWrappers PreconditionEisenstat::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -366,6 +386,8 @@ namespace PETScWrappers PreconditionICC::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -411,6 +433,8 @@ namespace PETScWrappers PreconditionILU::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -531,10 +555,12 @@ namespace PETScWrappers PreconditionBoomerAMG::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { +#ifdef PETSC_HAVE_HYPRE + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; -#ifdef PETSC_HAVE_HYPRE create_pc(); initialize (); @@ -542,7 +568,8 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); #else // PETSC_HAVE_HYPRE - (void)pc; + (void)matrix_; + (void)additional_data_; Assert (false, ExcMessage ("Your PETSc installation does not include a copy of " "the hypre package necessary for this preconditioner.")); @@ -582,6 +609,8 @@ namespace PETScWrappers PreconditionParaSails::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -677,6 +706,8 @@ namespace PETScWrappers PreconditionNone::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; @@ -722,6 +753,8 @@ namespace PETScWrappers PreconditionLU::initialize (const MatrixBase &matrix_, const AdditionalData &additional_data_) { + clear (); + matrix = static_cast(matrix_); additional_data = additional_data_; diff --git a/tests/petsc/reinit_preconditioner_01.cc b/tests/petsc/reinit_preconditioner_01.cc new file mode 100644 index 0000000000..efef00c4de --- /dev/null +++ b/tests/petsc/reinit_preconditioner_01.cc @@ -0,0 +1,107 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check re-initializing a preconditioner (parallel version) + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + +template +void test () +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + + deallog << numproc << std::endl; + + // each processor owns 2 indices and all + // are ghosting Element 1 (the second) + + IndexSet local_active(numproc*2); + local_active.add_range(myid*2,myid*2+2); + IndexSet local_relevant= local_active; + local_relevant.add_range(0,1); + + DynamicSparsityPattern csp (local_relevant); + + for (unsigned int i=0; i<2*numproc; ++i) + if (local_relevant.is_element(i)) + csp.add(i,i); + + csp.add(0,1); + + + PETScWrappers::MPI::SparseMatrix mat; + mat.reinit (local_active, local_active, csp, MPI_COMM_WORLD); + + Assert(mat.n()==numproc*2, ExcInternalError()); + Assert(mat.m()==numproc*2, ExcInternalError()); + + // set local values + mat.set(myid*2,myid*2, 1.0+myid*2.0); + mat.set(myid*2+1,myid*2+1, 2.0+myid*2.0); + + mat.compress(VectorOperation::insert); + + { + PETScWrappers::MPI::Vector src, dst; + src.reinit(local_active, MPI_COMM_WORLD); + dst.reinit(local_active, MPI_COMM_WORLD); + src(myid*2) = 1.0; + src.compress(VectorOperation::insert); + + PRE pre; + pre.initialize(mat); + pre.vmult(dst, src); + dst.print(deallog.get_file_stream()); + + mat.add(0,0,1.0); + mat.compress(VectorOperation::add); + + pre.initialize(mat); + pre.vmult(dst, src); + dst.print(deallog.get_file_stream()); + } + + if (myid==0) + deallog << "OK" << std::endl; +} + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + + test (); + test (); + //test (); + //test (); + test (); + //test (); + //test (); + //test (); + test (); + test (); + test (); +} diff --git a/tests/petsc/reinit_preconditioner_01.mpirun=3.output b/tests/petsc/reinit_preconditioner_01.mpirun=3.output new file mode 100644 index 0000000000..bb2ee5df7a --- /dev/null +++ b/tests/petsc/reinit_preconditioner_01.mpirun=3.output @@ -0,0 +1,65 @@ + +DEAL:0::3 +[Proc0 0-1] 1.000e+00 0.000e+00 +[Proc0 0-1] 2.500e-01 0.000e+00 +DEAL:0::OK +DEAL:0::3 +[Proc0 0-1] 1.000e+00 0.000e+00 +[Proc0 0-1] 2.500e-01 0.000e+00 +DEAL:0::OK +DEAL:0::3 +[Proc0 0-1] 1.000e+00 0.000e+00 +[Proc0 0-1] 4.000e+00 0.000e+00 +DEAL:0::OK +DEAL:0::3 +[Proc0 0-1] 1.000e+00 0.000e+00 +[Proc0 0-1] 2.500e-01 0.000e+00 +DEAL:0::OK +DEAL:0::3 +[Proc0 0-1] 1.000e+00 0.000e+00 +[Proc0 0-1] 2.500e-01 0.000e+00 +DEAL:0::OK +DEAL:0::3 +[Proc0 0-1] 1.000e+00 0.000e+00 +[Proc0 0-1] 1.000e+00 0.000e+00 +DEAL:0::OK + +DEAL:1::3 +[Proc1 2-3] 3.333e-01 0.000e+00 +[Proc1 2-3] 3.333e-01 0.000e+00 +DEAL:1::3 +[Proc1 2-3] 3.333e-01 0.000e+00 +[Proc1 2-3] 3.333e-01 0.000e+00 +DEAL:1::3 +[Proc1 2-3] 3.000e+00 0.000e+00 +[Proc1 2-3] 3.000e+00 0.000e+00 +DEAL:1::3 +[Proc1 2-3] 3.333e-01 0.000e+00 +[Proc1 2-3] 3.333e-01 0.000e+00 +DEAL:1::3 +[Proc1 2-3] 3.333e-01 0.000e+00 +[Proc1 2-3] 3.333e-01 0.000e+00 +DEAL:1::3 +[Proc1 2-3] 1.000e+00 0.000e+00 +[Proc1 2-3] 1.000e+00 0.000e+00 + + +DEAL:2::3 +[Proc2 4-5] 2.000e-01 0.000e+00 +[Proc2 4-5] 2.000e-01 0.000e+00 +DEAL:2::3 +[Proc2 4-5] 2.000e-01 0.000e+00 +[Proc2 4-5] 2.000e-01 0.000e+00 +DEAL:2::3 +[Proc2 4-5] 5.000e+00 0.000e+00 +[Proc2 4-5] 5.000e+00 0.000e+00 +DEAL:2::3 +[Proc2 4-5] 2.000e-01 0.000e+00 +[Proc2 4-5] 2.000e-01 0.000e+00 +DEAL:2::3 +[Proc2 4-5] 2.000e-01 0.000e+00 +[Proc2 4-5] 2.000e-01 0.000e+00 +DEAL:2::3 +[Proc2 4-5] 1.000e+00 0.000e+00 +[Proc2 4-5] 1.000e+00 0.000e+00 + diff --git a/tests/petsc/reinit_preconditioner_02.cc b/tests/petsc/reinit_preconditioner_02.cc new file mode 100644 index 0000000000..cbbe0c30fb --- /dev/null +++ b/tests/petsc/reinit_preconditioner_02.cc @@ -0,0 +1,94 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// check re-initializing a preconditioner (serial version) + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void test () +{ + DynamicSparsityPattern csp (5, 5); + + for (unsigned int i=0; i<5; ++i) + csp.add(i,i); + + csp.add(0,1); + csp.add(1,0); + + PETScWrappers::SparseMatrix mat; + mat.reinit (csp); + + for (unsigned int i=0; i<5; ++i) + mat.set(i,i, 1.0+i*2.0); + mat.set(0,1, 0.1); + mat.set(1,0, 0.1); + + mat.compress(VectorOperation::insert); + + { + PETScWrappers::Vector src, dst; + src.reinit(5); + dst.reinit(5); + src(0) = 1.0; + src(1) = 2.0; + src.compress(VectorOperation::insert); + + PRE pre; + pre.initialize(mat); + pre.vmult(dst, src); + dst.print(deallog.get_file_stream()); + + mat.add(0,0,1.0); + mat.compress(VectorOperation::add); + + pre.initialize(mat); + pre.vmult(dst, src); + dst.print(deallog.get_file_stream()); + } + + deallog << "OK" << std::endl; +} + + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + + test (); + test (); + test (); + test (); + // todo: this crashes test (); + test (); + test (); + test (); + test (); + test (); + test (); +} diff --git a/tests/petsc/reinit_preconditioner_02.output b/tests/petsc/reinit_preconditioner_02.output new file mode 100644 index 0000000000..93db3a6236 --- /dev/null +++ b/tests/petsc/reinit_preconditioner_02.output @@ -0,0 +1,31 @@ + +1.000e+00 6.667e-01 0.000e+00 0.000e+00 0.000e+00 +5.000e-01 6.667e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +9.365e-01 6.355e-01 0.000e+00 0.000e+00 0.000e+00 +4.674e-01 6.511e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +9.365e-01 6.355e-01 0.000e+00 0.000e+00 0.000e+00 +4.674e-01 6.511e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +9.365e-01 6.355e-01 0.000e+00 0.000e+00 0.000e+00 +4.674e-01 6.511e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +9.365e-01 6.355e-01 0.000e+00 0.000e+00 0.000e+00 +4.674e-01 6.511e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +9.365e-01 6.355e-01 0.000e+00 0.000e+00 0.000e+00 +4.674e-01 6.511e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +9.365e-01 6.355e-01 0.000e+00 0.000e+00 0.000e+00 +4.674e-01 6.511e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +9.367e-01 6.333e-01 0.000e+00 0.000e+00 0.000e+00 +4.675e-01 6.500e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +1.000e+00 6.667e-01 0.000e+00 0.000e+00 0.000e+00 +5.000e-01 6.667e-01 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK +1.000e+00 2.000e+00 0.000e+00 0.000e+00 0.000e+00 +1.000e+00 2.000e+00 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::OK -- 2.39.5