diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml
index 7bdfc69be..b1a4e778a 100644
--- a/.github/workflows/windows.yml
+++ b/.github/workflows/windows.yml
@@ -24,7 +24,7 @@ jobs:
 #           - {os: windows-2016, toolset: msvc, version: 14.16, cxxstd: 11} 
 #           - {os: windows-2019, toolset: msvc, version: 14.28, cxxstd: 11}
 #           - {os: windows-2019, toolset: msvc, version: 14.28, cxxstd: 17}
-           - {os: windows-2019, toolset: msvc, version: 14.28, cxxstd: latest}
+           - {os: windows-2019, toolset: msvc, version: 14.29, cxxstd: latest}
              
      steps:
       - uses: actions/checkout@v2
diff --git a/include/boost/numeric/ublas/tensor/extents/extents_functions.hpp b/include/boost/numeric/ublas/tensor/extents/extents_functions.hpp
index 149e12bde..cfa247b4d 100644
--- a/include/boost/numeric/ublas/tensor/extents/extents_functions.hpp
+++ b/include/boost/numeric/ublas/tensor/extents/extents_functions.hpp
@@ -77,6 +77,37 @@ template <class D>
 //         std::all_of(cbegin(e)+2,cend(e)    , [](auto a){return a==1UL;});
 }
 
+/** @brief Returns true if extents equals (m,[1,1,...,1]) with m>=1 */
+template <class D>
+[[nodiscard]] inline constexpr bool is_row_vector(extents_base<D> const& e)
+{
+  if (empty(e) || size(e) == 1  ) {return false;}
+
+  if(cbegin(e)[0] ==  1ul &&
+      cbegin(e)[1] >  1ul &&
+      std::all_of(cbegin(e)+2ul,cend  (e) , [](auto a){return a==1ul;})){
+    return true;
+  }
+
+  return false;
+}
+
+
+/** @brief Returns true if extents equals (m,[1,1,...,1]) with m>=1 */
+template <class D>
+[[nodiscard]] inline constexpr bool is_col_vector(extents_base<D> const& e)
+{
+  if (empty(e) || size(e) == 1  ) {return false;}
+
+  if(cbegin(e)[0] >  1ul &&
+      cbegin(e)[1] == 1ul &&
+      std::all_of(cbegin(e)+2ul,cend  (e) , [](auto a){return a==1ul;})){
+    return true;
+  }
+
+  return false;
+}
+
 /** @brief Returns true if (m,[n,1,...,1]) with m>=1 or n>=1 */
 template <class D>
 [[nodiscard]] inline constexpr bool is_matrix(extents_base<D> const& e)
diff --git a/include/boost/numeric/ublas/tensor/function/tensor_times_vector.hpp b/include/boost/numeric/ublas/tensor/function/tensor_times_vector.hpp
index 5702f97fb..611c5a82e 100644
--- a/include/boost/numeric/ublas/tensor/function/tensor_times_vector.hpp
+++ b/include/boost/numeric/ublas/tensor/function/tensor_times_vector.hpp
@@ -15,6 +15,7 @@
 #include <stdexcept>
 #include <type_traits>
 
+#include "../multiplication.hpp"
 #include "../extents.hpp"
 #include "../type_traits.hpp"
 #include "../tags.hpp"
@@ -49,6 +50,190 @@ using enable_ttv_if_extent_has_dynamic_rank = std::enable_if_t<is_dynamic_rank_v
 } // namespace detail
 
 
