#13: Add elementwise_reciprocal for tpetra and tpetra block

.gitignore

@@ -1,2 +1,5 @@
 .vscode
 generated_docs
+build
+install
+

docs/source/components/ops.rst

@@ -88,6 +88,7 @@ Use the links below to view details about each operation.
     ops/deep_copy
     ops/dot
     ops/elementwise_multiply
+    ops/elementwise_reciprocal
     ops/extent
     ops/fill
     ops/min_max

docs/source/components/ops/elementwise_reciprocal.rst

@@ -0,0 +1,38 @@
+.. include:: ../../mydefs.rst
+
+``elementwise_reciprocal``
+========================
+
+Header: ``<pressio/ops.hpp>``
+
+API
+---
+
+.. code-block:: cpp
+
+  namespace pressio { namespace ops{
+
+  template <class T1, class T2>
+  void elementwise_multiply(const T1 & z, T2 & y);
+
+  }} // end namespace pressio::ops
+
+Description
+-----------
+
+* Computes the elementwise reciprocal of ``z`` and stores it in ``y``
+
+* ``T1`` and ``T2`` must be one-dimensional containers with compatible scalar types. They may be:
+
+  * Eigen vectors or Pressio expressions built from Eigen rank-1 containers
+
+  * Kokkos rank-1 views or Pressio expressions built from Kokkos rank-1 containers
+
+  * Tpetra vectors or column expressions
+
+  * Tpetra block vectors or column expressions
+
+Notes
+-----
+
+* See the :doc:`ops homepage <../ops>` for a table of booleans to use when checking that ``T1`` and ``T2`` have the correct types.

include/pressio/ops.hpp

@@ -139,6 +139,7 @@ template<class ...> struct matching_extents;
 #include "ops/tpetra/ops_pow.hpp"
 #include "ops/tpetra/ops_rank1_update.hpp"
 #include "ops/tpetra/ops_elementwise_multiply.hpp"
+#include "ops/tpetra/ops_elementwise_reciprocal.hpp"
 #include "ops/tpetra/ops_multi_vector_update.hpp"
 #include "ops/tpetra/ops_level2.hpp"
 #include "ops/tpetra/ops_level3.hpp"
@@ -157,6 +158,7 @@ template<class ...> struct matching_extents;
 #include "ops/tpetra_block/ops_pow.hpp"
 #include "ops/tpetra_block/ops_rank1_update.hpp"
 #include "ops/tpetra_block/ops_elementwise_multiply.hpp"
+#include "ops/tpetra_block/ops_elementwise_reciprocal.hpp"
 #include "ops/tpetra_block/ops_multi_vector_update.hpp"
 #include "ops/tpetra_block/ops_level2.hpp"
 #include "ops/tpetra_block/ops_level3.hpp"

include/pressio/ops/tpetra/ops_elementwise_reciprocal.hpp

@@ -0,0 +1,78 @@
+/*
+//@HEADER
+// ************************************************************************
+//
+// ops_elementwise_reciprocal.hpp
+//                     		  Pressio
+//                             Copyright 2019
+//    National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+//
+// Under the terms of Contract DE-NA0003525 with NTESS, the
+// U.S. Government retains certain rights in this software.
+//
+// Pressio is licensed under BSD-3-Clause terms of use:
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+// contributors may be used to endorse or promote products derived
+// from this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Questions? Contact Francesco Rizzi (fnrizzi@sandia.gov)
+//
+// ************************************************************************
+//@HEADER
+*/
+
+#ifndef PRESSIOOPS_OPS_TPETRA_OPS_ELEMENTWISE_RECIPROCAL_HPP_
+#define PRESSIOOPS_OPS_TPETRA_OPS_ELEMENTWISE_RECIPROCAL_HPP_
+
+namespace pressio{ namespace ops{
+
+//----------------------------------------------------------------------
+// computing elementwise:  y = 1/z
+//----------------------------------------------------------------------
+template <typename T1, typename T2>
+std::enable_if_t<
+  // TPL/container specific
+     (   ::pressio::is_vector_tpetra<T1>::value
+      || ::pressio::is_expression_column_acting_on_tpetra<T1>::value)
+  && (   ::pressio::is_vector_tpetra<T2>::value
+      || ::pressio::is_expression_column_acting_on_tpetra<T2>::value)
+  // scalar compatibility
+  && ::pressio::all_have_traits_and_same_scalar<T1, T2>::value
+  && (std::is_floating_point<typename ::pressio::Traits<T1>::scalar_type>::value
+   || std::is_integral<typename ::pressio::Traits<T1>::scalar_type>::value)
+  >
+elementwise_reciprocal(const T1 & zin, T2 & yin)
+{
+  auto z = impl::get_native(zin);
+  auto y = impl::get_native(yin);
+  assert(::pressio::ops::extent(z, 0)==::pressio::ops::extent(y, 0));
+  y.reciprocal(z);
+}
+
+}}//end namespace pressio::ops
+#endif  // PRESSIOOPS_OPS_TPETRA_OPS_ELEMENTWISE_RECIPROCAL_HPP_

