]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use common vector interface for Chebyshev smoother 3767/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 11 Jan 2017 14:58:28 +0000 (15:58 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 11 Jan 2017 15:24:04 +0000 (16:24 +0100)
include/deal.II/lac/precondition.h
tests/lac/precondition_chebyshev_05.cc [new file with mode: 0644]
tests/lac/precondition_chebyshev_05.with_trilinos=true.with_lapack=true.output [new file with mode: 0644]

index 6af9b36195b8c31bc1a7f29604a2de6637ab4ec8..c201cf15cd0683fe7d6af47f26262bef4c8d0e19 100644 (file)
@@ -1903,7 +1903,6 @@ namespace internal
       (void)matrix;
       (void)preconditioner;
       AssertThrow(preconditioner.get() != NULL, ExcNotInitialized());
-      AssertDimension(preconditioner->m(), matrix.m());
     }
 
     template <typename MatrixType, typename VectorType>
@@ -1925,7 +1924,10 @@ namespace internal
           // Check if we can initialize from vector that then gets set to zero
           // as the matrix will own the memory
           preconditioner->reinit(diagonal_inverse);
-          diagonal_inverse.reinit(0);
+          {
+            VectorType empty_vector;
+            diagonal_inverse.reinit(empty_vector);
+          }
 
           // This part only works in serial
           if (preconditioner->m() != matrix.m())
@@ -2050,11 +2052,14 @@ PreconditionChebyshev<MatrixType,VectorType,PreconditionerType>::clear ()
   is_initialized = false;
   theta = delta = 1.0;
   matrix_ptr = 0;
-  data.matrix_diagonal_inverse.reinit(0);
+  {
+    VectorType empty_vector;
+    data.matrix_diagonal_inverse.reinit(empty_vector);
+    update1.reinit(empty_vector);
+    update2.reinit(empty_vector);
+    update3.reinit(empty_vector);
+  }
   data.preconditioner.reset();
-  update1.reinit(0);
-  update2.reinit(0);
-  update3.reinit(0);
 }
 
 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
diff --git a/tests/lac/precondition_chebyshev_05.cc b/tests/lac/precondition_chebyshev_05.cc
new file mode 100644 (file)
index 0000000..97ba24a
--- /dev/null
@@ -0,0 +1,100 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test PreconditionChebyshev on parallel Trilinos matrix and vector
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include "../testmatrix.h"
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iomanip>
+#include <cmath>
+
+
+
+int main(int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(4);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+
+  for (unsigned int size=4; size <= 16; size *= 2)
+    {
+      unsigned int dim = (size-1)*(size-1);
+
+      deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+      // Make matrix
+      FDMatrix testproblem(size, size);
+      SparsityPattern structure(dim, dim, 5);
+      testproblem.five_point_structure(structure);
+      structure.compress();
+      SparseMatrix<double>  A(structure);
+      testproblem.five_point(A);
+      TrilinosWrappers::SparseMatrix AA;
+      AA.reinit(A);
+
+      PreconditionChebyshev<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector, TrilinosWrappers::PreconditionJacobi> cheby;
+      PreconditionChebyshev<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector, TrilinosWrappers::PreconditionJacobi>::AdditionalData cheby_data;
+      cheby_data.preconditioner.reset(new TrilinosWrappers::PreconditionJacobi());
+      cheby_data.preconditioner->initialize(AA);
+      cheby_data.degree = 10;
+      cheby_data.smoothing_range = 40;
+      cheby.initialize(AA, cheby_data);
+
+      IndexSet set(dim);
+      set.add_range(0, dim);
+      TrilinosWrappers::MPI::Vector v, tmp1, tmp2;
+      v.reinit(set, MPI_COMM_WORLD);
+      tmp1.reinit(set, MPI_COMM_WORLD);
+      tmp2.reinit(set, MPI_COMM_WORLD);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          for (unsigned int j=0; j<dim; ++j)
+            v(j) = 1. * Testing::rand()/RAND_MAX;
+
+          AA.vmult (tmp1, v);
+          cheby_data.preconditioner->vmult (tmp2, tmp1);
+          tmp2 -= v;
+          const double ilu_residual = tmp2.l2_norm();
+
+          AA.vmult (tmp1, v);
+          cheby.vmult (tmp2, tmp1);
+          tmp2 -= v;
+          const double cheby_residual = tmp2.l2_norm();
+
+          deallog << "Residual step i=" << i << ":  "
+                  << " jacobi=" << ilu_residual
+                  << ", cheby=" << cheby_residual
+                  << std::endl;
+        }
+    }
+
+  return 0;
+}
diff --git a/tests/lac/precondition_chebyshev_05.with_trilinos=true.with_lapack=true.output b/tests/lac/precondition_chebyshev_05.with_trilinos=true.with_lapack=true.output
new file mode 100644 (file)
index 0000000..7102557
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL::Size 4 Unknowns 9
+DEAL::Residual step i=0:   jacobi=1.2946, cheby=0.0432
+DEAL::Residual step i=1:   jacobi=1.3005, cheby=0.0339
+DEAL::Residual step i=2:   jacobi=0.7640, cheby=0.0260
+DEAL::Size 8 Unknowns 49
+DEAL::Residual step i=0:   jacobi=3.3132, cheby=0.1465
+DEAL::Residual step i=1:   jacobi=3.5095, cheby=0.1507
+DEAL::Residual step i=2:   jacobi=2.9701, cheby=0.1604
+DEAL::Size 16 Unknowns 225
+DEAL::Residual step i=0:   jacobi=7.5744, cheby=3.4216
+DEAL::Residual step i=1:   jacobi=7.6760, cheby=3.4354
+DEAL::Residual step i=2:   jacobi=7.0389, cheby=3.1193

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.