From d012d2a55b30aab82c2991fcfd13cdb06707fe8b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 13 Oct 2021 20:14:29 -0600 Subject: [PATCH] Implement functions within instead of outside a namespace. --- include/deal.II/sundials/n_vector.templates.h | 1070 +++++++++-------- 1 file changed, 539 insertions(+), 531 deletions(-) diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index 338d13e17c..f48798e415 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -296,697 +296,705 @@ namespace SUNDIALS -template -SUNDIALS::internal::NVectorContent::NVectorContent() - : vector(typename VectorMemory::Pointer(mem)) - , is_const(false) -{} +// -------------- Implementation of the functions declared above --------- +namespace SUNDIALS +{ + namespace internal + { + template + NVectorContent::NVectorContent() + : vector(typename VectorMemory::Pointer(mem)) + , is_const(false) + {} -template -SUNDIALS::internal::NVectorContent::NVectorContent( - VectorType *vector) - : vector(vector, [](VectorType *) { /* not owning memory -> don't free*/ }) - , is_const(false) -{} + template + NVectorContent::NVectorContent(VectorType *vector) + : vector(vector, + [](VectorType *) { /* not owning memory -> don't free*/ }) + , is_const(false) + {} -template -SUNDIALS::internal::NVectorContent::NVectorContent( - const VectorType *vector) - : vector(const_cast(vector), - [](VectorType *) { /* not owning memory -> don't free*/ }) - , is_const(true) -{} + template + NVectorContent::NVectorContent(const VectorType *vector) + : vector(const_cast(vector), + [](VectorType *) { /* not owning memory -> don't free*/ }) + , is_const(true) + {} -template -VectorType * -SUNDIALS::internal::NVectorContent::get() -{ - AssertThrow( - !is_const, - ExcMessage("Tried to access a constant vector content as non-const." - " This most likely happened because a vector that is passed to a" - " NVectorView() call should not be const.\n" - "Alternatively, if you tried to access the vector, use" - " unwrap_nvector_const() for this vector.")); - return vector.get(); -} + template + VectorType * + NVectorContent::get() + { + AssertThrow( + !is_const, + ExcMessage( + "Tried to access a constant vector content as non-const." + " This most likely happened because a vector that is passed to a" + " NVectorView() call should not be const.\n" + "Alternatively, if you tried to access the vector, use" + " unwrap_nvector_const() for this vector.")); + return vector.get(); + } -template -const VectorType * -SUNDIALS::internal::NVectorContent::get() const -{ - return vector.get(); -} + template + const VectorType * + NVectorContent::get() const + { + return vector.get(); + } -template -SUNDIALS::internal::NVectorView -SUNDIALS::internal::make_nvector_view(VectorType &vector) -{ - return NVectorView(vector); -} + template + NVectorView + make_nvector_view(VectorType &vector) + { + return NVectorView(vector); + } -template -VectorType * -SUNDIALS::internal::unwrap_nvector(N_Vector v) -{ - Assert(v != nullptr, ExcInternalError()); - Assert(v->content != nullptr, ExcInternalError()); - auto *pContent = reinterpret_cast *>(v->content); - return pContent->get(); -} + template + VectorType * + unwrap_nvector(N_Vector v) + { + Assert(v != nullptr, ExcInternalError()); + Assert(v->content != nullptr, ExcInternalError()); + auto *pContent = + reinterpret_cast *>(v->content); + return pContent->get(); + } -template -const VectorType * -SUNDIALS::internal::unwrap_nvector_const(N_Vector v) -{ - Assert(v != nullptr, ExcInternalError()); - Assert(v->content != nullptr, ExcInternalError()); - const auto *pContent = - reinterpret_cast *>(v->content); - return pContent->get(); -} + template + const VectorType * + unwrap_nvector_const(N_Vector v) + { + Assert(v != nullptr, ExcInternalError()); + Assert(v->content != nullptr, ExcInternalError()); + const auto *pContent = + reinterpret_cast *>(v->content); + return pContent->get(); + } -template -SUNDIALS::internal::NVectorView::NVectorView(VectorType &vector) - : vector_ptr( - create_nvector( - new NVectorContent::type>( - &vector)), - [](N_Vector v) { N_VDestroy(v); }) -{} + template + NVectorView::NVectorView(VectorType &vector) + : vector_ptr( + create_nvector( + new NVectorContent::type>( + &vector)), + [](N_Vector v) { N_VDestroy(v); }) + {} -template -SUNDIALS::internal::NVectorView::operator N_Vector() const -{ - Assert(vector_ptr != nullptr, ExcNotInitialized()); - return vector_ptr.get(); -} + template + NVectorView::operator N_Vector() const + { + Assert(vector_ptr != nullptr, ExcNotInitialized()); + return vector_ptr.get(); + } -template -N_Vector -SUNDIALS::internal::NVectorView::operator->() const -{ - Assert(vector_ptr != nullptr, ExcNotInitialized()); - return vector_ptr.get(); -} + template + N_Vector + NVectorView::operator->() const + { + Assert(vector_ptr != nullptr, ExcNotInitialized()); + return vector_ptr.get(); + } -template -N_Vector -SUNDIALS::internal::create_nvector(NVectorContent *content) -{ - // Create an N_Vector with operators attached and empty content - N_Vector v = create_empty_nvector(); - Assert(v != nullptr, ExcInternalError()); - v->content = content; - Assert(v->content != nullptr, ExcInternalError()); - return (v); -} + template + N_Vector + create_nvector(NVectorContent *content) + { + // Create an N_Vector with operators attached and empty content + N_Vector v = create_empty_nvector(); + Assert(v != nullptr, ExcInternalError()); + v->content = content; + Assert(v->content != nullptr, ExcInternalError()); + return (v); + } -N_Vector_ID SUNDIALS::internal::NVectorOperations::get_vector_id(N_Vector) -{ - return SUNDIALS_NVEC_CUSTOM; -} + namespace NVectorOperations + { + N_Vector_ID get_vector_id(N_Vector) + { + return SUNDIALS_NVEC_CUSTOM; + } -N_Vector -SUNDIALS::internal::NVectorOperations::clone_empty(N_Vector w) -{ - Assert(w != nullptr, ExcInternalError()); - N_Vector v = N_VNewEmpty(); - Assert(v != nullptr, ExcInternalError()); - int status = N_VCopyOps(w, v); - Assert(status == 0, ExcInternalError()); - (void)status; + N_Vector + clone_empty(N_Vector w) + { + Assert(w != nullptr, ExcInternalError()); + N_Vector v = N_VNewEmpty(); + Assert(v != nullptr, ExcInternalError()); - return v; -} + int status = N_VCopyOps(w, v); + Assert(status == 0, ExcInternalError()); + (void)status; + return v; + } -template -N_Vector -SUNDIALS::internal::NVectorOperations::clone(N_Vector w) -{ - N_Vector v = clone_empty(w); - // the corresponding delete is called in destroy() - auto cloned = new NVectorContent(); - auto *w_dealii = unwrap_nvector_const(w); + template + N_Vector + clone(N_Vector w) + { + N_Vector v = clone_empty(w); - // reinit the cloned vector based on the layout of the source vector - cloned->get()->reinit(*w_dealii); - v->content = cloned; + // the corresponding delete is called in destroy() + auto cloned = new NVectorContent(); + auto *w_dealii = unwrap_nvector_const(w); - return v; -} + // reinit the cloned vector based on the layout of the source vector + cloned->get()->reinit(*w_dealii); + v->content = cloned; + return v; + } -template -void -SUNDIALS::internal::NVectorOperations::destroy(N_Vector v) -{ - // support destroying a nullptr because SUNDIALS vectors also do it - if (v == nullptr) - return; - if (v->content != nullptr) - { - auto *content = - reinterpret_cast *>(v->content); - // the NVectorContent knows if it owns the memory and will free - // correctly - delete content; - v->content = nullptr; - } + template + void + destroy(N_Vector v) + { + // support destroying a nullptr because SUNDIALS vectors also do it + if (v == nullptr) + return; - N_VFreeEmpty(v); -} + if (v->content != nullptr) + { + auto *content = + reinterpret_cast *>(v->content); + // the NVectorContent knows if it owns the memory and will free + // correctly + delete content; + v->content = nullptr; + } + N_VFreeEmpty(v); + } -template ::value, int>> -MPI_Comm -SUNDIALS::internal::NVectorOperations::get_communicator(N_Vector v) -{ - return unwrap_nvector_const(v)->block(0).get_mpi_communicator(); -} + template ::value, int>> + MPI_Comm + get_communicator(N_Vector v) + { + return unwrap_nvector_const(v) + ->block(0) + .get_mpi_communicator(); + } -template ::value, int>> -MPI_Comm -SUNDIALS::internal::NVectorOperations::get_communicator(N_Vector v) -{ - return unwrap_nvector_const(v)->get_mpi_communicator(); -} + template ::value, int>> + MPI_Comm + get_communicator(N_Vector v) + { + return unwrap_nvector_const(v)->get_mpi_communicator(); + } -template ::value, int>> -void * - SUNDIALS::internal::NVectorOperations::get_communicator_as_void_ptr(N_Vector) -{ - // required by SUNDIALS: MPI-unaware vectors should return the nullptr as comm - return nullptr; -} + template < + typename VectorType, + typename std::enable_if_t::value, int>> + void *get_communicator_as_void_ptr(N_Vector) + { + // required by SUNDIALS: MPI-unaware vectors should return the nullptr + // as comm + return nullptr; + } -template ::value, int>> -void * -SUNDIALS::internal::NVectorOperations::get_communicator_as_void_ptr(N_Vector v) -{ + + template < + typename VectorType, + typename std::enable_if_t::value, int>> + void * + get_communicator_as_void_ptr(N_Vector v) + { # ifndef DEAL_II_WITH_MPI - (void)v; - return nullptr; + (void)v; + return nullptr; # else - return get_communicator(v); + return get_communicator(v); # endif -} + } -template -sunindextype -SUNDIALS::internal::NVectorOperations::get_global_length(N_Vector v) -{ - return unwrap_nvector_const(v)->size(); -} + template + sunindextype + get_global_length(N_Vector v) + { + return unwrap_nvector_const(v)->size(); + } -template -void -SUNDIALS::internal::NVectorOperations::linear_sum(realtype a, - N_Vector x, - realtype b, - N_Vector y, - N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *y_dealii = unwrap_nvector_const(y); - auto *z_dealii = unwrap_nvector(z); - - if (z_dealii == x_dealii) - z_dealii->sadd(a, b, *y_dealii); - else if (z_dealii == y_dealii) - z_dealii->sadd(b, a, *x_dealii); - else - { - *z_dealii = 0; - z_dealii->add(a, *x_dealii, b, *y_dealii); - } -} + template + void + linear_sum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *y_dealii = unwrap_nvector_const(y); + auto *z_dealii = unwrap_nvector(z); + if (z_dealii == x_dealii) + z_dealii->sadd(a, b, *y_dealii); + else if (z_dealii == y_dealii) + z_dealii->sadd(b, a, *x_dealii); + else + { + *z_dealii = 0; + z_dealii->add(a, *x_dealii, b, *y_dealii); + } + } -template -realtype -SUNDIALS::internal::NVectorOperations::dot_product(N_Vector x, N_Vector y) -{ - return *unwrap_nvector_const(x) * - *unwrap_nvector_const(y); -} + template + realtype + dot_product(N_Vector x, N_Vector y) + { + return *unwrap_nvector_const(x) * + *unwrap_nvector_const(y); + } -template -realtype -SUNDIALS::internal::NVectorOperations::weighted_l2_norm(N_Vector x, N_Vector w) -{ - // TODO copy can be avoided by a custom kernel - VectorType tmp = *unwrap_nvector_const(x); - auto * w_dealii = unwrap_nvector_const(w); - tmp.scale(*w_dealii); - return tmp.l2_norm(); -} + template + realtype + weighted_l2_norm(N_Vector x, N_Vector w) + { + // TODO copy can be avoided by a custom kernel + VectorType tmp = *unwrap_nvector_const(x); + auto * w_dealii = unwrap_nvector_const(w); + tmp.scale(*w_dealii); + return tmp.l2_norm(); + } -template -realtype -SUNDIALS::internal::NVectorOperations::l1_norm(N_Vector x) -{ - return unwrap_nvector_const(x)->l1_norm(); -} + template + realtype + l1_norm(N_Vector x) + { + return unwrap_nvector_const(x)->l1_norm(); + } -template -void -SUNDIALS::internal::NVectorOperations::set_constant(realtype c, N_Vector v) -{ - auto *v_dealii = unwrap_nvector(v); - *v_dealii = c; -} + template + void + set_constant(realtype c, N_Vector v) + { + auto *v_dealii = unwrap_nvector(v); + *v_dealii = c; + } -template -void -SUNDIALS::internal::NVectorOperations::add_constant(N_Vector x, - realtype b, - N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *z_dealii = unwrap_nvector(z); - if (x_dealii == z_dealii) - z_dealii->add(b); - else - { - *z_dealii = *x_dealii; - z_dealii->add(b); - } -} + template + void + add_constant(N_Vector x, realtype b, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *z_dealii = unwrap_nvector(z); + if (x_dealii == z_dealii) + z_dealii->add(b); + else + { + *z_dealii = *x_dealii; + z_dealii->add(b); + } + } -template -void -SUNDIALS::internal::NVectorOperations::elementwise_product(N_Vector x, - N_Vector y, - N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *y_dealii = unwrap_nvector_const(y); - auto *z_dealii = unwrap_nvector(z); - if (z_dealii == x_dealii) - z_dealii->scale(*y_dealii); - else if (z_dealii == y_dealii) - z_dealii->scale(*x_dealii); - else - { - *z_dealii = *y_dealii; - z_dealii->scale(*x_dealii); - } -} - + template + void + elementwise_product(N_Vector x, N_Vector y, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *y_dealii = unwrap_nvector_const(y); + auto *z_dealii = unwrap_nvector(z); + if (z_dealii == x_dealii) + z_dealii->scale(*y_dealii); + else if (z_dealii == y_dealii) + z_dealii->scale(*x_dealii); + else + { + *z_dealii = *y_dealii; + z_dealii->scale(*x_dealii); + } + } -template -realtype -SUNDIALS::internal::NVectorOperations::weighted_rms_norm(N_Vector x, N_Vector w) -{ - // TODO copy can be avoided by a custom kernel - VectorType tmp = *unwrap_nvector_const(x); - auto * w_dealii = unwrap_nvector_const(w); - const auto n = tmp.size(); - tmp.scale(*w_dealii); - return tmp.l2_norm() / std::sqrt(n); -} + template + realtype + weighted_rms_norm(N_Vector x, N_Vector w) + { + // TODO copy can be avoided by a custom kernel + VectorType tmp = *unwrap_nvector_const(x); + auto * w_dealii = unwrap_nvector_const(w); + const auto n = tmp.size(); + tmp.scale(*w_dealii); + return tmp.l2_norm() / std::sqrt(n); + } -template -realtype -SUNDIALS::internal::NVectorOperations::weighted_rms_norm_mask(N_Vector x, - N_Vector w, - N_Vector mask) -{ - // TODO copy can be avoided by a custom kernel - VectorType tmp = *unwrap_nvector_const(x); - auto * w_dealii = unwrap_nvector_const(w); - auto * mask_dealii = unwrap_nvector_const(mask); - const auto n = tmp.size(); - tmp.scale(*w_dealii); - tmp.scale(*mask_dealii); - return tmp.l2_norm() / std::sqrt(n); -} + template + realtype + weighted_rms_norm_mask(N_Vector x, N_Vector w, N_Vector mask) + { + // TODO copy can be avoided by a custom kernel + VectorType tmp = *unwrap_nvector_const(x); + auto * w_dealii = unwrap_nvector_const(w); + auto * mask_dealii = unwrap_nvector_const(mask); + const auto n = tmp.size(); + tmp.scale(*w_dealii); + tmp.scale(*mask_dealii); + return tmp.l2_norm() / std::sqrt(n); + } -template -realtype -SUNDIALS::internal::NVectorOperations::max_norm(N_Vector x) -{ - return unwrap_nvector_const(x)->linfty_norm(); -} + template + realtype + max_norm(N_Vector x) + { + return unwrap_nvector_const(x)->linfty_norm(); + } -template ::value, int>> -realtype -SUNDIALS::internal::NVectorOperations::min_element(N_Vector x) -{ - auto *vector = unwrap_nvector_const(x); - return *std::min_element(vector->begin(), vector->end()); -} + template < + typename VectorType, + typename std::enable_if_t::value, int>> + realtype + min_element(N_Vector x) + { + auto *vector = unwrap_nvector_const(x); + return *std::min_element(vector->begin(), vector->end()); + } -template ::value && - !IsBlockVector::value, - int>> -realtype -SUNDIALS::internal::NVectorOperations::min_element(N_Vector x) -{ - auto *vector = unwrap_nvector_const(x); - const auto indexed_less_than = [&](const IndexSet::size_type idxa, - const IndexSet::size_type idxb) { - return (*vector)[idxa] < (*vector)[idxb]; - }; + template < + typename VectorType, + typename std::enable_if_t::value && + !IsBlockVector::value, + int>> + realtype + min_element(N_Vector x) + { + auto *vector = unwrap_nvector_const(x); - auto local_elements = vector->locally_owned_elements(); - const auto local_min = *std::min_element(local_elements.begin(), - local_elements.end(), - indexed_less_than); - return Utilities::MPI::min((*vector)[local_min], - get_communicator(x)); -} + const auto indexed_less_than = [&](const IndexSet::size_type idxa, + const IndexSet::size_type idxb) { + return (*vector)[idxa] < (*vector)[idxb]; + }; + auto local_elements = vector->locally_owned_elements(); + const auto local_min = *std::min_element(local_elements.begin(), + local_elements.end(), + indexed_less_than); + return Utilities::MPI::min((*vector)[local_min], + get_communicator(x)); + } -template ::value && - IsBlockVector::value, - int>> -realtype -SUNDIALS::internal::NVectorOperations::min_element(N_Vector x) -{ - auto *vector = unwrap_nvector_const(x); - // initialize local minimum to the largest possible value - auto proc_local_min = - std::numeric_limits::max(); - for (unsigned i = 0; i < vector->n_blocks(); ++i) - { - const auto indexed_less_than = [&](const IndexSet::size_type idxa, - const IndexSet::size_type idxb) { - return vector->block(i)[idxa] < vector->block(i)[idxb]; - }; - - auto local_elements = vector->block(i).locally_owned_elements(); - - const auto block_local_min_element = - std::min_element(local_elements.begin(), - local_elements.end(), - indexed_less_than); - - // guard against empty blocks on this processor - if (block_local_min_element != local_elements.end()) - proc_local_min = - std::min(proc_local_min, vector->block(i)[*block_local_min_element]); - } + template < + typename VectorType, + typename std::enable_if_t::value && + IsBlockVector::value, + int>> + realtype + min_element(N_Vector x) + { + auto *vector = unwrap_nvector_const(x); - return Utilities::MPI::min(proc_local_min, get_communicator(x)); -} + // initialize local minimum to the largest possible value + auto proc_local_min = + std::numeric_limits::max(); + for (unsigned i = 0; i < vector->n_blocks(); ++i) + { + const auto indexed_less_than = [&](const IndexSet::size_type idxa, + const IndexSet::size_type idxb) { + return vector->block(i)[idxa] < vector->block(i)[idxb]; + }; + auto local_elements = vector->block(i).locally_owned_elements(); -template -void -SUNDIALS::internal::NVectorOperations::scale(realtype c, N_Vector x, N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *z_dealii = unwrap_nvector(z); + const auto block_local_min_element = + std::min_element(local_elements.begin(), + local_elements.end(), + indexed_less_than); - if (x_dealii == z_dealii) - (*z_dealii) *= c; - else - z_dealii->sadd(0.0, c, *x_dealii); -} + // guard against empty blocks on this processor + if (block_local_min_element != local_elements.end()) + proc_local_min = + std::min(proc_local_min, + vector->block(i)[*block_local_min_element]); + } + return Utilities::MPI::min(proc_local_min, + get_communicator(x)); + } -template ::value, int>> -void -SUNDIALS::internal::NVectorOperations::elementwise_div(N_Vector x, - N_Vector y, - N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *y_dealii = unwrap_nvector_const(y); - auto *z_dealii = unwrap_nvector(z); - AssertDimension(x_dealii->size(), z_dealii->size()); - AssertDimension(x_dealii->size(), y_dealii->size()); + template + void + scale(realtype c, N_Vector x, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *z_dealii = unwrap_nvector(z); - auto x_ele = x_dealii->locally_owned_elements(); - for (const auto idx : x_ele) - { - (*z_dealii)[idx] = (*x_dealii)[idx] / (*y_dealii)[idx]; - } - z_dealii->compress(VectorOperation::insert); -} + if (x_dealii == z_dealii) + (*z_dealii) *= c; + else + z_dealii->sadd(0.0, c, *x_dealii); + } -template ::value, int>> -void -SUNDIALS::internal::NVectorOperations::elementwise_div(N_Vector x, - N_Vector y, - N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *y_dealii = unwrap_nvector_const(y); - auto *z_dealii = unwrap_nvector(z); + template ::value, int>> + void + elementwise_div(N_Vector x, N_Vector y, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *y_dealii = unwrap_nvector_const(y); + auto *z_dealii = unwrap_nvector(z); - AssertDimension(x_dealii->size(), z_dealii->size()); - AssertDimension(x_dealii->size(), y_dealii->size()); - AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks()); - AssertDimension(x_dealii->n_blocks(), y_dealii->n_blocks()); + AssertDimension(x_dealii->size(), z_dealii->size()); + AssertDimension(x_dealii->size(), y_dealii->size()); - for (unsigned i = 0; i < x_dealii->n_blocks(); ++i) - { - auto x_ele = x_dealii->block(i).locally_owned_elements(); - for (const auto idx : x_ele) - { - z_dealii->block(i)[idx] = - x_dealii->block(i)[idx] / y_dealii->block(i)[idx]; - } - } - z_dealii->compress(VectorOperation::insert); -} + auto x_ele = x_dealii->locally_owned_elements(); + for (const auto idx : x_ele) + { + (*z_dealii)[idx] = (*x_dealii)[idx] / (*y_dealii)[idx]; + } + z_dealii->compress(VectorOperation::insert); + } -template ::value, int>> -void -SUNDIALS::internal::NVectorOperations::elementwise_inv(N_Vector x, N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *z_dealii = unwrap_nvector(z); + template ::value, int>> + void + elementwise_div(N_Vector x, N_Vector y, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *y_dealii = unwrap_nvector_const(y); + auto *z_dealii = unwrap_nvector(z); + + AssertDimension(x_dealii->size(), z_dealii->size()); + AssertDimension(x_dealii->size(), y_dealii->size()); + AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks()); + AssertDimension(x_dealii->n_blocks(), y_dealii->n_blocks()); + + for (unsigned i = 0; i < x_dealii->n_blocks(); ++i) + { + auto x_ele = x_dealii->block(i).locally_owned_elements(); + for (const auto idx : x_ele) + { + z_dealii->block(i)[idx] = + x_dealii->block(i)[idx] / y_dealii->block(i)[idx]; + } + } + z_dealii->compress(VectorOperation::insert); + } - AssertDimension(x_dealii->size(), z_dealii->size()); - auto x_ele = x_dealii->locally_owned_elements(); - for (const auto idx : x_ele) - { - (*z_dealii)[idx] = 1.0 / (*x_dealii)[idx]; - } - z_dealii->compress(VectorOperation::insert); -} + template ::value, int>> + void + elementwise_inv(N_Vector x, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *z_dealii = unwrap_nvector(z); + AssertDimension(x_dealii->size(), z_dealii->size()); -template ::value, int>> -void -SUNDIALS::internal::NVectorOperations::elementwise_inv(N_Vector x, N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *z_dealii = unwrap_nvector(z); + auto x_ele = x_dealii->locally_owned_elements(); + for (const auto idx : x_ele) + { + (*z_dealii)[idx] = 1.0 / (*x_dealii)[idx]; + } + z_dealii->compress(VectorOperation::insert); + } - AssertDimension(x_dealii->size(), z_dealii->size()); - AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks()); - for (unsigned i = 0; i < x_dealii->n_blocks(); ++i) - { - auto x_ele = x_dealii->block(i).locally_owned_elements(); - for (const auto idx : x_ele) - { - z_dealii->block(i)[idx] = 1.0 / x_dealii->block(i)[idx]; - } - } - z_dealii->compress(VectorOperation::insert); -} + template ::value, int>> + void + elementwise_inv(N_Vector x, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *z_dealii = unwrap_nvector(z); + AssertDimension(x_dealii->size(), z_dealii->size()); + AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks()); + for (unsigned i = 0; i < x_dealii->n_blocks(); ++i) + { + auto x_ele = x_dealii->block(i).locally_owned_elements(); + for (const auto idx : x_ele) + { + z_dealii->block(i)[idx] = 1.0 / x_dealii->block(i)[idx]; + } + } -template ::value, int>> -void -SUNDIALS::internal::NVectorOperations::elementwise_abs(N_Vector x, N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *z_dealii = unwrap_nvector(z); + z_dealii->compress(VectorOperation::insert); + } - AssertDimension(x_dealii->size(), z_dealii->size()); - auto x_ele = x_dealii->locally_owned_elements(); - for (const auto idx : x_ele) - { - (*z_dealii)[idx] = std::fabs((*x_dealii)[idx]); - } - z_dealii->compress(VectorOperation::insert); -} + template ::value, int>> + void + elementwise_abs(N_Vector x, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *z_dealii = unwrap_nvector(z); + AssertDimension(x_dealii->size(), z_dealii->size()); + auto x_ele = x_dealii->locally_owned_elements(); + for (const auto idx : x_ele) + { + (*z_dealii)[idx] = std::fabs((*x_dealii)[idx]); + } -template ::value, int>> -void -SUNDIALS::internal::NVectorOperations::elementwise_abs(N_Vector x, N_Vector z) -{ - auto *x_dealii = unwrap_nvector_const(x); - auto *z_dealii = unwrap_nvector(z); + z_dealii->compress(VectorOperation::insert); + } - AssertDimension(x_dealii->size(), z_dealii->size()); - AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks()); - for (unsigned i = 0; i < x_dealii->n_blocks(); ++i) - { - auto x_ele = x_dealii->block(i).locally_owned_elements(); - for (const auto idx : x_ele) - { - z_dealii->block(i)[idx] = std::fabs(x_dealii->block(i)[idx]); - } - } - z_dealii->compress(VectorOperation::insert); -} + template ::value, int>> + void + elementwise_abs(N_Vector x, N_Vector z) + { + auto *x_dealii = unwrap_nvector_const(x); + auto *z_dealii = unwrap_nvector(z); + + AssertDimension(x_dealii->size(), z_dealii->size()); + AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks()); + + for (unsigned i = 0; i < x_dealii->n_blocks(); ++i) + { + auto x_ele = x_dealii->block(i).locally_owned_elements(); + for (const auto idx : x_ele) + { + z_dealii->block(i)[idx] = std::fabs(x_dealii->block(i)[idx]); + } + } + + z_dealii->compress(VectorOperation::insert); + } + } // namespace NVectorOperations -template -N_Vector -SUNDIALS::internal::create_empty_nvector() -{ - N_Vector v = N_VNewEmpty(); - Assert(v != nullptr, ExcInternalError()); - - /* constructors, destructors, and utility operations */ - v->ops->nvgetvectorid = &NVectorOperations::get_vector_id; - v->ops->nvclone = &NVectorOperations::clone; - v->ops->nvcloneempty = &NVectorOperations::clone_empty; - v->ops->nvdestroy = &NVectorOperations::destroy; - // v->ops->nvspace = undef; + template + N_Vector + create_empty_nvector() + { + N_Vector v = N_VNewEmpty(); + Assert(v != nullptr, ExcInternalError()); + + /* constructors, destructors, and utility operations */ + v->ops->nvgetvectorid = &NVectorOperations::get_vector_id; + v->ops->nvclone = &NVectorOperations::clone; + v->ops->nvcloneempty = &NVectorOperations::clone_empty; + v->ops->nvdestroy = &NVectorOperations::destroy; + // v->ops->nvspace = undef; # if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0) - v->ops->nvgetcommunicator = - NVectorOperations::get_communicator_as_void_ptr; - v->ops->nvgetlength = &NVectorOperations::get_global_length; + v->ops->nvgetcommunicator = + &NVectorOperations::get_communicator_as_void_ptr; + v->ops->nvgetlength = &NVectorOperations::get_global_length; # endif - /* standard vector operations */ - v->ops->nvlinearsum = &NVectorOperations::linear_sum; - v->ops->nvconst = &NVectorOperations::set_constant; - v->ops->nvprod = &NVectorOperations::elementwise_product; - v->ops->nvdiv = &NVectorOperations::elementwise_div; - v->ops->nvscale = &NVectorOperations::scale; - v->ops->nvabs = &NVectorOperations::elementwise_abs; - v->ops->nvinv = &NVectorOperations::elementwise_inv; - v->ops->nvaddconst = &NVectorOperations::add_constant; - v->ops->nvdotprod = &NVectorOperations::dot_product; - v->ops->nvmaxnorm = &NVectorOperations::max_norm; - v->ops->nvwrmsnorm = &NVectorOperations::weighted_rms_norm; - v->ops->nvmin = &NVectorOperations::min_element; - v->ops->nvwl2norm = &NVectorOperations::weighted_l2_norm; - v->ops->nvl1norm = &NVectorOperations::l1_norm; - v->ops->nvwrmsnormmask = - &NVectorOperations::weighted_rms_norm_mask; - // v->ops->nvcompare = undef; - // v->ops->nvinvtest = undef; - // v->ops->nvconstrmask = undef; - // v->ops->nvminquotient = undef; - - /* fused and vector array operations are disabled (NULL) by default */ - - /* local reduction operations */ - // v->ops->nvdotprodlocal = undef; - // v->ops->nvmaxnormlocal = undef; - // v->ops->nvminlocal = undef; - // v->ops->nvl1normlocal = undef; - // v->ops->nvinvtestlocal = undef; - // v->ops->nvconstrmasklocal = undef; - // v->ops->nvminquotientlocal = undef; - // v->ops->nvwsqrsumlocal = undef; - // v->ops->nvwsqrsummasklocal = undef; - return v; -} + /* standard vector operations */ + v->ops->nvlinearsum = &NVectorOperations::linear_sum; + v->ops->nvconst = &NVectorOperations::set_constant; + v->ops->nvprod = &NVectorOperations::elementwise_product; + v->ops->nvdiv = &NVectorOperations::elementwise_div; + v->ops->nvscale = &NVectorOperations::scale; + v->ops->nvabs = &NVectorOperations::elementwise_abs; + v->ops->nvinv = &NVectorOperations::elementwise_inv; + v->ops->nvaddconst = &NVectorOperations::add_constant; + v->ops->nvdotprod = &NVectorOperations::dot_product; + v->ops->nvmaxnorm = &NVectorOperations::max_norm; + v->ops->nvwrmsnorm = &NVectorOperations::weighted_rms_norm; + v->ops->nvmin = &NVectorOperations::min_element; + v->ops->nvwl2norm = &NVectorOperations::weighted_l2_norm; + v->ops->nvl1norm = &NVectorOperations::l1_norm; + v->ops->nvwrmsnormmask = + &NVectorOperations::weighted_rms_norm_mask; + // v->ops->nvcompare = undef; + // v->ops->nvinvtest = undef; + // v->ops->nvconstrmask = undef; + // v->ops->nvminquotient = undef; + + /* fused and vector array operations are disabled (NULL) by default */ + + /* local reduction operations */ + // v->ops->nvdotprodlocal = undef; + // v->ops->nvmaxnormlocal = undef; + // v->ops->nvminlocal = undef; + // v->ops->nvl1normlocal = undef; + // v->ops->nvinvtestlocal = undef; + // v->ops->nvconstrmasklocal = undef; + // v->ops->nvminquotientlocal = undef; + // v->ops->nvwsqrsumlocal = undef; + // v->ops->nvwsqrsummasklocal = undef; + return v; + } + + } // namespace internal +} // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE -- 2.39.5