introduced is_contiguous(Range) + use range_traits instead of directly reading Range::order

Author: evaleev

btas/optimize/contract.h

@@ -17,7 +17,7 @@ template<typename _T, class _TensorA, class _TensorB, class _TensorC,
 void contract_211(const _T& alpha, const _TensorA& A, const btas::DEFAULT::index<_UA>& aA, const _TensorB& B, const btas::DEFAULT::index<_UB>& aB,
                   const _T& beta, _TensorC& C, const btas::DEFAULT::index<_UC>& aC, const bool conjgA, const bool conjgB) {
   assert(aA.size() == 2 && aB.size() == 1 && aC.size() == 1);
-  assert(A.range().ordinal().contiguous() && B.range().ordinal().contiguous() && C.range().ordinal().contiguous());
+  assert(is_contiguous(A.range()) && is_contiguous(B.range()) && is_contiguous(C.range()));
   if (conjgB) throw std::logic_error("complex conjugation of 1-index tensors is not considered in contract_211");
 
   const bool notrans = aB[0] == aA[1];
@@ -36,7 +36,7 @@ void contract_222(const _T& alpha, const _TensorA& A, const btas::DEFAULT::index
                   const _T& beta, _TensorC& C, const btas::DEFAULT::index<_UC>& aC, const bool conjgA, const bool conjgB) {
   // TODO we do not consider complex matrices yet.
   assert(aA.size() == 2 && aB.size() == 2 && aC.size() == 2);
-  assert(A.range().ordinal().contiguous() && B.range().ordinal().contiguous() && C.range().ordinal().contiguous());
+  assert(is_contiguous(A.range()) && is_contiguous(B.range()) && is_contiguous(C.range()));
   if (std::find(aA.begin(), aA.end(), aC.front()) != aA.end()) {
     // then multiply A * B -> C
     const bool notransA = aA.front() == aC.front();
@@ -65,10 +65,10 @@ template<typename _T, class _TensorA, class _TensorB, class _TensorC,
 void contract_323(const _T& alpha, const _TensorA& A, const btas::DEFAULT::index<_UA>& aA, const _TensorB& B, const btas::DEFAULT::index<_UB>& aB,
                   const _T& beta, _TensorC& C, const btas::DEFAULT::index<_UC>& aC, const bool conjgA, const bool conjgB) {
   assert(aA.size() == 3 && aB.size() == 2 && aC.size() == 3);
-  assert(A.range().ordinal().contiguous() && B.range().ordinal().contiguous() && C.range().ordinal().contiguous());
+  assert(is_contiguous(A.range()) && is_contiguous(B.range()) && is_contiguous(C.range()));
   if (conjgA) throw std::logic_error("complex conjugation of 3-index tensors is not considered in contract_323");
 
-  // TODO this function is limited to special cases where one of three indices of A will be replaced in C. Permuation is not considered so far.
+  // TODO this function is limited to special cases where one of three indices of A will be replaced in C. Permutation is not considered so far.
   // first idenfity which indices to be rotated
   int irot = -1;
   for (int i = 0; i != 3; ++i)
@@ -128,7 +128,7 @@ template<typename _T, class _TensorA, class _TensorB, class _TensorC,
 void contract_332(const _T& alpha, const _TensorA& A, const btas::DEFAULT::index<_UA>& aA, const _TensorB& B, const btas::DEFAULT::index<_UB>& aB,
                   const _T& beta, _TensorC& C, const btas::DEFAULT::index<_UC>& aC, const bool conjgA, const bool conjgB) {
   assert(aA.size() == 3 && aB.size() == 3 && aC.size() == 2);
-  assert(A.range().ordinal().contiguous() && B.range().ordinal().contiguous() && C.range().ordinal().contiguous());
+  assert(is_contiguous(A.range()) && is_contiguous(B.range()) && is_contiguous(C.range()));
 
   const bool back2  = aA[0] == aB[0] && aA[1] == aB[1];
   const bool front2 = aA[1] == aB[1] && aA[2] == aB[2];

btas/range.h

@@ -1181,12 +1181,22 @@ namespace btas {
                           std::rbegin(r2_extent));
     }
 
+    /// Tests whether a range is contiguous, i.e. whether its ordinal values form a contiguous range
+
+    /// \param range a Range
+    /// \return true if \p range is contiguous
+    template <CBLAS_ORDER _Order,
+        typename _Index,
+        typename _Ordinal>
+    inline bool is_contiguous(const RangeNd<_Order, _Index, _Ordinal>& range) {
+      return range.ordinal().contiguous();
+    }
+
     /// Permutes a Range
 
-    /// permutes the axes using permutation \c p={p[0],p[1],...} specified in the preimage ("from") convention;
-    /// for example, after this call \c lobound()[p[i]] will return the value originally
-    /// returned by \c lobound()[i]
-    /// \param perm a sequence specifying from-permutation of the axes
+    /// permutes the dimensions using permutation \c p = {p[0], p[1], ... }; for example, if \c lobound() initially returned
+    /// {lb[0], lb[1], ... }, after this call \c lobound() will return {lb[p[0]], lb[p[1]], ...}.
+    /// \param perm an array specifying permutation of the dimensions
     template <CBLAS_ORDER _Order,
               typename _Index,
               typename _Ordinal,
@@ -1216,10 +1226,9 @@ namespace btas {
 
     /// Permutes a Range
 
-    /// permutes the axes using permutation \c p={p[0],p[1],...} specified in the preimage ("from") convention;
-    /// for example, after this call \c lobound()[p[i]] will return the value originally
-    /// returned by \c lobound()[i]
-    /// \param perm an initializer list specifying from-permutation of the axes
+    /// permutes the axes using permutation \c p = {p[0], p[1], ... }; for example, if \c lobound() initially returned
+    /// {lb[0], lb[1], ... }, after this call \c lobound() will return {lb[p[0]], lb[p[1]], ...} .
+    /// \param perm an array specifying permutation of the axes
     template <CBLAS_ORDER _Order,
               typename _Index,
               typename _Ordinal,

btas/tensor.h

@@ -736,9 +736,9 @@ namespace btas {
   bool operator==(const _Tensor1& t1, const _Tensor2& t2) {
       using std::cbegin;
       using std::cend;
-      if (t1.range().order == t2.range().order &&
-          t1.range().ordinal().contiguous() &&
-          t2.range().ordinal().contiguous()) // plain Tensor
+      if (btas::range_traits<std::decay_t<decltype(t1.range())>>::order == btas::range_traits<std::decay_t<decltype(t2.range())>>::order &&
+          is_contiguous(t1.range()) &&
+          is_contiguous(t2.range())) // plain Tensor
         return congruent(t1.range(), t2.range()) && std::equal(cbegin(t1.storage()),
                                                                cend(t1.storage()),
                                                                cbegin(t2.storage()));