Fixed initialization of ordering projection tensors

Signed-off-by: Dmitry I. Lyakh <quant4me@gmail.com>

src/numerics/functor_init_proj.cpp

 /** ExaTN::Numerics: Tensor Functor: Initialization of Ordering Projection tensors
-REVISION: 2021/02/16
+REVISION: 2021/09/21
 
 Copyright (C) 2018-2021 Dmitry I. Lyakh (Liakh)
 Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
@@ -35,7 +35,7 @@ int FunctorInitProj::apply(talsh::Tensor & local_tensor) //tensor slice (in gene
 
  assert(rank%2 == 0); //only even-order tensors
 
- auto init_proj = [&](auto * tensor_body){
+ auto init_proj_strong = [&](auto * tensor_body){
   std::vector<DimOffset> bas(rank);
   for(unsigned int i = 0; i < rank; ++i) bas[i] = offsets[i]; //tensor slice dimension base offsets
   std::vector<DimExtent> ext(rank);
@@ -53,29 +53,71 @@ int FunctorInitProj::apply(talsh::Tensor & local_tensor) //tensor slice (in gene
   return 0;
  };
 
+ auto init_proj_weak = [&](auto * tensor_body){
+  std::vector<DimOffset> bas(rank);
+  for(unsigned int i = 0; i < rank; ++i) bas[i] = offsets[i]; //tensor slice dimension base offsets
+  std::vector<DimExtent> ext(rank);
+  for(unsigned int i = 0; i < rank; ++i) ext[i] = extents[i]; //tensor slice dimension extents
+  TensorRange rng(bas,ext);
+  bool not_over = true;
+  while(not_over){
+   if(rng.nondecreasingOrderDiag()){
+    tensor_body[rng.localOffset()] = 1.0;
+   }else{
+    tensor_body[rng.localOffset()] = 0.0;
+   }
+   not_over = rng.next();
+  }
+  return 0;
+ };
+
  auto access_granted = false;
  {//Try REAL32:
   float * body;
   access_granted = local_tensor.getDataAccessHost(&body);
-  if(access_granted) return init_proj(body);
+  if(access_granted){
+   if(weak_ordering_){
+    return init_proj_weak(body);
+   }else{
+    return init_proj_strong(body);
+   }
+  }
  }
 
  {//Try REAL64:
   double * body;
   access_granted = local_tensor.getDataAccessHost(&body);
-  if(access_granted) return init_proj(body);
+  if(access_granted){
+   if(weak_ordering_){
+    return init_proj_weak(body);
+   }else{
+    return init_proj_strong(body);
+   }
+  }
  }
 
  {//Try COMPLEX32:
   std::complex<float> * body;
   access_granted = local_tensor.getDataAccessHost(&body);
-  if(access_granted) return init_proj(body);
+  if(access_granted){
+   if(weak_ordering_){
+    return init_proj_weak(body);
+   }else{
+    return init_proj_strong(body);
+   }
+  }
  }
 
  {//Try COMPLEX64:
   std::complex<double> * body;
   access_granted = local_tensor.getDataAccessHost(&body);
-  if(access_granted) return init_proj(body);
+  if(access_granted){
+   if(weak_ordering_){
+    return init_proj_weak(body);
+   }else{
+    return init_proj_strong(body);
+   }
+  }
  }
 
  std::cout << "#ERROR(exatn::numerics::FunctorInitProj): Unknown data kind in talsh::Tensor!" << std::endl;

src/numerics/functor_init_proj.hpp

 /** ExaTN::Numerics: Tensor Functor: Initialization of Ordering Projection tensors
-REVISION: 2021/02/16
+REVISION: 2021/09/21
 
 Copyright (C) 2018-2021 Dmitry I. Lyakh (Liakh)
 Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
@@ -8,7 +8,7 @@ Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
  (A) This tensor functor (method) is used to initialize an ordering projection tensor.
      An ordering projection tensor is an even-order tensor with the following elements:
      - Element = 1 if: The first half and the second half of the indices are the same and
-                       both are in a monotonically increasing order;
+                       both are in a monotonically increasing (or non-decreasing) order;
      - Elelent = 0 otherwise.
 **/
 