+namespace detail {
+template <class TC, class TA, class  V, class EC>
+inline auto scalar_scalar_prod(TA const &a, V const &b, EC const& nc_base)
+{
+  assert(ublas::is_scalar(a.extents()));
+  using tensor = TC;
+  using value = typename tensor::value_type;
+  using shape = typename tensor::extents_type;
+  return tensor(shape(nc_base),value(a[0]*b(0)));
+}
+
+template <class TC, class TA, class  V, class EC>
+inline auto vector_vector_prod(TA const &a, V const &b, EC& nc_base, std::size_t m)
+{
+  auto const& na = a.extents();
+
+  assert( ublas::is_vector(na));
+  assert(!ublas::is_scalar(na));
+  assert( ublas::size(na) > 1u);
+  assert(m > 0);
+
+  using tensor = TC;
+  using value = typename tensor::value_type;
+  using shape = typename tensor::extents_type;
+
+  auto const n1 = na[0];
+  auto const n2 = na[1];
+  auto const s = b.size();
+
+  // general
+  // [n1 n2 1 ... 1] xj [s 1] for any 1 <= j <= p with n1==1 or n2==1
+
+
+  // [n1 1 1 ... 1] x1 [n1 1] -> [1 1 1 ... 1]
+  // [1 n2 1 ... 1] x2 [n2 1] -> [1 1 1 ... 1]
+
+
+  assert(n1>1 || n2>1);
+
+  if( (n1>1u && m==1u)  || (n2>1u && m==2u) ){
+    if(m==1u) assert(n2==1u && n1==s);
+    if(m==2u) assert(n1==1u && n2==s);
+    auto cc = std::inner_product( a.begin(), a.end(), b.begin(), value(0) );
+    return tensor(shape(nc_base),value(cc));
+  }
+
+  // [n1 1 1 ... 1] xj [1 1] -> [n1 1 1 ... 1] with j != 1
+  // [1 n2 1 ... 1] xj [1 1] -> [1 n2 1 ... 1] with j != 2
+
+//if( (n1>1u && m!=1u) && (n2>0u && m!=2u) ){
+
+  if(n1>1u) assert(m!=1u);
+  if(n2>1u) assert(m!=2u);
+  assert(s==1u);
+
+  if(n1>1u) assert(n2==1u);
+  if(n2>1u) assert(n1==1u);
+
+  if(n1>1u) nc_base[0] = n1;
+  if(n2>1u) nc_base[1] = n2;
+
+  auto bb = b(0);
+  auto c = tensor(shape(nc_base));
+  std::transform(a.begin(),a.end(),c.begin(),[bb](auto aa){ return aa*bb; });
+  return c;
+//}
+
+
+}
+
+
+/** Computes a matrix-vector product.
+ *
+ *
+ *  @note assume stride 1 for specific dimensions and therefore requires refactoring for subtensor
+ *
+*/
+template <class TC, class TA, class  V, class EC>
+inline auto matrix_vector_prod(TA const &a, V const &b, EC& nc_base, std::size_t m)
+{
+  auto const& na = a.extents();
+
+  assert( ublas::is_matrix(na));
+  assert(!ublas::is_vector(na));
+  assert(!ublas::is_scalar(na));
+  assert( ublas::size(na) > 1u);
+  assert(m > 0);
+
+  using tensor = TC;
+  using shape  = typename tensor::extents_type;
+  using size_t = typename shape::value_type;
+
+  auto const n1 = na[0];
+  auto const n2 = na[1];
+  auto const s = b.size();
+
+  // general
+  // [n1 n2 1 ... 1] xj [s 1] for any 1 <= j <= p with either n1>1 and n2>1
+
+
+  // if [n1 n2 1 ... 1] xj [1 1] -> [n1 n2 1 ... 1] for j > 2
+  if(m > 2){
+    nc_base[0] = n1;
+    nc_base[1] = n2;
+    assert(s == 1);
+    auto c  = tensor(shape(nc_base));
+    auto const bb = b(0);
+    std::transform(a.begin(),a.end(), c.begin(), [bb](auto aa){return aa*bb;});
+    return c;
+  }
+
+
+  // [n1 n2 1 ... 1] x1 [n1 1] -> [n2 1 ... 1] -> vector-times-matrix
+  // [n1 n2 1 ... 1] x2 [n2 1] -> [n1 1 ... 1] -> matrix-times-vector
+
+  nc_base[0] = m==1 ? n2 : n1;
+
+  auto c  = tensor(shape(nc_base));
+  auto const& wa = a.strides();
+  auto const* bdata = &(b(0));
+
+  detail::recursive::mtv(m-1,n1,n2, c.data(), size_t(1), a.data(), wa[0], wa[1], bdata, size_t(1));
+
+  return c;
+}
+
+
+
+template <class TC, class TA, class  V, class EC>
+inline auto tensor_vector_prod(TA const &a, V const &b, EC& nc_base, std::size_t m)
+{
+  auto const& na = a.extents();
+
+  assert( ublas::is_tensor(na));
+  assert( ublas::size(na) > 1u);
+  assert(m > 0);
+
+  using tensor = TC;
+  using shape  = typename tensor::extents_type;
+  using layout = typename tensor::layout_type;
+
+  auto const pa = a.rank();
+  auto const nm = na[m-1];
+  auto const s = b.size();
+
+  auto nb = extents<2>{std::size_t(b.size()),std::size_t(1ul)};
+  auto wb = ublas::to_strides(nb,layout{} );
+
+  //TODO: Include an outer product when legacy vector becomes a new vector.
+
+  for (auto i = 0ul, j = 0ul; i < pa; ++i)
+    if (i != m - 1)
+      nc_base[j++] = na.at(i);
+
+  auto c  = tensor(shape(nc_base));
+
+  // [n1 n2 ... nm ... np] xm [1 1] -> [n1 n2 ... nm-1 nm+1 ... np]
+
+  if(s == 0){
+    assert(nm == 1);
+    auto const bb = b(0);
+    std::transform(a.begin(),a.end(), c.begin(), [bb](auto aa){return aa*bb;});
+    return c;
+  }
+
+
+  // if [n1 n2 n3 ... np] xm [nm 1] -> [n1 n2 ... nm-1 nm+1 ... np]
+
+  auto const& nc = c.extents();
+  auto const& wc = c.strides();
+  auto const& wa = a.strides();
+  auto const* bp = &(b(0));
+
+  ttv(m, pa,
+      c.data(), nc.data(), wc.data(),
+      a.data(), na.data(), wa.data(),
+      bp,       nb.data(), wb.data());
+
+  return c;
+}
+
+}//namespace detail
+
+
 /** @brief Computes the m-mode tensor-times-vector product
      *
      * Implements C[i1,...,im-1,im+1,...,ip] = A[i1,i2,...,ip] * b[im]
@@ -63,45 +248,49 @@ using enable_ttv_if_extent_has_dynamic_rank = std::enable_if_t<is_dynamic_rank_v
     */
 template <class TE, class  A, class  T = typename tensor_core< TE >::value,
           detail::enable_ttv_if_extent_has_dynamic_rank<TE> = true >
-inline decltype(auto) prod( tensor_core< TE > const &a, vector<T, A> const &b, const std::size_t m)
+inline auto prod( tensor_core< TE > const &a, vector<T, A> const &b, const std::size_t m)
 {
 
   using tensor            = tensor_core< TE >;
   using shape             = typename tensor::extents_type;
-  using value             = typename tensor::value_type;
-  using layout            = typename tensor::layout_type;
   using resize_tag        = typename tensor::resizable_tag;
 
-  auto const p = a.rank();
+  auto const pa = a.rank();
 
   static_assert(std::is_same_v<resize_tag,storage_resizable_container_tag>);
   static_assert(is_dynamic_v<shape>);
 
   if (m == 0ul)  throw std::length_error("error in boost::numeric::ublas::prod(ttv): contraction mode must be greater than zero.");
-  if (p < m)     throw std::length_error("error in boost::numeric::ublas::prod(ttv): rank of tensor must be greater than or equal to the contraction mode.");
+  if (pa < m)    throw std::length_error("error in boost::numeric::ublas::prod(ttv): rank of tensor must be greater than or equal to the contraction mode.");
   if (a.empty()) throw std::length_error("error in boost::numeric::ublas::prod(ttv): first argument tensor should not be empty.");
   if (b.empty()) throw std::length_error("error in boost::numeric::ublas::prod(ttv): second argument vector should not be empty.");
 
   auto const& na = a.extents();
-  auto nb = extents<2>{std::size_t(b.size()),std::size_t(1ul)};
-  auto wb = ublas::to_strides(nb,layout{} );
+
+  if(b.size() != na[m-1]) throw std::length_error("error in boost::numeric::ublas::prod(ttv): dimension mismatch of tensor and vector.");
 
   auto const sz = std::max( std::size_t(ublas::size(na)-1u), std::size_t(2) );
   auto nc_base = typename shape::base_type(sz,1);
 
-  for (auto i = 0ul, j = 0ul; i < p; ++i)
-    if (i != m - 1)
-      nc_base[j++] = na.at(i);
+  // output scalar tensor
+  if(ublas::is_scalar(na)){
+    return detail::scalar_scalar_prod<tensor>(a,b,nc_base);
+  }
+
+  // output scalar tensor or vector tensor
+  if (ublas::is_vector(na)){
+    return detail::vector_vector_prod<tensor>(a,b,nc_base,m);
+  }
+
+  // output scalar tensor or vector tensor
+  if (ublas::is_matrix(na)){
+    return detail::matrix_vector_prod<tensor>(a,b,nc_base,m);
+  }
+
+  assert(ublas::is_tensor(na));
+  return detail::tensor_vector_prod<tensor>(a,b,nc_base,m);
 
-  auto nc = shape(nc_base);
-  auto c = tensor( nc, value{} );
 
-  auto const* bb = &(b(0));
-  ttv(m, p,
-      c.data(), c.extents().data(), c.strides().data(),
-      a.data(), a.extents().data(), a.strides().data(),
-      bb,       nb.data(),          wb.data());
-  return c;
 }
 
 