include/pressio/ops/tpetra_block/ops_elementwise_reciprocal.hpp

@@ -0,0 +1,77 @@
+/*
+//@HEADER
+// ************************************************************************
+//
+// ops_elementwise_reciprocal.hpp
+//                     		  Pressio
+//                             Copyright 2019
+//    National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+//
+// Under the terms of Contract DE-NA0003525 with NTESS, the
+// U.S. Government retains certain rights in this software.
+//
+// Pressio is licensed under BSD-3-Clause terms of use:
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+// contributors may be used to endorse or promote products derived
+// from this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Questions? Contact Francesco Rizzi (fnrizzi@sandia.gov)
+//
+// ************************************************************************
+//@HEADER
+*/
+
+#ifndef PRESSIOOPS_OPS_TPETRA_BLOCK_OPS_ELEMENTWISE_RECIPROCAL_HPP_
+#define PRESSIOOPS_OPS_TPETRA_BLOCK_OPS_ELEMENTWISE_RECIPROCAL_HPP_
+
+namespace pressio{ namespace ops{
+
+//----------------------------------------------------------------------
+// computing elementwise:  y = 1/z
+//----------------------------------------------------------------------
+template <typename T1, typename T2>
+std::enable_if_t<
+     (   ::pressio::is_vector_tpetra_block<T1>::value
+      || ::pressio::is_expression_column_acting_on_tpetra_block<T1>::value)
+  && (   ::pressio::is_vector_tpetra_block<T2>::value
+      || ::pressio::is_expression_column_acting_on_tpetra_block<T2>::value)
+  // scalar compatibility
+  && ::pressio::all_have_traits_and_same_scalar<T1, T2>::value
+  && (std::is_floating_point<typename ::pressio::Traits<T1>::scalar_type>::value
+   || std::is_integral<typename ::pressio::Traits<T1>::scalar_type>::value)
+  >
+elementwise_reciprocal(const T1 & z, T2 & y)
+{
+  assert(extent(z,0)==extent(y,0));
+  auto z_tp = impl::get_underlying_tpetra_object(z);
+  auto y_tp = impl::get_underlying_tpetra_object(y);
+  y_tp.reciprocal(z_tp);
+}
+
+}}//end namespace pressio::ops
+#endif  // PRESSIOOPS_OPS_TPETRA_BLOCK_OPS_ELEMENTWISE_RECIPROCAL_HPP_

tests/ops/ops_tpetra_block_column_expr.cc

@@ -351,3 +351,21 @@ TEST_F(ops_tpetra_block, column_expr_elementwiseMultiply)
     EXPECT_DOUBLE_EQ(y_h(i,0), 12);
   }
 }
