From: Timo Heister Date: Sat, 5 Sep 2015 15:36:25 +0000 (-0400) Subject: split compilation units X-Git-Tag: v8.4.0-rc2~470^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=69f75a491cf534901254957e1c604d8cc8344808;p=dealii.git split compilation units --- diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 904e10d5cf..faa37b8b6e 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -4051,12 +4051,9 @@ void FESubfaceValues::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" diff --git a/source/fe/fe_values.inst.in b/source/fe/fe_values.inst.in index 04d8d1400b..8678b601a1 100644 --- a/source/fe/fe_values.inst.in +++ b/source/fe/fe_values.inst.in @@ -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; template class FEValues; template class FEValuesBase:: @@ -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::reinit( const TriaIterator, lda> >&); template void FEFaceValues::reinit( const TriaIterator, lda> >&, unsigned int); template void FESubfaceValues::reinit( const TriaIterator, 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 ::get_function_values @@ -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 (const dealii::VEC&, std::vector >::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::get_function_values @@ -214,7 +205,6 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi (const VEC&, const VectorSlice >&, std::vector > &, 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::get_function_values (const dealii::IndexSet&, std::vector::type> &) const; @@ -335,7 +324,6 @@ for (deal_II_dimension : DIMENSIONS) ::get_function_divergences (const dealii::IndexSet&, std::vector >::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::get_function_values (const IndexSet&, std::vector&) const; @@ -498,6 +485,4 @@ for (deal_II_dimension : DIMENSIONS) std::vector > &, bool) const; #endif -#endif - } diff --git a/source/fe/fe_values_inst2.cc b/source/fe/fe_values_inst2.cc index ad1aa2b829..9519d43e7d 100644 --- a/source/fe/fe_values_inst2.cc +++ b/source/fe/fe_values_inst2.cc @@ -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" diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index e5bfce3c0d..d9decfb112 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -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 diff --git a/source/lac/sparse_matrix.cc b/source/lac/sparse_matrix.cc index ea4ce40002..4db5f38645 100644 --- a/source/lac/sparse_matrix.cc +++ b/source/lac/sparse_matrix.cc @@ -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 index 0000000000..aaa7275f92 --- /dev/null +++ b/source/lac/sparse_matrix_inst2.cc @@ -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 +#include + +DEAL_II_NAMESPACE_OPEN + +#define SPLIT_INSTANTIATIONS_COUNT 2 +#define SPLIT_INSTANTIATIONS_INDEX 1 +#include "sparse_matrix.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/trilinos_precondition.cc b/source/lac/trilinos_precondition.cc index e7582f617e..9c644c439a 100644 --- a/source/lac/trilinos_precondition.cc +++ b/source/lac/trilinos_precondition.cc @@ -27,14 +27,6 @@ DEAL_II_DISABLE_EXTRA_DIAGNOSTICS # include # include # include -# include -# include - -#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0) -# include -# include -# include -#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 > &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(additional_data.smoother_sweeps)); - parameter_list.set("cycle applications", - static_cast(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(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 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(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 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 (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 - void - PreconditionAMG:: - initialize (const ::dealii::SparseMatrix &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(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 (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 &, - const AdditionalData &, const double, - const ::dealii::SparsityPattern *); - template void PreconditionAMG::initialize (const ::dealii::SparseMatrix &, - 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 > &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(additional_data.smoother_sweeps)); - parameter_list.set("cycle applications", - static_cast(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(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 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(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 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 rcp_matrix = Teuchos::rcpFromRef( - *(const_cast(&matrix))); - Teuchos::RCP > muelu_crs_matrix = - Teuchos::rcp(new Xpetra::EpetraCrsMatrix (rcp_matrix)); - Teuchos::RCP > muelu_matrix = - Teuchos::rcp(new Xpetra::CrsMatrixWrap (muelu_crs_matrix)); - - // Create the multigrid hierarchy using ML parameters. - Teuchos::RCP > hierarchy_factory; - hierarchy_factory = Teuchos::rcp( - new MueLu::MLParameterListInterpreter (muelu_parameters)); - Teuchos::RCP > 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 - void - PreconditionAMGMueLu:: - initialize (const ::dealii::SparseMatrix &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(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 &, - const AdditionalData &, const double, - const ::dealii::SparsityPattern *); - template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix &, - 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 index 0000000000..32bd5f82b4 --- /dev/null +++ b/source/lac/trilinos_precondition2.cc @@ -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 + +#ifdef DEAL_II_WITH_TRILINOS + +# include +# include +# include + +DEAL_II_DISABLE_EXTRA_DIAGNOSTICS +# include +# include +# include +# include +# include +# include +# include +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 > &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(additional_data.smoother_sweeps)); + parameter_list.set("cycle applications", + static_cast(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(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 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(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 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 (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 + void + PreconditionAMG:: + initialize (const ::dealii::SparseMatrix &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(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 (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 &, + const AdditionalData &, const double, + const ::dealii::SparsityPattern *); + template void PreconditionAMG::initialize (const ::dealii::SparseMatrix &, + 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 index 0000000000..d23be4a75b --- /dev/null +++ b/source/lac/trilinos_precondition3.cc @@ -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 + +#ifdef DEAL_II_WITH_TRILINOS + +# include +# include +# include + +DEAL_II_DISABLE_EXTRA_DIAGNOSTICS +# include +# include +# include +# include +# include + +#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0) +# include +# include +# include +#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 > &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(additional_data.smoother_sweeps)); + parameter_list.set("cycle applications", + static_cast(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(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 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(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 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 rcp_matrix = Teuchos::rcpFromRef( + *(const_cast(&matrix))); + Teuchos::RCP > muelu_crs_matrix = + Teuchos::rcp(new Xpetra::EpetraCrsMatrix (rcp_matrix)); + Teuchos::RCP > muelu_matrix = + Teuchos::rcp(new Xpetra::CrsMatrixWrap (muelu_crs_matrix)); + + // Create the multigrid hierarchy using ML parameters. + Teuchos::RCP > hierarchy_factory; + hierarchy_factory = Teuchos::rcp( + new MueLu::MLParameterListInterpreter (muelu_parameters)); + Teuchos::RCP > 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 + void + PreconditionAMGMueLu:: + initialize (const ::dealii::SparseMatrix &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(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 &, + const AdditionalData &, const double, + const ::dealii::SparsityPattern *); + template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix &, + const AdditionalData &, const double, + const ::dealii::SparsityPattern *); +#endif + + + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_TRILINOS diff --git a/source/numerics/CMakeLists.txt b/source/numerics/CMakeLists.txt index 17ab06352e..086368040c 100644 --- a/source/numerics/CMakeLists.txt +++ b/source/numerics/CMakeLists.txt @@ -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 diff --git a/source/numerics/error_estimator.cc b/source/numerics/error_estimator.cc index 93383b29cd..2b40c67995 100644 --- a/source/numerics/error_estimator.cc +++ b/source/numerics/error_estimator.cc @@ -1271,7 +1271,10 @@ void KellyErrorEstimator::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 index 0000000000..07f99ac5f2 --- /dev/null +++ b/source/numerics/error_estimator_inst2.cc @@ -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"