]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow PETScWrappers::Precondition* re-init
authorTimo Heister <timo.heister@gmail.com>
Mon, 15 Aug 2016 00:49:07 +0000 (20:49 -0400)
committerTimo Heister <timo.heister@gmail.com>
Mon, 15 Aug 2016 00:49:07 +0000 (20:49 -0400)
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
source/lac/petsc_precondition.cc
tests/petsc/reinit_preconditioner_01.cc [new file with mode: 0644]
tests/petsc/reinit_preconditioner_01.mpirun=3.output [new file with mode: 0644]
tests/petsc/reinit_preconditioner_02.cc [new file with mode: 0644]
tests/petsc/reinit_preconditioner_02.output [new file with mode: 0644]

index 957f539957033e7a651dfb358146132a69ac55e4..c872a7344ee25c90d56331480927e68cb7cd328b 100644 (file)
@@ -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
index a16898df31e2cbf50aa034a5df0fc551aecfd0d2..3f9d104c14564a06724b60ce3a53cb600169805f 100644 (file)
@@ -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<char *>(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<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -191,6 +203,8 @@ namespace PETScWrappers
   PreconditionBlockJacobi::initialize (const MatrixBase     &matrix_,
                                        const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -226,6 +240,8 @@ namespace PETScWrappers
   PreconditionSOR::initialize (const MatrixBase     &matrix_,
                                const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -271,6 +287,8 @@ namespace PETScWrappers
   PreconditionSSOR::initialize (const MatrixBase     &matrix_,
                                 const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -320,6 +338,8 @@ namespace PETScWrappers
   PreconditionEisenstat::initialize (const MatrixBase     &matrix_,
                                      const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -366,6 +386,8 @@ namespace PETScWrappers
   PreconditionICC::initialize (const MatrixBase     &matrix_,
                                const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -411,6 +433,8 @@ namespace PETScWrappers
   PreconditionILU::initialize (const MatrixBase     &matrix_,
                                const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(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<Mat>(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<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -677,6 +706,8 @@ namespace PETScWrappers
   PreconditionNone::initialize (const MatrixBase     &matrix_,
                                 const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(matrix_);
     additional_data = additional_data_;
 
@@ -722,6 +753,8 @@ namespace PETScWrappers
   PreconditionLU::initialize (const MatrixBase     &matrix_,
                               const AdditionalData &additional_data_)
   {
+    clear ();
+
     matrix = static_cast<Mat>(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 (file)
index 0000000..efef00c
--- /dev/null
@@ -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 <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/base/index_set.h>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+
+template <class PRE>
+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<PETScWrappers::PreconditionJacobi> ();
+  test<PETScWrappers::PreconditionBlockJacobi> ();
+  //test<PETScWrappers::PreconditionSOR> ();
+  //test<PETScWrappers::PreconditionSSOR> ();
+  test<PETScWrappers::PreconditionEisenstat> ();
+  //test<PETScWrappers::PreconditionICC> ();
+  //test<PETScWrappers::PreconditionILU> ();
+  //test<PETScWrappers::PreconditionLU> ();
+  test<PETScWrappers::PreconditionBoomerAMG> ();
+  test<PETScWrappers::PreconditionParaSails> ();
+  test<PETScWrappers::PreconditionNone> ();
+}
diff --git a/tests/petsc/reinit_preconditioner_01.mpirun=3.output b/tests/petsc/reinit_preconditioner_01.mpirun=3.output
new file mode 100644 (file)
index 0000000..bb2ee5d
--- /dev/null
@@ -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 (file)
index 0000000..cbbe0c3
--- /dev/null
@@ -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 <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/base/index_set.h>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+
+template <class PRE>
+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<PETScWrappers::PreconditionJacobi> ();
+  test<PETScWrappers::PreconditionBlockJacobi> ();
+  test<PETScWrappers::PreconditionSOR> ();
+  test<PETScWrappers::PreconditionSSOR> ();
+  // todo: this crashes test<PETScWrappers::PreconditionEisenstat> ();
+  test<PETScWrappers::PreconditionICC> ();
+  test<PETScWrappers::PreconditionILU> ();
+  test<PETScWrappers::PreconditionLU> ();
+  test<PETScWrappers::PreconditionBoomerAMG> ();
+  test<PETScWrappers::PreconditionParaSails> ();
+  test<PETScWrappers::PreconditionNone> ();
+}
diff --git a/tests/petsc/reinit_preconditioner_02.output b/tests/petsc/reinit_preconditioner_02.output
new file mode 100644 (file)
index 0000000..93db3a6
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.