+
+TEST_F(ops_tpetra_block, column_expr_elementwiseReciprocal)
+{
+  // computing elementwise:  y = 1/z
+
+  auto z = pressio::column(*myMv_, 1);
+  pressio::ops::fill(z, 2.);
+  auto z_h = z.native().getLocalViewHost(Tpetra::Access::ReadOnlyStruct());
+
+  auto y = pressio::column(*myMv_, 2);
+
+  pressio::ops::elementwise_reciprocal(z, y);
+
+  auto y_h = y.native().getLocalViewHost(Tpetra::Access::ReadOnlyStruct());
+  for (int i=0; i<localSize_*blockSize_; ++i){
+    EXPECT_DOUBLE_EQ(y_h(i,0), 1 / z_h(i,0));
+  }
+}

tests/ops/ops_tpetra_block_vector.cc

@@ -412,3 +412,20 @@ TEST_F(ops_tpetra_block, vector_elementwiseMultiply)
   }
 
 }
+
+TEST_F(ops_tpetra_block, vector_elementwiseReciprocal)
+{
+  // computing elementwise:  y = 1/z
+
+  auto z = pressio::ops::clone(*myVector_);
+  pressio::ops::fill(z, 2.);
+  auto z_h = z.getVectorView().getLocalViewHost(Tpetra::Access::ReadOnlyStruct());
+
+  auto y = pressio::ops::clone(*myVector_);
+
+  pressio::ops::elementwise_reciprocal(z, y);
+  auto y_h = y.getVectorView().getLocalViewHost(Tpetra::Access::ReadWriteStruct());
+  for (int i=0; i<localSize_*blockSize_; ++i){
+    EXPECT_DOUBLE_EQ(y_h(i,0), 1 / z_h(i,0));
+  }
+}

tests/ops/ops_tpetra_column_expr.cc

@@ -504,3 +504,30 @@ TEST_F(ops_tpetra, column_expr_elementwiseMultiply)
     EXPECT_DOUBLE_EQ(ch(i,0), gold + offset);
   }
 }
+
+TEST_F(ops_tpetra, column_expr_elementwiseReciprocal)
+{
+  // 1. create z and fill with data
+  mvec_t z(contigMap_, 3);
+  auto zh = z.getLocalViewHost(Tpetra::Access::ReadWriteStruct());
+  for (int i=0; i<localSize_; ++i) {
+    for (int j=0; j<3; ++j){
+      zh(i,j) = static_cast<ST>(i + j + 1);
+    }
+  }
+
+  // 2. create y (will store reciprocal of z)
+  mvec_t y(contigMap_, 3);
+
+  // 3. compute ewise reciprocal
+  auto z_expr = pressio::column(z, 0);
+  auto y_expr = pressio::column(y, 0);
+  pressio::ops::elementwise_reciprocal(z_expr, y_expr);
+
+  // 4. check correctness
+  auto yh = y.getLocalViewHost(Tpetra::Access::ReadOnlyStruct());
+  for (int i=0; i<localSize_; ++i){
+    const auto gold = 1 / zh(i, 0);
+    EXPECT_DOUBLE_EQ(yh(i,0), gold);
+  }
+}

tests/ops/ops_tpetra_vector.cc

@@ -473,3 +473,26 @@ TEST_F(ops_tpetra, vector_elementwiseMultiply)
     EXPECT_DOUBLE_EQ(ch(i,0), gold + offset);
   }
 }
+
+TEST_F(ops_tpetra, vector_elementwiseReciprocal)
+{
+  // 1. Create and fill vector z with test data
+  vec_t z(contigMap_);
+  auto zh = z.getLocalViewHost(Tpetra::Access::ReadWriteStruct());
+  for (int i = 0; i < localSize_; ++i) {
+    zh(i, 0) = static_cast<ST>(i + 1);
+  }
+
+  // 2. Create vector y, which will store the elementwise reciprocal of z
+  vec_t y(contigMap_);
+
+  // 3. Compute elementwise reciprocal: y = 1 / z
+  pressio::ops::elementwise_reciprocal(z, y);
+
+  // 4. Check correctness of the reciprocal
+  auto yh = y.getLocalViewHost(Tpetra::Access::ReadOnlyStruct());
+  for (int i = 0; i < localSize_; ++i) {
+    const auto expected_value = 1.0 / zh(i, 0);
+    EXPECT_DOUBLE_EQ(yh(i, 0), expected_value);
+  }
+}