]> https://gitweb.dealii.org/ - dealii.git/commitdiff
split compilation units
authorTimo Heister <timo.heister@gmail.com>
Sat, 5 Sep 2015 15:36:25 +0000 (11:36 -0400)
committerTimo Heister <timo.heister@gmail.com>
Sat, 5 Sep 2015 15:36:25 +0000 (11:36 -0400)
12 files changed:
source/fe/fe_values.cc
source/fe/fe_values.inst.in
source/fe/fe_values_inst2.cc
source/lac/CMakeLists.txt
source/lac/sparse_matrix.cc
source/lac/sparse_matrix_inst2.cc [new file with mode: 0644]
source/lac/trilinos_precondition.cc
source/lac/trilinos_precondition2.cc [new file with mode: 0644]
source/lac/trilinos_precondition3.cc [new file with mode: 0644]
source/numerics/CMakeLists.txt
source/numerics/error_estimator.cc
source/numerics/error_estimator_inst2.cc [new file with mode: 0644]

index 904e10d5cff6f6d5cbc716910e5a98417311b9ce..faa37b8b6e2c4c60e7d0fe78b60baf3b69d04ec6 100644 (file)
@@ -4051,12 +4051,9 @@ void FESubfaceValues<dim,spacedim>::do_reinit (const unsigned int face_no,
 
 
 /*------------------------------- Explicit Instantiations -------------*/
-#ifdef FE_VALUES_INSTANTIATE_PART_TWO
-#define DIM_A 3
-#define DIM_B 3
-#else
-#define DIM_A 1
-#define DIM_B 2
+#define SPLIT_INSTANTIATIONS_COUNT 2
+#ifndef SPLIT_INSTANTIATIONS_INDEX
+#define SPLIT_INSTANTIATIONS_INDEX 0
 #endif
 #include "fe_values.inst"
 
index 04d8d1400be49aa8eba0714a6caf0e8f6d9f5a63..8678b601a1f6ba10a4b48acf0422aac7cd01de07 100644 (file)
@@ -21,8 +21,6 @@
 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
-#if (deal_II_space_dimension == DIM_A) || (deal_II_space_dimension == DIM_B)
-
     template class FEValuesBase<deal_II_dimension,deal_II_space_dimension>;
     template class FEValues<deal_II_dimension,deal_II_space_dimension>;
     template class FEValuesBase<deal_II_dimension,deal_II_space_dimension>::
@@ -40,7 +38,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
       template class SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>;
       template class Tensor<2, deal_II_dimension, deal_II_space_dimension>;
       \}
-#endif
 #endif
   }
 
@@ -48,15 +45,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 for (dof_handler : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; lda : BOOL)
   {
 #if deal_II_dimension <= deal_II_space_dimension
-#if (deal_II_space_dimension == DIM_A) || (deal_II_space_dimension == DIM_B)
-
     template void FEValues<deal_II_dimension,deal_II_space_dimension>::reinit(
     const TriaIterator<DoFCellAccessor<dof_handler<deal_II_dimension,deal_II_space_dimension>, lda> >&);
     template void FEFaceValues<deal_II_dimension,deal_II_space_dimension>::reinit(
     const TriaIterator<DoFCellAccessor<dof_handler<deal_II_dimension,deal_II_space_dimension>, lda> >&, unsigned int);
     template void FESubfaceValues<deal_II_dimension,deal_II_space_dimension>::reinit(
     const TriaIterator<DoFCellAccessor<dof_handler<deal_II_dimension,deal_II_space_dimension>, lda> >&, unsigned int, unsigned int);
-#endif
 #endif
   }
 
@@ -65,7 +59,6 @@ for (dof_handler : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; deal_II
 for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
-#if (deal_II_space_dimension == DIM_A) || (deal_II_space_dimension == DIM_B)
     template
       void FEValuesViews::Scalar<deal_II_dimension, deal_II_space_dimension>
       ::get_function_values<dealii::VEC>
@@ -129,7 +122,6 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi
       void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension>
       ::get_function_divergences<dealii::VEC>
       (const dealii::VEC&, std::vector<ProductType<dealii::VEC::value_type,dealii::Tensor<1,deal_II_space_dimension> >::type>&) const;
-#endif
 #endif
   }
 
@@ -138,7 +130,6 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi
 for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
-#if (deal_II_space_dimension == DIM_A) || (deal_II_space_dimension == DIM_B)
 
     template
       void FEValuesBase<deal_II_dimension,deal_II_space_dimension>::get_function_values<VEC>
@@ -214,7 +205,6 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi
       (const VEC&, const VectorSlice<const std::vector<types::global_dof_index> >&,
        std::vector<std::vector<VEC::value_type> > &, bool) const;
 
-#endif
 #endif
   }
 
@@ -223,7 +213,6 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi
 
 for (deal_II_dimension : DIMENSIONS)
   {
-#if (deal_II_dimension == DIM_A) || (deal_II_dimension == DIM_B)
     template
         void FEValuesViews::Scalar<deal_II_dimension>::get_function_values<dealii::IndexSet>
         (const dealii::IndexSet&, std::vector<ProductType<IndexSet::value_type,double>::type> &) const;
@@ -335,7 +324,6 @@ for (deal_II_dimension : DIMENSIONS)
         ::get_function_divergences<dealii::IndexSet>
     (const dealii::IndexSet&, std::vector<ProductType<IndexSet::value_type,dealii::Tensor<1,deal_II_dimension+1> >::type>&) const;
 
-#endif
 #endif
   }
 
@@ -345,7 +333,6 @@ for (deal_II_dimension : DIMENSIONS)
 
 for (deal_II_dimension : DIMENSIONS)
   {
-#if (deal_II_dimension == DIM_A) || (deal_II_dimension == DIM_B)
     template
     void FEValuesBase<deal_II_dimension>::get_function_values<IndexSet>
       (const IndexSet&, std::vector<IndexSet::value_type>&) const;
@@ -498,6 +485,4 @@ for (deal_II_dimension : DIMENSIONS)
        std::vector<std::vector<IndexSet::value_type> > &, bool) const;
 
 #endif
-#endif
-
   }
index ad1aa2b8298f40f8b2bd8a728e3732cef2a13108..9519d43e7d9508bc55ee98736135c89fb92d67a8 100644 (file)
@@ -16,6 +16,5 @@
 // This file compiles the second half of the instantiations from fe_values.cc
 // to get the memory consumption below 1.5gb with gcc.
 