@@ -143,7 +332,6 @@ inline auto prod( tensor_core< TE > const &a, vector<T, A> const &b, const std::
   constexpr auto p  = std::tuple_size_v<shape>;
   constexpr auto sz = std::max(std::size_t(std::tuple_size_v<shape>-1U),std::size_t(2));
 
-  using shape_b   = ublas::extents<2>;
   using shape_c   = ublas::extents<sz>;
   using tensor_c  = tensor_core<tensor_engine<shape_c,layout,container>>;
 
@@ -158,21 +346,25 @@ inline auto prod( tensor_core< TE > const &a, vector<T, A> const &b, const std::
 
   auto nc_base = typename shape_c::base_type{};
   std::fill(nc_base.begin(), nc_base.end(),std::size_t(1));
-  for (auto i = 0ul, j = 0ul; i < p; ++i)
-    if (i != m - 1)
-      nc_base[j++] = na.at(i);
 
-  auto nc = shape_c(std::move(nc_base));
-  auto nb = shape_b{b.size(),1UL};
-  auto wb = ublas::to_strides(nb,layout{});
-  auto c  = tensor_c( std::move(nc) );
-  auto const* bb = &(b(0));
 
-  ttv(m, p,
-      c.data(), c.extents().data(), c.strides().data(),
-      a.data(), a.extents().data(), a.strides().data(),
-      bb,       nb.data(),          wb.data() );
-  return c;
+  // output scalar tensor
+  if(ublas::is_scalar(na)){
+    return detail::scalar_scalar_prod<tensor_c>(a,b,nc_base);
+  }
+
+  // output scalar tensor or vector tensor
+  if (ublas::is_vector(na)){
+    return detail::vector_vector_prod<tensor_c>(a,b,nc_base,m);
+  }
+
+  // output scalar tensor or vector tensor
+  if (ublas::is_matrix(na)){
+    return detail::matrix_vector_prod<tensor_c>(a,b,nc_base,m);
+  }
+
+  assert(ublas::is_tensor(na));
+  return detail::tensor_vector_prod<tensor_c>(a,b,nc_base,m);
 }
 
 
@@ -201,7 +393,7 @@ inline auto prod( tensor_core< TE > const &a, vector<T, A> const &b)
   using shape         = typename tensor::extents;
   using layout        = typename tensor::layout;
   using shape_b       = extents<2>;
-  using shape_c       = remove_element_t<m,shape>;
+  using shape_c       = remove_element_t<m,shape>; // this is wrong
   using container_c   = rebind_storage_size_t<shape_c,container>;
   using tensor_c      = tensor_core<tensor_engine<shape_c,layout,container_c>>;
 
diff --git a/include/boost/numeric/ublas/tensor/multiplication.hpp b/include/boost/numeric/ublas/tensor/multiplication.hpp
index e2d94f7be..ea7901814 100644
--- a/include/boost/numeric/ublas/tensor/multiplication.hpp
+++ b/include/boost/numeric/ublas/tensor/multiplication.hpp
@@ -337,33 +337,44 @@ void ttv0(SizeType const r,
 
 /** @brief Computes the matrix-times-vector product
  *
- * Implements C[i1] = sum(A[i1,i2] * b[i2]) or C[i2] = sum(A[i1,i2] * b[i1])
- *
- * @note is used in function ttv
- *
- * @param[in]  m  zero-based contraction mode with m=0 or m=1
- * @param[out] c  pointer to the output tensor C
- * @param[in]  nc pointer to the extents of tensor C
- * @param[in]  wc pointer to the strides of tensor C
- * @param[in]  a  pointer to the first input tensor A
- * @param[in]  na pointer to the extents of input tensor A
- * @param[in]  wa pointer to the strides of input tensor A
- * @param[in]  b  pointer to the second input tensor B
+ * Implements C[i1] = sum(A[i1,i2] * B[i2]) if k = 1 or C[i2] = sum(A[i1,i2] * B[i1]) if k = 0
+ *
+ * [m,n] = size(A(..,:,..,:,..))
+ * [m]   = size(C(..,:,..))
+ * [n]   = size(B(..,:,..))
+ *
+ *
+ * @param[in]  k     if k = 0
+ * @param[in]  m     number of rows of A
+ * @param[in]  n     number of columns of A
+ * @param[out] c     pointer to C
+ * @param[in]  wc    m-th (k=1) or n-th (k=0) stride for C
+ * @param[in]  a     pointer to A
+ * @param[in]  wa_m  m-th stride for A
+ * @param[in]  wa_n  n-th stride for A
+ * @param[in]  b     pointer to B
+ * @param[in]  wb    n-th (k=1) or m-th (k=0) stride for B
 */
 template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
-void mtv(SizeType const m,
-         PointerOut c, SizeType const*const /*unsed*/, SizeType const*const wc,
-         PointerIn1 a, SizeType const*const na       , SizeType const*const wa,
-         PointerIn2 b)
+void mtv(SizeType const k,
+         SizeType const m,
+         SizeType const n,
+         PointerOut c, SizeType const wc,
+         PointerIn1 a, SizeType const wa_m, SizeType const wa_n,
+         PointerIn2 b, SizeType const wb)
 {
-  // decides whether matrix multiplied with vector or vector multiplied with matrix
-  const auto o = (m == 0) ? 1 : 0;
 
-  for(auto io = 0u; io < na[o]; c += wc[o], a += wa[o], ++io) {
+  auto const wa_x = k==0 ? wa_n : wa_m;
+  auto const wa_y = k==0 ? wa_m : wa_n;
+
+  auto const x = k==0 ? n : m;
+  auto const y = k==0 ? m : n;
+
+  for(auto ix = 0u; ix < x; c += wc, a += wa_x, ++ix) {
     auto c1 = c;
     auto a1 = a;
     auto b1 = b;
-    for(auto im = 0u; im < na[m]; a1 += wa[m], ++b1, ++im) {
+    for(auto iy = 0u; iy < y; a1 += wa_y, b1 += wb, ++iy) {
       *c1 += *a1 * *b1;
     }
   }
@@ -603,10 +614,12 @@ namespace boost::numeric::ublas {
  *   C[i1,i2,...,im-1,im+1,...,ip] = sum(A[i1,i2,...,im,...,ip] * b[im]) for m>1 and
  *   C[i2,...,ip]                  = sum(A[i1,...,ip]           * b[i1]) for m=1
  *
- * @note calls detail::ttv, detail::ttv0 or detail::mtv
+ * @note calls detail::ttv, detail::ttv0 for p == 1 or p == 2 use ublas::inner or ublas::mtv or ublas::vtm
+ *
+ *
  *
  * @param[in]  m  contraction mode with 0 < m <= p
- * @param[in]  p  number of dimensions (rank) of the first input tensor with p > 0
+ * @param[in]  p  number of dimensions (rank) of the first input tensor with p > 2
  * @param[out] c  pointer to the output tensor with rank p-1
  * @param[in]  nc pointer to the extents of tensor c
  * @param[in]  wc pointer to the strides of tensor c
@@ -621,50 +634,25 @@ template <class PointerOut, class PointerIn1, class PointerIn2, class SizeType>
 void ttv(SizeType const m, SizeType const p,
          PointerOut c,       SizeType const*const nc, SizeType const*const wc,
          const PointerIn1 a, SizeType const*const na, SizeType const*const wa,
-         const PointerIn2 b, SizeType const*const nb, SizeType const*const wb)
+         const PointerIn2 b, SizeType const*const nb, SizeType const*const /*unused*/)
 {
-  static_assert( std::is_pointer<PointerOut>::value && std::is_pointer<PointerIn1>::value & std::is_pointer<PointerIn2>::value,
+  static_assert( std::is_pointer<PointerOut>::value && std::is_pointer<PointerIn1>::value && std::is_pointer<PointerIn2>::value,
                 "Static error in boost::numeric::ublas::ttv: Argument types for pointers are not pointer types.");
 
-  if( m == 0){
-    throw std::length_error("Error in boost::numeric::ublas::ttv: Contraction mode must be greater than zero.");
-  }
+  assert(m != 0);
+  assert(p >= m);
+  assert(p >= 2);
 
-  if( p < m ){
-    throw std::length_error("Error in boost::numeric::ublas::ttv: Rank must be greater equal the modus.");
-  }
-  if( p == 0){
-    throw std::length_error("Error in boost::numeric::ublas::ttv: Rank must be greater than zero.");
-  }
-  if(c == nullptr || a == nullptr || b == nullptr){
-    throw std::length_error("Error in boost::numeric::ublas::ttv: Pointers shall not be null pointers.");
-  }
   for(auto i = 0u; i < m-1; ++i){
-    if(na[i] != nc[i]){
-      throw std::length_error("Error in boost::numeric::ublas::ttv: Extents (except of dimension mode) of A and C must be equal.");
-    }
-  }
-
-  const auto max = std::max(nb[0], nb[1]);
-  if(  na[m-1] != max){
-    throw std::length_error("Error in boost::numeric::ublas::ttv: Extent of dimension mode of A and b must be equal.");
+    assert(na[i] == nc[i]);
   }
 
+  assert(na[m-1] == std::max(nb[0], nb[1]));
 
-  if((m != 1) && (p > 2)){
-    detail::recursive::ttv(m-1, p-1, p-2, c, nc, wc,    a, na, wa,   b);
-  }
-  else if ((m == 1) && (p > 2)){
+  if(m == 1)
     detail::recursive::ttv0(p-1, c, nc, wc,  a, na, wa,   b);
-  }
-  else if( p == 2 ){
-    detail::recursive::mtv(m-1, c, nc, wc,  a, na, wa,   b);
-  }
-  else /*if( p == 1 )*/{
-    auto v = std::remove_pointer_t<std::remove_cv_t<PointerOut>>{};
-    *c = detail::recursive::inner(SizeType(0), na, a, wa, b, wb, v);
-  }
-
+  else
+    detail::recursive::ttv (m-1, p-1, p-2, c, nc, wc,    a, na, wa,   b);
 }
 
 
diff --git a/include/boost/numeric/ublas/tensor/subtensor.hpp b/include/boost/numeric/ublas/tensor/subtensor.hpp
index e8a150d16..a119300b7 100644
--- a/include/boost/numeric/ublas/tensor/subtensor.hpp
+++ b/include/boost/numeric/ublas/tensor/subtensor.hpp
@@ -129,10 +129,10 @@ class subtensor <tag::sliced, tensor_dynamic<T,F>>
 	subtensor(tensor_type& t, span_types&& ... spans)
 			: super_type     ()
       , spans_         (detail::generate_span_vector<span_type>(t.extents(),std::forward<span_types>(spans)...))
-      , extents_       (detail::compute_extents(spans_))
+      , extents_       (detail::to_extents(spans_))
     , strides_         (ublas::to_strides(extents_,layout_type{}))
-      , span_strides_  (detail::compute_span_strides(t.strides(),spans_))
-      , data_          {t.data() + detail::compute_offset(t.strides(), spans_)}
+      , span_strides_  (detail::to_span_strides(t.strides(),spans_))
+      , data_          {t.data() + detail::to_offset(t.strides(), spans_)}
 	{
 //		if( m == nullptr)
 //			throw std::length_error("Error in tensor_view<T>::tensor_view : multi_array_type is nullptr.");
diff --git a/include/boost/numeric/ublas/tensor/subtensor_utility.hpp b/include/boost/numeric/ublas/tensor/subtensor_utility.hpp
index 4c38be404..4b203c138 100644
--- a/include/boost/numeric/ublas/tensor/subtensor_utility.hpp
+++ b/include/boost/numeric/ublas/tensor/subtensor_utility.hpp
@@ -37,10 +37,10 @@ namespace boost::numeric::ublas::detail {
  * @param[in] spans vector of spans of the subtensor
 */
 template<integral size_type, class Spans>
-auto compute_span_strides(std::vector<size_type> const& strides, Spans const& spans)
+auto to_span_strides(std::vector<size_type> const& strides, Spans const& spans)
 {
   if(strides.size() != spans.size())
-    throw std::runtime_error("Error in boost::numeric::ublas::subtensor::compute_span_strides(): tensor strides.size() != spans.size()");
+    throw std::runtime_error("Error in boost::numeric::ublas::subtensor::to_span_strides(): tensor strides.size() != spans.size()");
 
   auto span_strides = std::vector<size_type>(spans.size());
 
@@ -60,7 +60,7 @@ auto compute_span_strides(std::vector<size_type> const& strides, Spans const& sp
  * @param[in] spans vector of spans of the subtensor
 */
 template<integral Size, class Spans>
-auto compute_offset(std::vector<Size> const& strides, Spans const& spans)
+auto to_offset(std::vector<Size> const& strides, Spans const& spans)
 {
   if(strides.size() != spans.size())
     throw std::runtime_error("Error in boost::numeric::ublas::subtensor::offset(): tensor strides.size() != spans.size()");
@@ -77,7 +77,7 @@ auto compute_offset(std::vector<Size> const& strides, Spans const& spans)
  * @param[in] spans vector of spans of the subtensor
  */
 template<class spans_type>
-auto compute_extents(spans_type const& spans)
+auto to_extents(spans_type const& spans)
 {
   using extents_t  = extents<>;
   using base_type  = typename extents_t::base_type;
diff --git a/include/boost/numeric/ublas/tensor/tensor/tensor_static_rank.hpp b/include/boost/numeric/ublas/tensor/tensor/tensor_static_rank.hpp
index fbb5074db..3a9205480 100644
--- a/include/boost/numeric/ublas/tensor/tensor/tensor_static_rank.hpp
+++ b/include/boost/numeric/ublas/tensor/tensor/tensor_static_rank.hpp
@@ -111,6 +111,21 @@ template<class V, class L, std::size_t N>
   {
   }
 
+  /** @brief Constructs a tensor_core with a \c shape and initial value
+     *
+     * @code auto t = tensor<float>(extents<>{4,3,2},5); @endcode
+     *
+     * @param i initial tensor_core with this value
+     */
+  inline tensor_core (extents_type e, value_type i)
+    : tensor_expression_type<self_type>{}
+    , _extents(std::move(e))
+    , _strides(to_strides(_extents,layout_type{}))
+    , _container(product(_extents))
+  {
+    std::fill(begin(),end(),i);
+  }
+
   /** @brief Constructs a tensor_core with a \c shape and initiates it with one-dimensional data
      *
      * @code auto t = tensor<float,3>(extents<>{3,4,2},std::vector<float>(3*4*2,1.f)); @endcode
diff --git a/test/tensor/Jamfile b/test/tensor/Jamfile
index 9ee07cd49..ffeb7173c 100644
--- a/test/tensor/Jamfile
+++ b/test/tensor/Jamfile
@@ -32,12 +32,9 @@ explicit unit_test_framework ;
 
 test-suite boost-ublas-tensor-test
     :
-<<<<<<< HEAD
     [ run test_access.cpp
           test_algorithms.cpp
           test_einstein_notation.cpp
-          test_subtensor.cpp
-          test_subtensor_utility.cpp
           test_expression.cpp
           test_expression_evaluation.cpp
           test_extents_dynamic.cpp
@@ -65,6 +62,8 @@ test-suite boost-ublas-tensor-test
           test_static_tensor.cpp
           test_static_tensor_matrix_vector.cpp
           test_strides.cpp
+          test_subtensor.cpp
+          test_subtensor_utility.cpp          
           test_tensor.cpp
           test_tensor_matrix_vector.cpp
           unit_test_framework
diff --git a/test/tensor/test_access.cpp b/test/tensor/test_access.cpp
index dd0b08607..35fcd5938 100644
--- a/test/tensor/test_access.cpp
+++ b/test/tensor/test_access.cpp
@@ -209,8 +209,8 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_compute_single_index_static_rank, layout_
     auto const& w    = ub::to_strides(n,layout_t{});
     auto const& i    = std::get<I>(multi_index);
     auto const& jref = std::get<I>(index);
-    constexpr auto r = std::get<I>(ranks);
     mp::mp_for_each<mp::mp_iota_c<std::size(i)>>( [&]( auto K ) {
+      constexpr auto r = std::get<I>(ranks);
       auto const& ii = std::get<K>(i);
       auto const  j  = ub::detail::compute_single_index<r>(ii.begin(), ii.end() , w.begin());
       BOOST_CHECK(j < prodn(n));
@@ -226,35 +226,27 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_compute_multi_index, layout_t,  layout_ty
   namespace mp = boost::mp11;
 
   constexpr auto is_first_order = std::is_same_v<layout_t,ub::first_order>;
-    constexpr auto const& index = is_first_order ? indexf : indexl;
-
-    for(auto k = 0u; k < index.size(); ++k){
-        auto const& n = shapes[k];
-        auto const& w    = ub::to_strides(n,layout_t{});
-        auto const& iref = multi_index[k];
-        auto const& jref = index[k];
-        for(auto kk = 0u; kk < iref.size(); ++kk){
-            auto const  jj = jref[kk];
-            auto const& ii = iref[kk];
-            auto i = multi_index_t(w.size());
-            ub::detail::compute_multi_index(jj, w.begin(), w.end(), i.begin(), layout_t{});
-//            if constexpr ( is_first_order )
-//                detail::compute_multi_index_first(jj, w.begin(), w.end(), i.begin());
-//            else
-//                detail::compute_multi_index_last (jj, w.begin(), w.end(), i.begin());
-
-//            std::cout << "j= " << jj << std::endl;
-//            std::cout << "i= [ "; for(auto iii : i) std::cout << iii << " "; std::cout << "];" << std::endl;
-//            std::cout << "ii_ref = [ "; for(auto iii : ii) std::cout << iii << " "; std::cout << "];" << std::endl;
-//            std::cout << "n= [ "; for(auto iii : n) std::cout << iii << " "; std::cout << "];" << std::endl;
-//            std::cout << "w= [ "; for(auto iii : w) std::cout << iii << " "; std::cout << "];" << std::endl;
-//            std::cout << std::endl;
-
-
+  constexpr auto const& index = is_first_order ? indexf : indexl;
 
-            BOOST_CHECK ( std::equal(i.begin(),i.end(),ii.begin()) ) ;
-        }
+  for(auto k = 0u; k < index.size(); ++k){
+    auto const& n = shapes[k];
+    auto const& w    = ub::to_strides(n,layout_t{});
+    auto const& iref = multi_index[k];
+    auto const& jref = index[k];
+    for(auto kk = 0u; kk < iref.size(); ++kk){
+      auto const  jj = jref[kk];
+      auto const& ii = iref[kk];
+      auto i = multi_index_t(w.size());
+      ub::detail::compute_multi_index(jj, w.begin(), w.end(), i.begin(), layout_t{});
+//      std::cout << "j= " << jj << std::endl;
+//      std::cout << "i= [ "; for(auto iii : i) std::cout << iii << " "; std::cout << "];" << std::endl;
+//      std::cout << "ii_ref = [ "; for(auto iii : ii) std::cout << iii << " "; std::cout << "];" << std::endl;
+//      std::cout << "n= [ "; for(auto iii : n) std::cout << iii << " "; std::cout << "];" << std::endl;
+//      std::cout << "w= [ "; for(auto iii : w) std::cout << iii << " "; std::cout << "];" << std::endl;
+//      std::cout << std::endl;
+      BOOST_CHECK ( std::equal(i.begin(),i.end(),ii.begin()) ) ;
     }
+  }
 }
 
 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_compute_multi_index_static_rank, layout_t,  layout_types, fixture )
@@ -270,12 +262,12 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_compute_multi_index_static_rank, layout_t
         auto const& n    = std::get<I>(shapes);
         auto const& iref = std::get<I>(multi_index);
         auto const& jref = std::get<I>(index);
-        auto const& w    = ub::to_strides(n,layout_t{});
-        constexpr auto r = std::get<I>(ranks);
+        auto const& w    = ub::to_strides(n,layout_t{});        
         mp::mp_for_each<mp::mp_iota_c<std::size(iref)>>( [&]( auto K ) {
             auto const  jj = std::get<K>(jref);
             auto const& ii = std::get<K>(iref);
             auto         i = multi_index_t(w.size());
+            constexpr auto r = std::get<I>(ranks);
             ub::detail::compute_multi_index<r>(jj, w.begin(), w.end(), i.begin(), layout_t{});
             BOOST_CHECK ( std::equal(i.begin(),i.end(),ii.begin()) ) ;
         });
diff --git a/test/tensor/test_functions.cpp b/test/tensor/test_functions.cpp
index 920fdac9a..4e8e4f464 100644
--- a/test/tensor/test_functions.cpp
+++ b/test/tensor/test_functions.cpp
@@ -1,4 +1,4 @@
-//
+//
 //  Copyright (c) 2018-2020, Cem Bassoy, cem.bassoy@gmail.com
 //  Copyright (c) 2019-2020, Amit Singh, amitsingh19975@gmail.com
 //
@@ -25,127 +25,126 @@
 BOOST_AUTO_TEST_SUITE ( test_tensor_functions)
 
 
-using test_types = zip<int,float,std::complex<float>>::with_t<boost::numeric::ublas::layout::first_order, boost::numeric::ublas::layout::last_order>;
+using test_types = zip<float>::with_t<boost::numeric::ublas::layout::first_order, boost::numeric::ublas::layout::last_order>;
 
+// int,float,std::complex<float>
 //using test_types = zip<int>::with_t<boost::numeric::ublas::layout::first_order>;
 
 
 struct fixture
 {
-    using dynamic_extents_type = boost::numeric::ublas::extents<>;
-    fixture()
-      : extents {
-          dynamic_extents_type{1,1}, // 1
-          dynamic_extents_type{2,3}, // 2
-          dynamic_extents_type{2,3,1}, // 3
-          dynamic_extents_type{4,2,3}, // 4
-          dynamic_extents_type{4,2,3,5}} // 5
-    {
-    }
+  using shape = boost::numeric::ublas::extents<>;
 
-    std::vector<dynamic_extents_type> extents;
+  const std::vector<shape> extents
+    {
+      shape{1,1}, // 1
+      shape{2,1}, // 2
+      shape{1,2}, // 3
+      shape{2,3}, // 4
+      shape{2,3,1}, // 5
+      shape{1,2,3}, // 6
+      shape{3,1,2}, // 7
+      shape{4,2,3}, // 8
+      shape{4,2,3,1}, // 9
+      shape{4,2,3,5} // 10
+    };
 };
 
 
 
 
-BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_vector, value,  test_types, fixture )
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_vector, pair,  test_types, fixture )
 {
-    namespace ublas = boost::numeric::ublas;
-    using value_type   = typename value::first_type;
-    using layout_type  = typename value::second_type;
-    using tensor_type  = ublas::tensor_dynamic<value_type,layout_type>;
-    using vector_type  = typename tensor_type::vector_type;
-
-
-    for(auto const& n : extents){
+  namespace ublas = boost::numeric::ublas;
+  using value   = typename pair::first_type;
+  using layout  = typename pair::second_type;
+  using tensor  = ublas::tensor_dynamic<value,layout>;
+  using vector  = typename tensor::vector_type;
 
-        auto a = tensor_type(n, value_type{2});
 
-        std::cout << "a=" << a << std::endl;
+  for(auto const& n : extents){
 
-        for(auto m = 0u; m < ublas::size(n); ++m){
+    auto v = value(2);
+    auto a = tensor(n, v);
 
-            auto b = vector_type  (n[m], value_type{1} );
+    for(auto m = 0u; m < ublas::size(n); ++m){
 
-            std::cout << "b=" << tensor_type(b) << std::endl;
+      auto b = vector  (n[m], value{1} );
 
-            std::cout << "m=" << m << std::endl;
+      auto c = ublas::prod(a, b, m+1);
 
-            auto c = ublas::prod(a, b, m+1);
+      auto vv = v * value(n[m]);
 
-            std::cout << "c=" << tensor_type(c) << std::endl;
+      BOOST_CHECK( std::all_of( c.begin(), c.end() , [vv](auto cc){ return cc == vv; } ) );
 
-            for(auto i = 0u; i < c.size(); ++i)
-                BOOST_CHECK_EQUAL( c[i] , value_type( static_cast< inner_type_t<value_type> >(n[m]) ) * a[i] );
-
-        }
     }
-  auto n = extents[4];
-  auto a = tensor_type(n, value_type{2});
-  auto b = vector_type(n[0], value_type{1});
+  }
+//  auto n = extents[4];
+//  auto a = tensor_type(n, value_type{2});
+//  auto b = vector_type(n[0], value_type{1});
 
-  auto empty = vector_type{};
+//  auto empty = vector_type{};
 
-  BOOST_CHECK_THROW(prod(a, b, 0), std::length_error);
-  BOOST_CHECK_THROW(prod(a, b, 9), std::length_error);
-  BOOST_CHECK_THROW(prod(a, empty, 2), std::length_error);
+//  BOOST_CHECK_THROW(prod(a, b, 0), std::length_error);
+//  BOOST_CHECK_THROW(prod(a, b, 9), std::length_error);
+//  BOOST_CHECK_THROW(prod(a, empty, 2), std::length_error);
 
 }
 
 BOOST_AUTO_TEST_CASE( test_tensor_prod_vector_exception )
 {
-//    namespace ublas = boost::numeric::ublas;
-//    using value_type   = float;
-//    using layout_type  = ublas::layout::first_order;
-//    using d_tensor_type  = ublas::tensor_dynamic<value_type,layout_type>;
-//    using vector_type  = typename d_tensor_type::vector_type;
+    namespace ublas = boost::numeric::ublas;
+    using value  = float;
+    using layout  = ublas::layout::first_order;
+    using tensor  = ublas::tensor_dynamic<value,layout>;
+    using vector  = typename tensor::vector_type;
 
-//    auto t1 = d_tensor_type{ublas::extents<>{},1.f};
-//    auto v1 = vector_type{3,value_type{1}};
+    auto t1 = tensor{ublas::extents<>{},1.f};
+    auto v1 = vector{3,value{1}};
 
-//    BOOST_REQUIRE_THROW(prod(t1,v1,0),std::length_error);
-//    BOOST_REQUIRE_THROW(prod(t1,v1,1),std::length_error);
-//    BOOST_REQUIRE_THROW(prod(t1,v1,3),std::length_error);
+    BOOST_REQUIRE_THROW(prod(t1,v1,0),std::length_error);
+    BOOST_REQUIRE_THROW(prod(t1,v1,1),std::length_error);
+    BOOST_REQUIRE_THROW(prod(t1,v1,3),std::length_error);
 }
 
 
 
 
-BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_matrix, value,  test_types, fixture )
+BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_matrix, pair,  test_types, fixture )
 {
-    namespace ublas = boost::numeric::ublas;
-    using value_type   = typename value::first_type;
-    using layout_type  = typename value::second_type;
-    using tensor_type  = ublas::tensor_dynamic<value_type,layout_type>;
-    using matrix_type  = typename tensor_type::matrix_type;
+  namespace ublas = boost::numeric::ublas;
+  using value   = typename pair::first_type;
+  using layout  = typename pair::second_type;
+  using tensor  = ublas::tensor_dynamic<value,layout>;
+  using matrix  = typename tensor::matrix_type;
 
 
-    for(auto const& n : extents) {
+  for(auto const& n : extents) {
 
-        auto a = tensor_type(n, value_type{2});
+    auto v = value{2};
+    auto a = tensor(n, v);
 
-        for(auto m = 0u; m < ublas::size(n); ++m){
+    for(auto m = 0u; m < ublas::size(n); ++m){
 
-            auto b  = matrix_type  ( n[m], n[m], value_type{1} );
+      auto b  = matrix  ( n[m], n[m], value{1} );
 
-            auto c = ublas::prod(a, b, m+1);
+      auto c = ublas::prod(a, b, m+1);
 
-            for(auto i = 0u; i < c.size(); ++i)
-                BOOST_CHECK_EQUAL( c[i] , value_type( static_cast< inner_type_t<value_type> >(n[m]) ) * a[i] );
+      auto vv = v * value(n[m]);
+      BOOST_CHECK( std::all_of( c.begin(), c.end() , [vv](auto cc){ return cc == vv; } ) );
 
-        }
     }
+  }
 
-  auto n = extents[4];
-  auto a = tensor_type(n, value_type{2});
-  auto b = matrix_type(n[0], n[0], value_type{1});
+  //  auto n = extents[4];
+  //  auto a = tensor_type(n, value_type{2});
+  //  auto b = matrix_type(n[0], n[0], value_type{1});
 
-  auto empty = matrix_type{};
+  //  auto empty = matrix_type{};
 
-  BOOST_CHECK_THROW(prod(a, b, 0), std::length_error);
-  BOOST_CHECK_THROW(prod(a, b, 9), std::length_error);
-  BOOST_CHECK_THROW(prod(a, empty, 2), std::invalid_argument);
+  //  BOOST_CHECK_THROW(prod(a, b, 0), std::length_error);
+  //  BOOST_CHECK_THROW(prod(a, b, 9), std::length_error);
+  //  BOOST_CHECK_THROW(prod(a, empty, 2), std::invalid_argument);
 }
 
 BOOST_AUTO_TEST_CASE( test_tensor_prod_matrix_exception )
@@ -196,8 +195,8 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_tensor_1, value,  test_types,
             for(auto i = 0ul; i < q; ++i)
                 acc *= value_type( static_cast< inner_type_t<value_type> >( a.extents().at(phi.at(i)-1) ) );
 
-            for(auto i = 0ul; i < c.size(); ++i)
-                BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
+            auto vv = acc * a[0] * b[0];
+            BOOST_CHECK( std::all_of(c.begin(), c.end() , [vv](auto cc){ return cc == vv; } ) );
 
         }
     }
@@ -296,8 +295,13 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_tensor_2, value,  test_types,
                 for(auto i = 0ul; i < q; ++i)
                     acc *= value_type( static_cast< inner_type_t<value_type> >( a.extents().at(phia.at(i)-1) ) );
 
-                for(auto i = 0ul; i < c.size(); ++i)
-                    BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
+
+                auto vv = acc * a[0] * b[0];
+                BOOST_CHECK( std::all_of(c.begin(), c.end() , [vv](auto cc){ return cc == vv; } ) );
+
+
+//                for(auto i = 0ul; i < c.size(); ++i)
+//                    BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
 
             }
 
diff --git a/test/tensor/test_multiplication.cpp b/test/tensor/test_multiplication.cpp
index c5ed51e5f..0f9700915 100644
--- a/test/tensor/test_multiplication.cpp
+++ b/test/tensor/test_multiplication.cpp
@@ -25,7 +25,9 @@
 BOOST_AUTO_TEST_SUITE (test_tensor_contraction)
 
 
-using test_types = zip<int,float,std::complex<float>>::with_t<boost::numeric::ublas::layout::first_order, boost::numeric::ublas::layout::last_order>;
+//using test_types = zip<int,float,std::complex<float>>::with_t<boost::numeric::ublas::layout::first_order, boost::numeric::ublas::layout::last_order>;
+
+using test_types = zip<float>::with_t<boost::numeric::ublas::layout::first_order, boost::numeric::ublas::layout::last_order>;
 
 //using test_types = zip<int>::with_t<boost::numeric::ublas::layout::first_order>;
 
@@ -62,38 +64,38 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_mtv, value,  test_types, fixture )
 
   for(auto const& na : extents) {
 
-    if(ublas::size(na) > 2)
+    if(ublas::is_scalar(na) || ublas::is_vector(na) || ublas::is_tensor(na))
       continue;
 
+    auto const n1 = na[0];
+    auto const n2 = na[1];
+
     auto a = vector_t(ublas::product(na), value_t{2});
     auto wa = ublas::to_strides(na,layout_t{});
-    for(auto m = std::size_t{0}; m < ublas::size(na); ++m){
+
+    for(auto m = std::size_t{0}; m < 2; ++m){
+
       auto nb = extents_t {na[m],std::size_t{1}};
       auto wb = ublas::to_strides(nb,layout_t{});
       auto b  = vector_t  (ublas::product(nb), value_t{1} );
 
-      auto nc_base = extents_base_t(std::max(std::size_t{ublas::size(na)-1u}, std::size_t{2}), 1);
+      // [n1 n2 1 ... 1] x1 [n1 1] -> [n2 1 ... 1]
+      // [n1 n2 1 ... 1] x2 [n2 1] -> [n1 1 ... 1]
 
-      for(auto i = 0ul, j = 0ul; i < ublas::size(na); ++i)
-        if(i != m)
-          nc_base[j++] = na[i];
+      auto nc_base = extents_base_t(std::max(std::size_t(ublas::size(na)-1u), std::size_t{2}), 1);
+      nc_base[0] = m==0 ? n2 : n1;
 
       auto nc = extents_t (nc_base);
       auto wc = ublas::to_strides(nc,layout_t{});
       auto c  = vector_t  (ublas::product(nc), value_t{0});
 
-      ublas::detail::recursive::mtv(
-        m,
-        c.data(), nc.data(), wc.data(),
-        a.data(), na.data(), wa.data(),
-        b.data());
+      assert( (m==0u) || (m==1u));
+
+      ublas::detail::recursive::mtv(m, n1,n2, c.data(), size_t(1), a.data(), wa[0], wa[1], b.data(), size_t(1));
 
       auto v = value_t(na[m]);
       BOOST_CHECK(std::equal(c.begin(),c.end(),a.begin(), [v](auto cc, auto aa){return cc == v*aa;}));
 
-//      for(auto i = 0u; i < c.size(); ++i)
-//        BOOST_CHECK_EQUAL( c[i] , value_t( static_cast< inner_type_t<value_t> >(na[m]) ) * a[i] );
-
     }
   }
 }
@@ -132,10 +134,6 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_mtm, value,  test_types, fixture )
     auto v = value_t(na[1])*a[0];
     BOOST_CHECK(std::all_of(c.begin(),c.end(), [v](auto cc){return cc == v;}));
 
-//    for(auto i = 0u; i < c.size(); ++i)
-//      BOOST_CHECK_EQUAL( c[i] , value_t( static_cast< inner_type_t<value_t> >(na[1]) ) * a[0] );
-
-
   }
 }
 
@@ -176,8 +174,6 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttv, value,  test_types, fixture )
       auto v = value_t(na[m]);
       BOOST_CHECK(std::equal(c.begin(),c.end(),a.begin(), [v](auto cc, auto aa){return cc == v*aa;}));
 
-//      for(auto i = 0u; i < c.size(); ++i)
-//        BOOST_CHECK_EQUAL( c[i] , value_t(na[m]) * a[i] );
 
     }
   }
diff --git a/test/tensor/test_subtensor.cpp b/test/tensor/test_subtensor.cpp
index f3db7334d..d091f4375 100644
--- a/test/tensor/test_subtensor.cpp
+++ b/test/tensor/test_subtensor.cpp
@@ -128,8 +128,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( subtensor_ctor2_test, value,  test_types )
     auto A    = tensor_type{1,2};
     auto Asub = subtensor_type( A, 0, 1  );
 
-    BOOST_CHECK_EQUAL( Asub.span_strides().at(0), 1 );
-    BOOST_CHECK_EQUAL( Asub.span_strides().at(1), 1 );
+    BOOST_CHECK_EQUAL( Asub.span_strides().at(0), A.strides().at(0) );
+    BOOST_CHECK_EQUAL( Asub.span_strides().at(1), A.strides().at(1) );
 
     BOOST_CHECK_EQUAL( Asub.strides().at(0), 1 );
     BOOST_CHECK_EQUAL( Asub.strides().at(1), 1 );
diff --git a/test/tensor/test_subtensor_utility.cpp b/test/tensor/test_subtensor_utility.cpp
index b96fb5b57..c1fab9fa1 100644
--- a/test/tensor/test_subtensor_utility.cpp
+++ b/test/tensor/test_subtensor_utility.cpp
@@ -293,12 +293,12 @@ BOOST_FIXTURE_TEST_CASE( extents_test, fixture_span_vector_shape )
 {
   namespace ublas = boost::numeric::ublas;
 
-  BOOST_CHECK ( std::equal( ublas::begin(std::get<0>(reference_)), ublas::begin(std::get<0>(reference_)),  ublas::begin(ublas::detail::compute_extents( std::get<0>(span_vectors_) ) ) ) );
-  BOOST_CHECK ( std::equal( ublas::begin(std::get<1>(reference_)), ublas::begin(std::get<1>(reference_)),  ublas::begin(ublas::detail::compute_extents( std::get<1>(span_vectors_) ) ) ) );
-  BOOST_CHECK ( std::equal( ublas::begin(std::get<2>(reference_)), ublas::begin(std::get<2>(reference_)),  ublas::begin(ublas::detail::compute_extents( std::get<2>(span_vectors_) ) ) ) );
-  BOOST_CHECK ( std::equal( ublas::begin(std::get<3>(reference_)), ublas::begin(std::get<3>(reference_)),  ublas::begin(ublas::detail::compute_extents( std::get<3>(span_vectors_) ) ) ) );
-  BOOST_CHECK ( std::equal( ublas::begin(std::get<4>(reference_)), ublas::begin(std::get<4>(reference_)),  ublas::begin(ublas::detail::compute_extents( std::get<4>(span_vectors_) ) ) ) );
-  BOOST_CHECK ( std::equal( ublas::begin(std::get<5>(reference_)), ublas::begin(std::get<5>(reference_)),  ublas::begin(ublas::detail::compute_extents( std::get<5>(span_vectors_) ) ) ) );
+  BOOST_CHECK ( std::equal( ublas::begin(std::get<0>(reference_)), ublas::begin(std::get<0>(reference_)),  ublas::begin(ublas::detail::to_extents( std::get<0>(span_vectors_) ) ) ) );
+  BOOST_CHECK ( std::equal( ublas::begin(std::get<1>(reference_)), ublas::begin(std::get<1>(reference_)),  ublas::begin(ublas::detail::to_extents( std::get<1>(span_vectors_) ) ) ) );
+  BOOST_CHECK ( std::equal( ublas::begin(std::get<2>(reference_)), ublas::begin(std::get<2>(reference_)),  ublas::begin(ublas::detail::to_extents( std::get<2>(span_vectors_) ) ) ) );
+  BOOST_CHECK ( std::equal( ublas::begin(std::get<3>(reference_)), ublas::begin(std::get<3>(reference_)),  ublas::begin(ublas::detail::to_extents( std::get<3>(span_vectors_) ) ) ) );
+  BOOST_CHECK ( std::equal( ublas::begin(std::get<4>(reference_)), ublas::begin(std::get<4>(reference_)),  ublas::begin(ublas::detail::to_extents( std::get<4>(span_vectors_) ) ) ) );
+  BOOST_CHECK ( std::equal( ublas::begin(std::get<5>(reference_)), ublas::begin(std::get<5>(reference_)),  ublas::begin(ublas::detail::to_extents( std::get<5>(span_vectors_) ) ) ) );
 
 }
 
@@ -314,35 +314,35 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( offset_test, layout, test_types, fixture_span_
 	{
     auto s = std::get<0>(span_vectors_);
     auto w = ublas::to_strides( std::get<0>(extents_), layout{} );
-    auto o = ublas::detail::compute_offset(w,s);
+    auto o = ublas::detail::to_offset(w,s);
 		BOOST_CHECK_EQUAL( o, 0  );
 	}
 
 	{
     auto s = std::get<1>(span_vectors_);
     auto w = ublas::to_strides( std::get<1>(extents_), layout{} );
-    auto o = ublas::detail::compute_offset(w,s);
+    auto o = ublas::detail::to_offset(w,s);
 		BOOST_CHECK_EQUAL( o, 0  );
 	}
 
 	{
     auto s = std::get<2>(span_vectors_);
     auto w = ublas::to_strides( std::get<2>(extents_), layout{} );
-    auto o = ublas::detail::compute_offset(w,s);
+    auto o = ublas::detail::to_offset(w,s);
 		BOOST_CHECK_EQUAL( o, 0  );
 	}
 
 	{
     auto s = std::get<3>(span_vectors_);
     auto w = ublas::to_strides( std::get<3>(extents_), layout{} );
-    auto o = ublas::detail::compute_offset(w,s);
+    auto o = ublas::detail::to_offset(w,s);
 		BOOST_CHECK_EQUAL( o, s[0].first()*w[0] + s[1].first()*w[1]  );
 	}
 
 	{
     auto s = std::get<4>(span_vectors_);
     auto w = ublas::to_strides( std::get<4>(extents_), layout{} );
-    auto o = ublas::detail::compute_offset(w,s);
+    auto o = ublas::detail::to_offset(w,s);
 		BOOST_CHECK_EQUAL( o, s[0].first()*w[0] + s[1].first()*w[1] + s[2].first()*w[2]   );
 	}
 
@@ -350,7 +350,7 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( offset_test, layout, test_types, fixture_span_
 	{
     auto s = std::get<5>(span_vectors_);
     auto w = ublas::to_strides( std::get<5>(extents_), layout{} );
-    auto o = ublas::detail::compute_offset(w,s);
+    auto o = ublas::detail::to_offset(w,s);
     BOOST_CHECK_EQUAL( o, s[0].first()*w[0] + s[1].first()*w[1] + s[2].first()*w[2] +  s[3].first()*w[3] );
 	}