@@ -33,7 +33,9 @@ namespace numerics{
 class FunctorInitProj: public talsh::TensorFunctor<Identifiable>{
 public:
 
- FunctorInitProj() = default;
+ FunctorInitProj(bool weak_ordering = false): //strong VS weak index ordering
+  weak_ordering_(weak_ordering)
+ {}
 
  virtual ~FunctorInitProj() = default;
 
@@ -61,6 +63,7 @@ public:
 
 private:
 
+ bool weak_ordering_;
 };
 
 } //namespace numerics

src/numerics/tensor_range.hpp

 /** ExaTN::Numerics: Tensor range
-REVISION: 2021/02/16
+REVISION: 2021/09/21
 
 Copyright (C) 2018-2021 Dmitry I. Lyakh (Liakh)
 Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
@@ -108,6 +108,14 @@ public:
      each is in a monotonically decreasing order. **/
  inline bool decreasingOrderDiag() const;
 
+ /** Tests whether both halves of the indices of the current multi-index are the same and
+     each is in a monotonically non-decreasing order. **/
+ inline bool nondecreasingOrderDiag() const;
+
+ /** Tests whether both halves of the indices of the current multi-index are the same and
+     each is in a monotonically non-increasing order. **/
+ inline bool nonincreasingOrderDiag() const;
+
  /** Returns the flat offset produced by the current multi-index value per se. **/
  inline DimOffset localOffset() const; //little endian
 
@@ -353,7 +361,7 @@ inline bool TensorRange::increasingOrderDiag() const
  if(indeed && half_size > 0){
   for(int i = 1; i < half_size; ++i){
    if(((bases_[i] + mlndx_[i]) <= (bases_[i-1] + mlndx_[i-1])) ||
-      ((bases_[i] + mlndx_[i]) == (bases_[half_size+i] + mlndx_[half_size+i]))){
+      ((bases_[i] + mlndx_[i]) != (bases_[half_size+i] + mlndx_[half_size+i]))){
     indeed = false;
     break;
    }
@@ -371,7 +379,43 @@ inline bool TensorRange::decreasingOrderDiag() const
  if(indeed && half_size > 0){
   for(int i = 1; i < half_size; ++i){
    if(((bases_[i] + mlndx_[i]) >= (bases_[i-1] + mlndx_[i-1])) ||
-      ((bases_[i] + mlndx_[i]) == (bases_[half_size+i] + mlndx_[half_size+i]))){
+      ((bases_[i] + mlndx_[i]) != (bases_[half_size+i] + mlndx_[half_size+i]))){
+    indeed = false;
+    break;
+   }
+  }
+  indeed = indeed && ((bases_[0] + mlndx_[0]) == (bases_[half_size] + mlndx_[half_size]));
+ }
+ return indeed;
+}
+
+
+inline bool TensorRange::nondecreasingOrderDiag() const
+{
+ bool indeed = (mlndx_.size() % 2 == 0);
+ int half_size = mlndx_.size() / 2;
+ if(indeed && half_size > 0){
+  for(int i = 1; i < half_size; ++i){
+   if(((bases_[i] + mlndx_[i]) < (bases_[i-1] + mlndx_[i-1])) ||
+      ((bases_[i] + mlndx_[i]) != (bases_[half_size+i] + mlndx_[half_size+i]))){
+    indeed = false;
+    break;
+   }
+  }
+  indeed = indeed && ((bases_[0] + mlndx_[0]) == (bases_[half_size] + mlndx_[half_size]));
+ }
+ return indeed;
+}
+
+
+inline bool TensorRange::nonincreasingOrderDiag() const
+{
+ bool indeed = (mlndx_.size() % 2 == 0);
+ int half_size = mlndx_.size() / 2;
+ if(indeed && half_size > 0){
+  for(int i = 1; i < half_size; ++i){
+   if(((bases_[i] + mlndx_[i]) > (bases_[i-1] + mlndx_[i-1])) ||
+      ((bases_[i] + mlndx_[i]) != (bases_[half_size+i] + mlndx_[half_size+i]))){
     indeed = false;
     break;
    }