-#define FE_VALUES_INSTANTIATE_PART_TWO
-
+#define SPLIT_INSTANTIATIONS_INDEX 1
 #include "fe_values.cc"
index e5bfce3c0d4cb4ebae47da1bc88cfc2679a65199..d9decfb1122199666c322b12dbbd9201d3effe30 100644 (file)
@@ -54,6 +54,7 @@ SET(_src
   sparse_direct.cc
   sparse_ilu.cc
   sparse_matrix.cc
+  sparse_matrix_inst2.cc
   sparse_matrix_ez.cc
   sparse_mic.cc
   sparse_vanka.cc
@@ -65,6 +66,8 @@ SET(_src
   trilinos_block_vector.cc
   trilinos_parallel_block_vector.cc
   trilinos_precondition.cc
+  trilinos_precondition2.cc
+  trilinos_precondition3.cc
   trilinos_solver.cc
   trilinos_sparse_matrix.cc
   trilinos_sparsity_pattern.cc
index ea4ce4000283023307d4fbc87d91ca0ea2cb00c3..4db5f3864518817084e2e105aa19d7b9df5c663e 100644 (file)
@@ -18,6 +18,8 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+#define SPLIT_INSTANTIATIONS_COUNT 2
+#define SPLIT_INSTANTIATIONS_INDEX 0
 #include "sparse_matrix.inst"
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/source/lac/sparse_matrix_inst2.cc b/source/lac/sparse_matrix_inst2.cc
new file mode 100644 (file)
index 0000000..aaa7275
--- /dev/null
@@ -0,0 +1,26 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/lac/sparse_matrix.templates.h>
+#include <deal.II/lac/block_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#define SPLIT_INSTANTIATIONS_COUNT 2
+#define SPLIT_INSTANTIATIONS_INDEX 1
+#include "sparse_matrix.inst"
+
+DEAL_II_NAMESPACE_CLOSE
index e7582f617ee095e407199ef0ceed09a6e3db9c87..9c644c439ac45e1e28e7a5c4e820ba246f5e1e59 100644 (file)
@@ -27,14 +27,6 @@ DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
 #  include <Teuchos_ParameterList.hpp>
 #  include <Teuchos_RCP.hpp>
 #  include <Epetra_MultiVector.h>
-#  include <ml_include.h>
-#  include <ml_MultiLevelPreconditioner.h>
-
-#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
-#  include <MueLu.hpp>
-#  include <MueLu_EpetraOperator.hpp>
-#  include <MueLu_MLParameterListInterpreter.hpp>
-#endif
 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 
 DEAL_II_NAMESPACE_OPEN
@@ -709,516 +701,10 @@ namespace TrilinosWrappers
 
 
 
-  /* -------------------------- PreconditionAMG -------------------------- */
-
-  PreconditionAMG::AdditionalData::
-  AdditionalData (const bool                             elliptic,
-                  const bool                             higher_order_elements,
-                  const unsigned int                     n_cycles,
-                  const bool                             w_cycle,
-                  const double                           aggregation_threshold,
-                  const std::vector<std::vector<bool> > &constant_modes,
-                  const unsigned int                     smoother_sweeps,
-                  const unsigned int                     smoother_overlap,
-                  const bool                             output_details,
-                  const char                            *smoother_type,
-                  const char                            *coarse_type)
-    :
-    elliptic (elliptic),
-    higher_order_elements (higher_order_elements),
-    n_cycles (n_cycles),
-    w_cycle (w_cycle),
-    aggregation_threshold (aggregation_threshold),
-    constant_modes (constant_modes),
-    smoother_sweeps (smoother_sweeps),
-    smoother_overlap (smoother_overlap),
-    output_details (output_details),
-    smoother_type (smoother_type),
-    coarse_type (coarse_type)
-  {}
-
-
-  PreconditionAMG::~PreconditionAMG()
-  {
-    preconditioner.reset();
-    trilinos_matrix.reset();
-  }
-
-
-
-  void
-  PreconditionAMG:: initialize (const SparseMatrix   &matrix,
-                                const AdditionalData &additional_data)
-  {
-    initialize(matrix.trilinos_matrix(), additional_data);
-  }
-
-
-
-  void
-  PreconditionAMG:: initialize (const Epetra_RowMatrix &matrix,
-                                const AdditionalData   &additional_data)
-  {
-    // Build the AMG preconditioner.
-    Teuchos::ParameterList parameter_list;
-
-    if (additional_data.elliptic == true)
-      {
-        ML_Epetra::SetDefaults("SA",parameter_list);
-
-        // uncoupled mode can give a lot of warnings or even fail when there
-        // are too many entries per row and aggreggation gets complicated, but
-        // MIS does not work if too few elements are located on one
-        // processor. work around these warnings by choosing the different
-        // strategies in different situations: for low order, always use the
-        // standard choice uncoupled. if higher order, right now we also just
-        // use Uncoupled, but we should be aware that maybe MIS might be
-        // needed
-        if (additional_data.higher_order_elements)
-          parameter_list.set("aggregation: type", "Uncoupled");
-      }
-    else
-      {
-        ML_Epetra::SetDefaults("NSSA",parameter_list);
-        parameter_list.set("aggregation: type", "Uncoupled");
-        parameter_list.set("aggregation: block scaling", true);
-      }
-
-    parameter_list.set("smoother: type", additional_data.smoother_type);
-    parameter_list.set("coarse: type", additional_data.coarse_type);
-
-    parameter_list.set("smoother: sweeps",
-                       static_cast<int>(additional_data.smoother_sweeps));
-    parameter_list.set("cycle applications",
-                       static_cast<int>(additional_data.n_cycles));
-    if (additional_data.w_cycle == true)
-      parameter_list.set("prec type", "MGW");
-    else
-      parameter_list.set("prec type", "MGV");
-
-    parameter_list.set("smoother: Chebyshev alpha",10.);
-    parameter_list.set("smoother: ifpack overlap",
-                       static_cast<int>(additional_data.smoother_overlap));
-    parameter_list.set("aggregation: threshold",
-                       additional_data.aggregation_threshold);
-    parameter_list.set("coarse: max size", 2000);
-
-    if (additional_data.output_details)
-      parameter_list.set("ML output", 10);
-    else
-      parameter_list.set("ML output", 0);
-
-    const Epetra_Map &domain_map = matrix.OperatorDomainMap();
-
-    const size_type constant_modes_dimension =
-      additional_data.constant_modes.size();
-    Epetra_MultiVector distributed_constant_modes (domain_map,
-                                                   constant_modes_dimension > 0 ?
-                                                   constant_modes_dimension : 1);
-    std::vector<double> dummy (constant_modes_dimension);
-
-    if (constant_modes_dimension > 0)
-      {
-        const size_type n_rows = n_global_rows(matrix);
-        const bool constant_modes_are_global =
-          additional_data.constant_modes[0].size() == n_rows;
-        const size_type n_relevant_rows =
-          constant_modes_are_global ? n_rows : additional_data.constant_modes[0].size();
-        const size_type my_size = domain_map.NumMyElements();
-        if (constant_modes_are_global == false)
-          Assert (n_relevant_rows == my_size,
-                  ExcDimensionMismatch(n_relevant_rows, my_size));
-        Assert (n_rows ==
-                static_cast<size_type>(global_length(distributed_constant_modes)),
-                ExcDimensionMismatch(n_rows,
-                                     global_length(distributed_constant_modes)));
-
-        (void)n_relevant_rows;
-        (void)global_length;
-
-        // Reshape null space as a contiguous vector of doubles so that
-        // Trilinos can read from it.
-        for (size_type d=0; d<constant_modes_dimension; ++d)
-          for (size_type row=0; row<my_size; ++row)
-            {
-              TrilinosWrappers::types::int_type global_row_id =
-                constant_modes_are_global ? gid(domain_map,row) : row;
-              distributed_constant_modes[d][row] =
-                additional_data.constant_modes[d][global_row_id];
-            }
-
-        parameter_list.set("null space: type", "pre-computed");
-        parameter_list.set("null space: dimension",
-                           distributed_constant_modes.NumVectors());
-        if (my_size > 0)
-          parameter_list.set("null space: vectors",
-                             distributed_constant_modes.Values());
-        // We need to set a valid pointer to data even if there is no data on
-        // the current processor. Therefore, pass a dummy in that case
-        else
-          parameter_list.set("null space: vectors",
-                             &dummy[0]);
-      }
-
-    initialize (matrix, parameter_list);
-
-    if (additional_data.output_details)
-      {
-        ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
-          dynamic_cast<ML_Epetra::MultiLevelPreconditioner *> (preconditioner.get());
-        Assert (multilevel_operator != 0,
-                ExcMessage ("Preconditioner setup failed."));
-        multilevel_operator->PrintUnused(0);
-      }
-  }
-
-
-
-  void
-  PreconditionAMG::initialize (const SparseMatrix           &matrix,
-                               const Teuchos::ParameterList &ml_parameters)
-  {
-    initialize(matrix.trilinos_matrix(), ml_parameters);
-  }
-
-
-
-  void
-  PreconditionAMG::initialize (const Epetra_RowMatrix       &matrix,
-                               const Teuchos::ParameterList &ml_parameters)
-  {
-    preconditioner.reset ();
-    preconditioner.reset (new ML_Epetra::MultiLevelPreconditioner
-                          (matrix, ml_parameters));
-  }
-
-
-
-  template <typename number>
-  void
-  PreconditionAMG::
-  initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
-              const AdditionalData                 &additional_data,
-              const double                          drop_tolerance,
-              const ::dealii::SparsityPattern      *use_this_sparsity)
-  {
-    preconditioner.reset();
-    const size_type n_rows = deal_ii_sparse_matrix.m();
-
-    // Init Epetra Matrix using an
-    // equidistributed map; avoid
-    // storing the nonzero
-    // elements.
-    vector_distributor.reset (new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n_rows),
-                                             0, communicator));
-
-    if (trilinos_matrix.get() == 0)
-      trilinos_matrix.reset (new SparseMatrix());
-
-    trilinos_matrix->reinit (*vector_distributor, *vector_distributor,
-                             deal_ii_sparse_matrix, drop_tolerance, true,
-                             use_this_sparsity);
-
-    initialize (*trilinos_matrix, additional_data);
-  }
-
-
-
-  void PreconditionAMG::reinit ()
-  {
-    ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
-      dynamic_cast<ML_Epetra::MultiLevelPreconditioner *> (preconditioner.get());
-    multilevel_operator->ReComputePreconditioner();
-  }
-
-
-
-  void PreconditionAMG::clear ()
-  {
-    PreconditionBase::clear();
-    trilinos_matrix.reset();
-  }
-
-
-
-  PreconditionAMG::size_type
-  PreconditionAMG::memory_consumption() const
-  {
-    unsigned int memory = sizeof(this);
-
-    // todo: find a way to read out ML's data
-    // sizes
-    if (trilinos_matrix.get() != 0)
-      memory += trilinos_matrix->memory_consumption();
-    return memory;
-  }
-
-
-
 
-  // explicit instantiations
-  template void PreconditionAMG::initialize (const ::dealii::SparseMatrix<double> &,
-                                             const AdditionalData &, const double,
-                                             const ::dealii::SparsityPattern *);
-  template void PreconditionAMG::initialize (const ::dealii::SparseMatrix<float> &,
-                                             const AdditionalData &, const double,
-                                             const ::dealii::SparsityPattern *);
 
 
 
-#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
-  /* -------------------------- PreconditionAMGMueLu --------------------- */
-
-  PreconditionAMGMueLu::AdditionalData::
-  AdditionalData (const bool                             elliptic,
-                  const unsigned int                     n_cycles,
-                  const bool                             w_cycle,
-                  const double                           aggregation_threshold,
-                  const std::vector<std::vector<bool> > &constant_modes,
-                  const unsigned int                     smoother_sweeps,
-                  const unsigned int                     smoother_overlap,
-                  const bool                             output_details,
-                  const char                            *smoother_type,
-                  const char                            *coarse_type)
-    :
-    elliptic (elliptic),
-    n_cycles (n_cycles),
-    w_cycle (w_cycle),
-    aggregation_threshold (aggregation_threshold),
-    constant_modes (constant_modes),
-    smoother_sweeps (smoother_sweeps),
-    smoother_overlap (smoother_overlap),
-    output_details (output_details),
-    smoother_type (smoother_type),
-    coarse_type (coarse_type)
-  {}
-
-
-  PreconditionAMGMueLu::~PreconditionAMGMueLu()
-  {
-    preconditioner.reset();
-    trilinos_matrix.reset();
-  }
-
-
-
-  void
-  PreconditionAMGMueLu::initialize (const SparseMatrix   &matrix,
-                                    const AdditionalData &additional_data)
-  {
-    initialize(matrix.trilinos_matrix(), additional_data);
-  }
-
-
-
-  void
-  PreconditionAMGMueLu::initialize (const Epetra_CrsMatrix &matrix,
-                                    const AdditionalData   &additional_data)
-  {
-    // Build the AMG preconditioner.
-    Teuchos::ParameterList parameter_list;
-
-    if (additional_data.elliptic == true)
-      ML_Epetra::SetDefaults("SA",parameter_list);
-    else
-      {
-        ML_Epetra::SetDefaults("NSSA",parameter_list);
-        parameter_list.set("aggregation: block scaling", true);
-      }
-    // MIS does not exist anymore, only choice are uncoupled and coupled. When using
-    // uncoupled, aggregates cannot span multiple processes. When using coupled
-    // aggregates can span multiple processes.
-    parameter_list.set("aggregation: type", "Uncoupled");
-
-    parameter_list.set("smoother: type", additional_data.smoother_type);
-    parameter_list.set("coarse: type", additional_data.coarse_type);
-
-    parameter_list.set("smoother: sweeps",
-                       static_cast<int>(additional_data.smoother_sweeps));
-    parameter_list.set("cycle applications",
-                       static_cast<int>(additional_data.n_cycles));
-    if (additional_data.w_cycle == true)
-      parameter_list.set("prec type", "MGW");
-    else
-      parameter_list.set("prec type", "MGV");
-
-    parameter_list.set("smoother: Chebyshev alpha",10.);
-    parameter_list.set("smoother: ifpack overlap",
-                       static_cast<int>(additional_data.smoother_overlap));
-    parameter_list.set("aggregation: threshold",
-                       additional_data.aggregation_threshold);
-    parameter_list.set("coarse: max size", 2000);
-
-    if (additional_data.output_details)
-      parameter_list.set("ML output", 10);
-    else
-      parameter_list.set("ML output", 0);
-
-    const Epetra_Map &domain_map = matrix.OperatorDomainMap();
-
-    const size_type constant_modes_dimension =
-      additional_data.constant_modes.size();
-    Epetra_MultiVector distributed_constant_modes (domain_map,
-                                                   constant_modes_dimension > 0 ?
-                                                   constant_modes_dimension : 1);
-    std::vector<double> dummy (constant_modes_dimension);
-
-    if (constant_modes_dimension > 0)
-      {
-        const size_type n_rows = n_global_rows(matrix);
-        const bool constant_modes_are_global =
-          additional_data.constant_modes[0].size() == n_rows;
-        const size_type n_relevant_rows =
-          constant_modes_are_global ? n_rows : additional_data.constant_modes[0].size();
-        const size_type my_size = domain_map.NumMyElements();
-        if (constant_modes_are_global == false)
-          Assert (n_relevant_rows == my_size,
-                  ExcDimensionMismatch(n_relevant_rows, my_size));
-        Assert (n_rows ==
-                static_cast<size_type>(global_length(distributed_constant_modes)),
-                ExcDimensionMismatch(n_rows,
-                                     global_length(distributed_constant_modes)));
-
-        (void)n_relevant_rows;
-        (void)global_length;
-
-        // Reshape null space as a contiguous vector of doubles so that
-        // Trilinos can read from it.
-        for (size_type d=0; d<constant_modes_dimension; ++d)
-          for (size_type row=0; row<my_size; ++row)
-            {
-              TrilinosWrappers::types::int_type global_row_id =
-                constant_modes_are_global ? gid(domain_map,row) : row;
-              distributed_constant_modes[d][row] =
-                additional_data.constant_modes[d][global_row_id];
-            }
-
-        parameter_list.set("null space: type", "pre-computed");
-        parameter_list.set("null space: dimension",
-                           distributed_constant_modes.NumVectors());
-        if (my_size > 0)
-          parameter_list.set("null space: vectors",
-                             distributed_constant_modes.Values());
-        // We need to set a valid pointer to data even if there is no data on
-        // the current processor. Therefore, pass a dummy in that case
-        else
-          parameter_list.set("null space: vectors",
-                             &dummy[0]);
-      }
-
-    initialize (matrix, parameter_list);
-  }
-
-
-
-  void
-  PreconditionAMGMueLu::initialize (const SparseMatrix           &matrix,
-                                    Teuchos::ParameterList &muelu_parameters)
-  {
-    initialize(matrix.trilinos_matrix(), muelu_parameters);
-  }
-
-
-
-  void
-  PreconditionAMGMueLu::initialize (const Epetra_CrsMatrix       &matrix,
-                                    Teuchos::ParameterList &muelu_parameters)
-  {
-    // We cannot use MueLu::CreateEpetraOperator directly because, we cannot
-    // transfer ownership of MueLu::EpetraOperator from Teuchos::RCP to
-    // std::shared_ptr.
-
-    // For now, just use serial node, i.e. no multithreaing or GPU.
-    typedef KokkosClassic::DefaultNode::DefaultNodeType node;
-    preconditioner.reset ();
-
-    // Cast matrix into a MueLu::Matrix. The constness needs to be cast away.
-    // MueLu uses Teuchos::RCP which are Trilinos version of std::shared_ptr.
-    Teuchos::RCP<Epetra_CrsMatrix> rcp_matrix = Teuchos::rcpFromRef(
-                                                  *(const_cast<Epetra_CrsMatrix *>(&matrix)));
-    Teuchos::RCP<Xpetra::CrsMatrix<double,int,int,node> > muelu_crs_matrix =
-      Teuchos::rcp(new Xpetra::EpetraCrsMatrix (rcp_matrix));
-    Teuchos::RCP<Xpetra::Matrix<double,int,int,node> > muelu_matrix =
-      Teuchos::rcp(new Xpetra::CrsMatrixWrap<double,int,int,node> (muelu_crs_matrix));
-
-    // Create the multigrid hierarchy using ML parameters.
-    Teuchos::RCP<MueLu::HierarchyManager<double,int,int,node> > hierarchy_factory;
-    hierarchy_factory = Teuchos::rcp(
-                          new MueLu::MLParameterListInterpreter<double,int,int,node> (muelu_parameters));
-    Teuchos::RCP<MueLu::Hierarchy<double,int,int,node> > hierarchy = hierarchy_factory->CreateHierarchy();
-    hierarchy->GetLevel(0)->Set("A",muelu_matrix);
-    hierarchy_factory->SetupHierarchy(*hierarchy);
-
-    // MueLu::EpetraOperator is just a wrapper around a "standard"
-    // Epetra_Operator.
-    preconditioner.reset(new MueLu::EpetraOperator(hierarchy));
-  }
-
-
-
-  template <typename number>
-  void
-  PreconditionAMGMueLu::
-  initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
-              const AdditionalData                 &additional_data,
-              const double                          drop_tolerance,
-              const ::dealii::SparsityPattern      *use_this_sparsity)
-  {
-    preconditioner.reset();
-    const size_type n_rows = deal_ii_sparse_matrix.m();
-
-    // Init Epetra Matrix using an
-    // equidistributed map; avoid
-    // storing the nonzero
-    // elements.
-    vector_distributor.reset (new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n_rows),
-                                             0, communicator));
-
-    if (trilinos_matrix.get() == 0)
-      trilinos_matrix.reset (new SparseMatrix());
-
-    trilinos_matrix->reinit (*vector_distributor, *vector_distributor,
-                             deal_ii_sparse_matrix, drop_tolerance, true,
-                             use_this_sparsity);
-
-    initialize (*trilinos_matrix, additional_data);
-  }
-
-
-
-  void PreconditionAMGMueLu::clear ()
-  {
-    PreconditionBase::clear();
-    trilinos_matrix.reset();
-  }
-
-
-
-  PreconditionAMGMueLu::size_type
-  PreconditionAMGMueLu::memory_consumption() const
-  {
-    unsigned int memory = sizeof(this);
-
-    // todo: find a way to read out ML's data
-    // sizes
-    if (trilinos_matrix.get() != 0)
-      memory += trilinos_matrix->memory_consumption();
-    return memory;
-  }
-
-
-
-
-  // explicit instantiations
-  template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix<double> &,
-                                                  const AdditionalData &, const double,
-                                                  const ::dealii::SparsityPattern *);
-  template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix<float> &,
-                                                  const AdditionalData &, const double,
-                                                  const ::dealii::SparsityPattern *);
-#endif
-
-
 
   /* -------------------------- PreconditionIdentity --------------------- */
 
diff --git a/source/lac/trilinos_precondition2.cc b/source/lac/trilinos_precondition2.cc
new file mode 100644 (file)
index 0000000..32bd5f8
--- /dev/null
@@ -0,0 +1,339 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/trilinos_precondition.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#  include <deal.II/lac/vector.h>
+#  include <deal.II/lac/sparse_matrix.h>
+#  include <deal.II/lac/trilinos_sparse_matrix.h>
+
+DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+#  include <Ifpack.h>
+#  include <Ifpack_Chebyshev.h>
+#  include <Teuchos_ParameterList.hpp>
+#  include <Teuchos_RCP.hpp>
+#  include <Epetra_MultiVector.h>
+#  include <ml_include.h>
+#  include <ml_MultiLevelPreconditioner.h>
+DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+  namespace
+  {
+#ifndef DEAL_II_WITH_64BIT_INDICES
+    int n_global_rows (const Epetra_RowMatrix &matrix)
+    {
+      return matrix.NumGlobalRows();
+    }
+
+    int global_length (const Epetra_MultiVector &vector)
+    {
+      return vector.GlobalLength();
+    }
+
+    int gid(const Epetra_Map &map, unsigned int i)
+    {
+      return map.GID(i);
+    }
+#else
+    long long int n_global_rows (const Epetra_RowMatrix &matrix)
+    {
+      return matrix.NumGlobalRows64();
+    }
+
+    long long int global_length (const Epetra_MultiVector &vector)
+    {
+      return vector.GlobalLength64();
+    }
+
+    long long int gid(const Epetra_Map &map, dealii::types::global_dof_index i)
+    {
+      return map.GID64(i);
+    }
+#endif
+  }
+
+
+
+  /* -------------------------- PreconditionAMG -------------------------- */
+
+  PreconditionAMG::AdditionalData::
+  AdditionalData (const bool                             elliptic,
+                  const bool                             higher_order_elements,
+                  const unsigned int                     n_cycles,
+                  const bool                             w_cycle,
+                  const double                           aggregation_threshold,
+                  const std::vector<std::vector<bool> > &constant_modes,
+                  const unsigned int                     smoother_sweeps,
+                  const unsigned int                     smoother_overlap,
+                  const bool                             output_details,
+                  const char                            *smoother_type,
+                  const char                            *coarse_type)
+    :
+    elliptic (elliptic),
+    higher_order_elements (higher_order_elements),
+    n_cycles (n_cycles),
+    w_cycle (w_cycle),
+    aggregation_threshold (aggregation_threshold),
+    constant_modes (constant_modes),
+    smoother_sweeps (smoother_sweeps),
+    smoother_overlap (smoother_overlap),
+    output_details (output_details),
+    smoother_type (smoother_type),
+    coarse_type (coarse_type)
+  {}
+
+
+  PreconditionAMG::~PreconditionAMG()
+  {
+    preconditioner.reset();
+    trilinos_matrix.reset();
+  }
+
+
+
+  void
+  PreconditionAMG:: initialize (const SparseMatrix   &matrix,
+                                const AdditionalData &additional_data)
+  {
+    initialize(matrix.trilinos_matrix(), additional_data);
+  }
+
+
+
+  void
+  PreconditionAMG:: initialize (const Epetra_RowMatrix &matrix,
+                                const AdditionalData   &additional_data)
+  {
+    // Build the AMG preconditioner.
+    Teuchos::ParameterList parameter_list;
+
+    if (additional_data.elliptic == true)
+      {
+        ML_Epetra::SetDefaults("SA",parameter_list);
+
+        // uncoupled mode can give a lot of warnings or even fail when there
+        // are too many entries per row and aggreggation gets complicated, but
+        // MIS does not work if too few elements are located on one
+        // processor. work around these warnings by choosing the different
+        // strategies in different situations: for low order, always use the
+        // standard choice uncoupled. if higher order, right now we also just
+        // use Uncoupled, but we should be aware that maybe MIS might be
+        // needed
+        if (additional_data.higher_order_elements)
+          parameter_list.set("aggregation: type", "Uncoupled");
+      }
+    else
+      {
+        ML_Epetra::SetDefaults("NSSA",parameter_list);
+        parameter_list.set("aggregation: type", "Uncoupled");
+        parameter_list.set("aggregation: block scaling", true);
+      }
+
+    parameter_list.set("smoother: type", additional_data.smoother_type);
+    parameter_list.set("coarse: type", additional_data.coarse_type);
+
+    parameter_list.set("smoother: sweeps",
+                       static_cast<int>(additional_data.smoother_sweeps));
+    parameter_list.set("cycle applications",
+                       static_cast<int>(additional_data.n_cycles));
+    if (additional_data.w_cycle == true)
+      parameter_list.set("prec type", "MGW");
+    else
+      parameter_list.set("prec type", "MGV");
+
+    parameter_list.set("smoother: Chebyshev alpha",10.);
+    parameter_list.set("smoother: ifpack overlap",
+                       static_cast<int>(additional_data.smoother_overlap));
+    parameter_list.set("aggregation: threshold",
+                       additional_data.aggregation_threshold);
+    parameter_list.set("coarse: max size", 2000);
+
+    if (additional_data.output_details)
+      parameter_list.set("ML output", 10);
+    else
+      parameter_list.set("ML output", 0);
+
+    const Epetra_Map &domain_map = matrix.OperatorDomainMap();
+
+    const size_type constant_modes_dimension =
+      additional_data.constant_modes.size();
+    Epetra_MultiVector distributed_constant_modes (domain_map,
+                                                   constant_modes_dimension > 0 ?
+                                                   constant_modes_dimension : 1);
+    std::vector<double> dummy (constant_modes_dimension);
+
+    if (constant_modes_dimension > 0)
+      {
+        const size_type n_rows = n_global_rows(matrix);
+        const bool constant_modes_are_global =
+          additional_data.constant_modes[0].size() == n_rows;
+        const size_type n_relevant_rows =
+          constant_modes_are_global ? n_rows : additional_data.constant_modes[0].size();
+        const size_type my_size = domain_map.NumMyElements();
+        if (constant_modes_are_global == false)
+          Assert (n_relevant_rows == my_size,
+                  ExcDimensionMismatch(n_relevant_rows, my_size));
+        Assert (n_rows ==
+                static_cast<size_type>(global_length(distributed_constant_modes)),
+                ExcDimensionMismatch(n_rows,
+                                     global_length(distributed_constant_modes)));
+
+        (void)n_relevant_rows;
+        (void)global_length;
+
+        // Reshape null space as a contiguous vector of doubles so that
+        // Trilinos can read from it.
+        for (size_type d=0; d<constant_modes_dimension; ++d)
+          for (size_type row=0; row<my_size; ++row)
+            {
+              TrilinosWrappers::types::int_type global_row_id =
+                constant_modes_are_global ? gid(domain_map,row) : row;
+              distributed_constant_modes[d][row] =
+                additional_data.constant_modes[d][global_row_id];
+            }
+
+        parameter_list.set("null space: type", "pre-computed");
+        parameter_list.set("null space: dimension",
+                           distributed_constant_modes.NumVectors());
+        if (my_size > 0)
+          parameter_list.set("null space: vectors",
+                             distributed_constant_modes.Values());
+        // We need to set a valid pointer to data even if there is no data on
+        // the current processor. Therefore, pass a dummy in that case
+        else
+          parameter_list.set("null space: vectors",
+                             &dummy[0]);
+      }
+
+    initialize (matrix, parameter_list);
+
+    if (additional_data.output_details)
+      {
+        ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
+          dynamic_cast<ML_Epetra::MultiLevelPreconditioner *> (preconditioner.get());
+        Assert (multilevel_operator != 0,
+                ExcMessage ("Preconditioner setup failed."));
+        multilevel_operator->PrintUnused(0);
+      }
+  }
+
+
+
+  void
+  PreconditionAMG::initialize (const SparseMatrix           &matrix,
+                               const Teuchos::ParameterList &ml_parameters)
+  {
+    initialize(matrix.trilinos_matrix(), ml_parameters);
+  }
+
+
+
+  void
+  PreconditionAMG::initialize (const Epetra_RowMatrix       &matrix,
+                               const Teuchos::ParameterList &ml_parameters)
+  {
+    preconditioner.reset ();
+    preconditioner.reset (new ML_Epetra::MultiLevelPreconditioner
+                          (matrix, ml_parameters));
+  }
+
+
+
+  template <typename number>
+  void
+  PreconditionAMG::
+  initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
+              const AdditionalData                 &additional_data,
+              const double                          drop_tolerance,
+              const ::dealii::SparsityPattern      *use_this_sparsity)
+  {
+    preconditioner.reset();
+    const size_type n_rows = deal_ii_sparse_matrix.m();
+
+    // Init Epetra Matrix using an
+    // equidistributed map; avoid
+    // storing the nonzero
+    // elements.
+    vector_distributor.reset (new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n_rows),
+                                             0, communicator));
+
+    if (trilinos_matrix.get() == 0)
+      trilinos_matrix.reset (new SparseMatrix());
+
+    trilinos_matrix->reinit (*vector_distributor, *vector_distributor,
+                             deal_ii_sparse_matrix, drop_tolerance, true,
+                             use_this_sparsity);
+
+    initialize (*trilinos_matrix, additional_data);
+  }
+
+
+
+  void PreconditionAMG::reinit ()
+  {
+    ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
+      dynamic_cast<ML_Epetra::MultiLevelPreconditioner *> (preconditioner.get());
+    multilevel_operator->ReComputePreconditioner();
+  }
+
+
+
+  void PreconditionAMG::clear ()
+  {
+    PreconditionBase::clear();
+    trilinos_matrix.reset();
+  }
+
+
+
+  PreconditionAMG::size_type
+  PreconditionAMG::memory_consumption() const
+  {
+    unsigned int memory = sizeof(this);
+
+    // todo: find a way to read out ML's data
+    // sizes
+    if (trilinos_matrix.get() != 0)
+      memory += trilinos_matrix->memory_consumption();
+    return memory;
+  }
+
+
+
+
+  // explicit instantiations
+  template void PreconditionAMG::initialize (const ::dealii::SparseMatrix<double> &,
+                                             const AdditionalData &, const double,
+                                             const ::dealii::SparsityPattern *);
+  template void PreconditionAMG::initialize (const ::dealii::SparseMatrix<float> &,
+                                             const AdditionalData &, const double,
+                                             const ::dealii::SparsityPattern *);
+
+
+
+
+
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_TRILINOS
diff --git a/source/lac/trilinos_precondition3.cc b/source/lac/trilinos_precondition3.cc
new file mode 100644 (file)
index 0000000..d23be4a
--- /dev/null
@@ -0,0 +1,338 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/trilinos_precondition.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#  include <deal.II/lac/vector.h>
+#  include <deal.II/lac/sparse_matrix.h>
+#  include <deal.II/lac/trilinos_sparse_matrix.h>
+
+DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+#  include <Teuchos_ParameterList.hpp>
+#  include <Teuchos_RCP.hpp>
+#  include <Epetra_MultiVector.h>
+#  include <ml_include.h>
+#  include <ml_MultiLevelPreconditioner.h>
+
+#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
+#  include <MueLu.hpp>
+#  include <MueLu_EpetraOperator.hpp>
+#  include <MueLu_MLParameterListInterpreter.hpp>
+#endif
+DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+  namespace
+  {
+#ifndef DEAL_II_WITH_64BIT_INDICES
+    int n_global_rows (const Epetra_RowMatrix &matrix)
+    {
+      return matrix.NumGlobalRows();
+    }
+
+    int global_length (const Epetra_MultiVector &vector)
+    {
+      return vector.GlobalLength();
+    }
+
+    int gid(const Epetra_Map &map, unsigned int i)
+    {
+      return map.GID(i);
+    }
+#else
+    long long int n_global_rows (const Epetra_RowMatrix &matrix)
+    {
+      return matrix.NumGlobalRows64();
+    }
+
+    long long int global_length (const Epetra_MultiVector &vector)
+    {
+      return vector.GlobalLength64();
+    }
+
+    long long int gid(const Epetra_Map &map, dealii::types::global_dof_index i)
+    {
+      return map.GID64(i);
+    }
+#endif
+  }
+
+
+
+
+
+#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
+  /* -------------------------- PreconditionAMGMueLu --------------------- */
+
+  PreconditionAMGMueLu::AdditionalData::
+  AdditionalData (const bool                             elliptic,
+                  const unsigned int                     n_cycles,
+                  const bool                             w_cycle,
+                  const double                           aggregation_threshold,
+                  const std::vector<std::vector<bool> > &constant_modes,
+                  const unsigned int                     smoother_sweeps,
+                  const unsigned int                     smoother_overlap,
+                  const bool                             output_details,
+                  const char                            *smoother_type,
+                  const char                            *coarse_type)
+    :
+    elliptic (elliptic),
+    n_cycles (n_cycles),
+    w_cycle (w_cycle),
+    aggregation_threshold (aggregation_threshold),
+    constant_modes (constant_modes),
+    smoother_sweeps (smoother_sweeps),
+    smoother_overlap (smoother_overlap),
+    output_details (output_details),
+    smoother_type (smoother_type),
+    coarse_type (coarse_type)
+  {}
+
+
+  PreconditionAMGMueLu::~PreconditionAMGMueLu()
+  {
+    preconditioner.reset();
+    trilinos_matrix.reset();
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const SparseMatrix   &matrix,
+                                    const AdditionalData &additional_data)
+  {
+    initialize(matrix.trilinos_matrix(), additional_data);
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const Epetra_CrsMatrix &matrix,
+                                    const AdditionalData   &additional_data)
+  {
+    // Build the AMG preconditioner.
+    Teuchos::ParameterList parameter_list;
+
+    if (additional_data.elliptic == true)
+      ML_Epetra::SetDefaults("SA",parameter_list);
+    else
+      {
+        ML_Epetra::SetDefaults("NSSA",parameter_list);
+        parameter_list.set("aggregation: block scaling", true);
+      }
+    // MIS does not exist anymore, only choice are uncoupled and coupled. When using
+    // uncoupled, aggregates cannot span multiple processes. When using coupled
+    // aggregates can span multiple processes.
+    parameter_list.set("aggregation: type", "Uncoupled");
+
+    parameter_list.set("smoother: type", additional_data.smoother_type);
+    parameter_list.set("coarse: type", additional_data.coarse_type);
+
+    parameter_list.set("smoother: sweeps",
+                       static_cast<int>(additional_data.smoother_sweeps));
+    parameter_list.set("cycle applications",
+                       static_cast<int>(additional_data.n_cycles));
+    if (additional_data.w_cycle == true)
+      parameter_list.set("prec type", "MGW");
+    else
+      parameter_list.set("prec type", "MGV");
+
+    parameter_list.set("smoother: Chebyshev alpha",10.);
+    parameter_list.set("smoother: ifpack overlap",
+                       static_cast<int>(additional_data.smoother_overlap));
+    parameter_list.set("aggregation: threshold",
+                       additional_data.aggregation_threshold);
+    parameter_list.set("coarse: max size", 2000);
+
+    if (additional_data.output_details)
+      parameter_list.set("ML output", 10);
+    else
+      parameter_list.set("ML output", 0);
+
+    const Epetra_Map &domain_map = matrix.OperatorDomainMap();
+
+    const size_type constant_modes_dimension =
+      additional_data.constant_modes.size();
+    Epetra_MultiVector distributed_constant_modes (domain_map,
+                                                   constant_modes_dimension > 0 ?
+                                                   constant_modes_dimension : 1);
+    std::vector<double> dummy (constant_modes_dimension);
+
+    if (constant_modes_dimension > 0)
+      {
+        const size_type n_rows = n_global_rows(matrix);
+        const bool constant_modes_are_global =
+          additional_data.constant_modes[0].size() == n_rows;
+        const size_type n_relevant_rows =
+          constant_modes_are_global ? n_rows : additional_data.constant_modes[0].size();
+        const size_type my_size = domain_map.NumMyElements();
+        if (constant_modes_are_global == false)
+          Assert (n_relevant_rows == my_size,
+                  ExcDimensionMismatch(n_relevant_rows, my_size));
+        Assert (n_rows ==
+                static_cast<size_type>(global_length(distributed_constant_modes)),
+                ExcDimensionMismatch(n_rows,
+                                     global_length(distributed_constant_modes)));
+
+        (void)n_relevant_rows;
+        (void)global_length;
+
+        // Reshape null space as a contiguous vector of doubles so that
+        // Trilinos can read from it.
+        for (size_type d=0; d<constant_modes_dimension; ++d)
+          for (size_type row=0; row<my_size; ++row)
+            {
+              TrilinosWrappers::types::int_type global_row_id =
+                constant_modes_are_global ? gid(domain_map,row) : row;
+              distributed_constant_modes[d][row] =
+                additional_data.constant_modes[d][global_row_id];
+            }
+
+        parameter_list.set("null space: type", "pre-computed");
+        parameter_list.set("null space: dimension",
+                           distributed_constant_modes.NumVectors());
+        if (my_size > 0)
+          parameter_list.set("null space: vectors",
+                             distributed_constant_modes.Values());
+        // We need to set a valid pointer to data even if there is no data on
+        // the current processor. Therefore, pass a dummy in that case
+        else
+          parameter_list.set("null space: vectors",
+                             &dummy[0]);
+      }
+
+    initialize (matrix, parameter_list);
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const SparseMatrix           &matrix,
+                                    Teuchos::ParameterList &muelu_parameters)
+  {
+    initialize(matrix.trilinos_matrix(), muelu_parameters);
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const Epetra_CrsMatrix       &matrix,
+                                    Teuchos::ParameterList &muelu_parameters)
+  {
+    // We cannot use MueLu::CreateEpetraOperator directly because, we cannot
+    // transfer ownership of MueLu::EpetraOperator from Teuchos::RCP to
+    // std::shared_ptr.
+
+    // For now, just use serial node, i.e. no multithreaing or GPU.
+    typedef KokkosClassic::DefaultNode::DefaultNodeType node;
+    preconditioner.reset ();
+
+    // Cast matrix into a MueLu::Matrix. The constness needs to be cast away.
+    // MueLu uses Teuchos::RCP which are Trilinos version of std::shared_ptr.
+    Teuchos::RCP<Epetra_CrsMatrix> rcp_matrix = Teuchos::rcpFromRef(
+                                                  *(const_cast<Epetra_CrsMatrix *>(&matrix)));
+    Teuchos::RCP<Xpetra::CrsMatrix<double,int,int,node> > muelu_crs_matrix =
+      Teuchos::rcp(new Xpetra::EpetraCrsMatrix (rcp_matrix));
+    Teuchos::RCP<Xpetra::Matrix<double,int,int,node> > muelu_matrix =
+      Teuchos::rcp(new Xpetra::CrsMatrixWrap<double,int,int,node> (muelu_crs_matrix));
+
+    // Create the multigrid hierarchy using ML parameters.
+    Teuchos::RCP<MueLu::HierarchyManager<double,int,int,node> > hierarchy_factory;
+    hierarchy_factory = Teuchos::rcp(
+                          new MueLu::MLParameterListInterpreter<double,int,int,node> (muelu_parameters));
+    Teuchos::RCP<MueLu::Hierarchy<double,int,int,node> > hierarchy = hierarchy_factory->CreateHierarchy();
+    hierarchy->GetLevel(0)->Set("A",muelu_matrix);
+    hierarchy_factory->SetupHierarchy(*hierarchy);
+
+    // MueLu::EpetraOperator is just a wrapper around a "standard"
+    // Epetra_Operator.
+    preconditioner.reset(new MueLu::EpetraOperator(hierarchy));
+  }
+
+
+
+  template <typename number>
+  void
+  PreconditionAMGMueLu::
+  initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
+              const AdditionalData                 &additional_data,
+              const double                          drop_tolerance,
+              const ::dealii::SparsityPattern      *use_this_sparsity)
+  {
+    preconditioner.reset();
+    const size_type n_rows = deal_ii_sparse_matrix.m();
+
+    // Init Epetra Matrix using an
+    // equidistributed map; avoid
+    // storing the nonzero
+    // elements.
+    vector_distributor.reset (new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n_rows),
+                                             0, communicator));
+
+    if (trilinos_matrix.get() == 0)
+      trilinos_matrix.reset (new SparseMatrix());
+
+    trilinos_matrix->reinit (*vector_distributor, *vector_distributor,
+                             deal_ii_sparse_matrix, drop_tolerance, true,
+                             use_this_sparsity);
+
+    initialize (*trilinos_matrix, additional_data);
+  }
+
+
+
+  void PreconditionAMGMueLu::clear ()
+  {
+    PreconditionBase::clear();
+    trilinos_matrix.reset();
+  }
+
+
+
+  PreconditionAMGMueLu::size_type
+  PreconditionAMGMueLu::memory_consumption() const
+  {
+    unsigned int memory = sizeof(this);
+
+    // todo: find a way to read out ML's data
+    // sizes
+    if (trilinos_matrix.get() != 0)
+      memory += trilinos_matrix->memory_consumption();
+    return memory;
+  }
+
+
+
+  // explicit instantiations
+  template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix<double> &,
+                                                  const AdditionalData &, const double,
+                                                  const ::dealii::SparsityPattern *);
+  template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix<float> &,
+                                                  const AdditionalData &, const double,
+                                                  const ::dealii::SparsityPattern *);
+#endif
+
+
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_TRILINOS
index 17ab06352ed1113daf756fd938d76a14aadf8cf1..086368040c8ab7cb027056b467d8ba6c73a5d3f1 100644 (file)
@@ -26,6 +26,7 @@ SET(_src
   dof_output_operator.cc
   error_estimator_1d.cc
   error_estimator.cc
+  error_estimator_inst2.cc
   fe_field_function.cc
   histogram.cc
   matrix_tools.cc
index 93383b29cda9cdd22fd91349e21ab3ec33e3326c..2b40c679956af52f15c78cde4075bf97b35c5d83 100644 (file)
@@ -1271,7 +1271,10 @@ void KellyErrorEstimator<dim, spacedim>::estimate (const DH
 
 
 
-// explicit instantiations
+#define SPLIT_INSTANTIATIONS_COUNT 2
+#ifndef SPLIT_INSTANTIATIONS_INDEX
+#define SPLIT_INSTANTIATIONS_INDEX 0
+#endif
 #include "error_estimator.inst"
 
 
diff --git a/source/numerics/error_estimator_inst2.cc b/source/numerics/error_estimator_inst2.cc
new file mode 100644 (file)
index 0000000..07f99ac
--- /dev/null
@@ -0,0 +1,20 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+// This file compiles the second half of the instantiations from error_estimator.cc
+// to get the memory consumption below 1.5gb with gcc.
+
+#define SPLIT_INSTANTIATIONS_INDEX 1
+#include "error_estimator.cc"

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.