1. Added capability (via constructor) to repack tensor data into a different...

1. Added capability (via constructor) to repack tensor data into a different packed symmetric layout without doing the normal symmetrization permutations that happen during sum
2. Fixed bug associated with C["ij"]=A["ij"]*B["ijkl"] where C and A are SY and B is fully NS.
3. Getting (2) to work in parallel required resolving a bug associated with desymmetrization (in the two unfold_broken_Sym functions), which were changing symmetry without changing sym_table, and as a result yielding extra mapping restrictions.

The bug fix to (3) might improve performance generally.

Makefile

@@ -6,7 +6,7 @@ EXAMPLES = dft dft_3D gemm gemm_4D scalar trace weigh_4D subworld_gemm \
            permute_multiworld strassen slice_gemm ccsd sparse_permuted_slice 
 
 TESTS = test_suite pgemm_test nonsq_pgemm_test diag_sym sym3 readwrite_test \
-        ccsdt_t3_to_t2 ccsdt_map_test multi_tsr_sym diag_ctr readall_test 
+        ccsdt_t3_to_t2 ccsdt_map_test multi_tsr_sym diag_ctr readall_test  sy_times_ns repack
 
 BENCHMARKS = nonsq_pgemm_bench bench_contraction bench_nosym_transp bench_redistribution
 

examples/dft.cxx

@@ -45,7 +45,7 @@ int test_dft(int64_t n,
 
   /*DFT.contract(std::complex<double> (1.0, 0.0), DFT, "ij", IDFT, "jk", 
                std::complex<double> (0.0, 0.0), "ik");*/
-  DFT["ik"] = DFT["ij"]*IDFT["jk"];
+  DFT["ik"] = .5*DFT["ij"]*IDFT["jk"];
 
  
   DFT.read_local(&np, &idx, &data);

src/contraction/contraction.cxx

@@ -507,21 +507,28 @@ namespace CTF_int {
             C->order, idx_C,
             &num_tot, &idx_arr);
 
+    int nA_sym[A->order];
+    if (new_contraction != NULL)
+      memcpy(nA_sym, nA->sym, sizeof(int)*nA->order);
     for (i=0; i<A->order; i++){
       if (A->sym[i] != NS){
         iA = idx_A[i];
         if (idx_arr[3*iA+1] != -1){
           if (B->sym[idx_arr[3*iA+1]] == NS ||
               idx_A[i+1] != idx_B[idx_arr[3*iA+1]+1]){
-            if (new_contraction != NULL)
-              nA->sym[i] = NS;
+            if (new_contraction != NULL){
+              nA_sym[i] = NS;
+              nA->set_sym(nA_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i;
           }
         } else {
           if (idx_arr[3*idx_A[i+1]+1] != -1){
-            if (new_contraction != NULL)
-              nA->sym[i] = NS;
+            if (new_contraction != NULL){
+              nA_sym[i] = NS;
+              nA->set_sym(nA_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i;
           }       
@@ -529,15 +536,19 @@ namespace CTF_int {
         if (idx_arr[3*iA+2] != -1){
           if (C->sym[idx_arr[3*iA+2]] == NS ||
               idx_A[i+1] != idx_C[idx_arr[3*iA+2]+1]){
-            if (new_contraction != NULL)
-              nA->sym[i] = NS;
+            if (new_contraction != NULL){
+              nA_sym[i] = NS;
+              nA->set_sym(nA_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i;
           }
         } else {
           if (idx_arr[3*idx_A[i+1]+2] != -1){
-            if (new_contraction != NULL)
-              nA->sym[i] = NS;
+            if (new_contraction != NULL){
+              nA_sym[i] = NS;
+              nA->set_sym(nA_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i;
           }       
@@ -546,21 +557,28 @@ namespace CTF_int {
     }
 
    
+    int nB_sym[B->order];
+    if (new_contraction != NULL)
+      memcpy(nB_sym, nB->sym, sizeof(int)*nB->order);
     for (i=0; i<B->order; i++){
       if (B->sym[i] != NS){
         iB = idx_B[i];
         if (idx_arr[3*iB+0] != -1){
           if (A->sym[idx_arr[3*iB+0]] == NS ||
               idx_B[i+1] != idx_A[idx_arr[3*iB+0]+1]){
-            if (new_contraction != NULL)
-              nB->sym[i] = NS;
+            if (new_contraction != NULL){
+              nB_sym[i] = NS;
+              nB->set_sym(nB_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i+1;
           }
         } else {
           if (idx_arr[3*idx_B[i+1]+0] != -1){
-            if (new_contraction != NULL)
-              nB->sym[i] = NS;
+            if (new_contraction != NULL){
+              nB_sym[i] = NS;
+              nB->set_sym(nB_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i+1;
           }       
@@ -568,15 +586,19 @@ namespace CTF_int {
         if (idx_arr[3*iB+2] != -1){
           if (C->sym[idx_arr[3*iB+2]] == NS || 
               idx_B[i+1] != idx_C[idx_arr[3*iB+2]+1]){
-            if (new_contraction != NULL)
-              nB->sym[i] = NS;
+            if (new_contraction != NULL){
+              nB_sym[i] = NS;
+              nB->set_sym(nB_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i+1;
           }
         } else {
           if (idx_arr[3*idx_B[i+1]+2] != -1){
-            if (new_contraction != NULL)
-              nB->sym[i] = NS;
+            if (new_contraction != NULL){
+              nB_sym[i] = NS;
+              nB->set_sym(nB_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i+1;
           }       
@@ -603,35 +625,46 @@ namespace CTF_int {
       }
     }
     if (!is_preserv){
+      int nC_sym[C->order];
+      if (new_contraction != NULL)
+        memcpy(nC_sym, nC->sym, sizeof(int)*nC->order);
       for (i=0; i<C->order; i++){
         if (C->sym[i] != NS){
           iC = idx_C[i];
           if (idx_arr[3*iC+1] != -1){
             if (B->sym[idx_arr[3*iC+1]] == NS ||
                 idx_C[i+1] != idx_B[idx_arr[3*iC+1]+1]){
-              if (new_contraction != NULL)
-                nC->sym[i] = NS;
+              if (new_contraction != NULL){
+                nC_sym[i] = NS;
+                nC->set_sym(nC_sym);
+              }
               CTF_int::cdealloc(idx_arr); 
               return 3*i+2;
             }
           } else if (idx_arr[3*idx_C[i+1]+1] != -1){
-            if (new_contraction != NULL)
-              nC->sym[i] = NS;
+            if (new_contraction != NULL){
+              nC_sym[i] = NS;
+              nC->set_sym(nC_sym);
+            }
             CTF_int::cdealloc(idx_arr); 
             return 3*i+2;
           }       
           if (idx_arr[3*iC+0] != -1){
             if (A->sym[idx_arr[3*iC+0]] == NS ||
                 idx_C[i+1] != idx_A[idx_arr[3*iC+0]+1]){
-              if (new_contraction != NULL)
-                nC->sym[i] = NS;
+              if (new_contraction != NULL){
+                nC_sym[i] = NS;
+                nC->set_sym(nC_sym);
+              }
               CTF_int::cdealloc(idx_arr); 
               return 3*i+2;
             }
           } else if (idx_arr[3*iC+0] == -1){
             if (idx_arr[3*idx_C[i+1]] != -1){
-              if (new_contraction != NULL)
-                nC->sym[i] = NS;
+              if (new_contraction != NULL){
+                nC_sym[i] = NS;
+                nC->set_sym(nC_sym);
+              }
               CTF_int::cdealloc(idx_arr); 
               return 3*i+2;
             }       
@@ -645,8 +678,10 @@ namespace CTF_int {
         iA2 = idx_A[i+1];
         if (idx_arr[3*iA+2] == -1 &&
             idx_arr[3*iA2+2] == -1){
-          if (new_contraction != NULL)
-            nA->sym[i] = NS;
+          if (new_contraction != NULL){
+            nA_sym[i] = NS;
+            nA->set_sym(nA_sym);
+          }
           CTF_int::cdealloc(idx_arr); 
           return 3*i;
         }
@@ -658,8 +693,10 @@ namespace CTF_int {
         iB2 = idx_B[i+1];
         if (idx_arr[3*iB+2] == -1 &&
             idx_arr[3*iB2+2] == -1){
-          if (new_contraction != NULL)
-            nB->sym[i] = NS;
+          if (new_contraction != NULL){
+            nB_sym[i] = NS;
+            nB->set_sym(nB_sym);
+          }
           CTF_int::cdealloc(idx_arr); 
           return 3*i+1;
         }
@@ -1487,8 +1524,10 @@ namespace CTF_int {
     tC->topo = topo;
     
     /* Map the weigh indices of A, B, and C*/
-    ret = map_weigh_indices(idx_arr, idx_weigh, num_tot, num_weigh, topo, tA, tB, tC);
+
+
     int stat;
+    ret = map_weigh_indices(idx_arr, idx_weigh, num_tot, num_weigh, topo, tA, tB, tC);
     do {
       if (ret == NEGATIVE) {
         stat = ret;

src/interface/tensor.cxx

@@ -10,16 +10,6 @@ namespace CTF {
   template<typename dtype, bool is_ord>
   Tensor<dtype, is_ord>::Tensor() : CTF_int::tensor() { }
 
-  template<typename dtype, bool is_ord>
-  Tensor<dtype, is_ord>::Tensor(tensor const & A,
-                                bool           copy)
-    : CTF_int::tensor(&A, copy) { }
-
-  template<typename dtype, bool is_ord>
-  Tensor<dtype, is_ord>::Tensor(tensor const & A,
-                                World &        world_)
-    : CTF_int::tensor(A.sr, A.order, A.lens, A.sym, A.wrld, 1, A.name, A.profile) { }
-
 
   template<typename dtype, bool is_ord>
   Tensor<dtype, is_ord>::Tensor(int                       order,
@@ -79,6 +69,21 @@ namespace CTF {
 #endif
   }
 
+  template<typename dtype, bool is_ord>
+  Tensor<dtype, is_ord>::Tensor(tensor const & A,
+                                bool           copy)
+    : CTF_int::tensor(&A, copy) { }
+
+  template<typename dtype, bool is_ord>
+  Tensor<dtype, is_ord>::Tensor(tensor const & A,
+                                World &        world_)
+    : CTF_int::tensor(A.sr, A.order, A.lens, A.sym, world_, 1, A.name, A.profile) { }
+    //: CTF_int::tensor(A.sr, A.order, A.lens, A.sym, A.wrld, 1, A.name, A.profile) { }
+
+  template<typename dtype, bool is_ord>
+  Tensor<dtype, is_ord>::Tensor(tensor &    A,
+                                int const * new_sym)
+    : CTF_int::tensor(&A, new_sym){ }
 
 
 

src/interface/tensor.h

@@ -136,22 +136,6 @@ namespace CTF {
        */
       Tensor();
 
-      /**
-       * \brief copies a tensor (setting data to zero or copying A)
-       * \param[in] A tensor to copy
-       * \param[in] copy whether to copy the data of A into the new tensor
-       */
-      Tensor(tensor const & A,
-             bool           copy = true);
-      
-      /**
-       * \brief creates a zeroed out copy (data not copied) of a tensor in a different world
-       * \param[in] A tensor whose characteristics to copy
-       * \param[in] world_ a world for the tensor we are creating to live in, can be different from A
-       */
-      Tensor(tensor const & A,
-             World &        wrld);
-
       /**
        * \brief defines tensor filled with zeros on the default algstrct
        * \param[in] order_ number of dimensions of tensor
@@ -204,6 +188,35 @@ namespace CTF {
              char const *              name=NULL,
              int                       profile=0);
 
+
+      /**
+       * \brief copies a tensor (setting data to zero or copying A)
+       * \param[in] A tensor to copy
+       * \param[in] copy whether to copy the data of A into the new tensor
+       */
+      Tensor(tensor const & A,
+             bool           copy = true);
+
+      /**
+       * \brief repacks the tensor other to a different symmetry 
+       *        (assumes existing data contains the symmetry and keeps only values with indices in increasing order)
+       * WARN: LIMITATION: new_sym must cause unidirectional structural changes, i.e. {NS,NS}->{SY,NS} OK, {SY,NS}->{NS,NS} OK, {NS,NS,SY,NS}->{SY,NS,NS,NS} NOT OK!
+       * \param[in] A tensor to copy
+       * \param[in] new_sym new symmetry array (replaces this->sym)
+       */
+      Tensor(tensor &    A,
+             int const * new_sym);
+      
+      /**
+       * \brief creates a zeroed out copy (data not copied) of a tensor in a different world
+       * \param[in] A tensor whose characteristics to copy
+       * \param[in] world_ a world for the tensor we are creating to live in, can be different from A
+       */
+      Tensor(tensor const & A,
+             World &        wrld);
+
+
+
       /**
        * \brief defines tensor filled with zeros on the default algstrct on a user-specified distributed layout
        * \param[in] order_ number of dimensions of tensor

src/summation/summation.cxx

@@ -1307,10 +1307,16 @@ namespace CTF_int {
     if (nnew_sum != NULL && sidx != -1){
       if(sidx%2 == 0){
         new_sum->A = new tensor(A, 0, 0);
-        new_sum->A->sym[sidx/2] = NS;
+        int nA_sym[A->order];
+        memcpy(nA_sym, new_sum->A->sym, sizeof(int)*new_sum->A->order);
+        nA_sym[sidx/2] = NS;
+        new_sum->A->set_sym(nA_sym);
       } else {
         new_sum->B = new tensor(B, 0, 0);
-        new_sum->B->sym[sidx/2] = NS;
+        int nB_sym[B->order];
+        memcpy(nB_sym, new_sum->B->sym, sizeof(int)*new_sum->B->order);
+        nB_sym[sidx/2] = NS;
+        new_sum->B->set_sym(nB_sym);
       }
     }
     CTF_int::cdealloc(idx_arr);

src/symmetry/symmetrization.cxx

@@ -285,7 +285,7 @@ namespace CTF_int {
         idx_map_B[i] = i;
       }
       
-      if (scal_diag){
+      if (0){//scal_diag){
         //FIXME: this is not robust when doing e.g. {SY, SY, SY, NS} -> {SY, NS, SY, NS}
         for (i=-num_sy_neg-1; i<num_sy; i++){
           if (i==-1) continue;

src/tensor/untyped_tensor.cxx

@@ -91,6 +91,56 @@ namespace CTF_int {
     } else data = NULL;
 
   }
+ 
+  tensor::tensor(tensor * other, int const * new_sym){
+    char * nname = (char*)alloc(strlen(other->name) + 2);
+    char d[] = "\'";
+    strcpy(nname, other->name);
+    strcat(nname, d);
+    if (other->wrld->rank == 0) {
+      DPRINTF(1,"Repacking tensor %s into %s\n",other->name);
+    }
+
+    bool has_chng=false, less_sym=false, more_sym=false;
+    for (int i=0; i<other->order; i++){
+      if (other->sym[i] != new_sym[i]){
+        has_chng = true;
+        if (other->sym[i] == NS){
+          assert(!less_sym); 
+          more_sym = true;
+        }
+        if (new_sym[i] == NS){
+          assert(!more_sym); 
+          less_sym = true;
+        }
+      }
+    }
+  
+    this->has_zero_edge_len = other->has_zero_edge_len;
+
+  
+    if (!less_sym && !more_sym){
+      this->init(other->sr, other->order, other->lens,
+                 new_sym, other->wrld, 0, nname,
+                 other->profile);
+      copy_tensor_data(other);
+      if (has_chng)
+        zero_out_padding();
+    } else {
+      this->init(other->sr, other->order, other->lens,
+                 new_sym, other->wrld, 1, nname,
+                 other->profile);
+      int idx[order];
+      for (int j=0; j<order; j++){
+        idx[j] = j;
+      }
+      summation ts(other, idx, sr->mulid(), this, idx, sr->addid());
+      ts.sum_tensors(true);
+    }
+    cdealloc(nname);
+  }
+
+
 
   void tensor::copy_tensor_data(tensor const * other){
     //FIXME: do not unfold
@@ -163,7 +213,6 @@ namespace CTF_int {
     memcpy(this->pad_edge_len, other->pad_edge_len, sizeof(int)*other->order);
     memcpy(this->padding, other->padding, sizeof(int)*other->order);
     memcpy(this->sym, other->sym, sizeof(int)*other->order);
-    memcpy(this->sym_table, other->sym_table, sizeof(int)*other->order*other->order);
     this->is_mapped = other->is_mapped;
     this->is_cyclic = other->is_cyclic;
     this->topo      = other->topo;
@@ -223,13 +272,8 @@ namespace CTF_int {
     this->pad_edge_len = (int*)CTF_int::alloc(order*sizeof(int));
     memcpy(this->pad_edge_len, lens, order*sizeof(int));
     this->sym      = (int*)CTF_int::alloc(order*sizeof(int));
-    if (sym_ == NULL)
-      std::fill(this->sym, this->sym+order, NS);
-    else
-      memcpy(this->sym, sym_, order*sizeof(int));
-  
-    this->sym_table = (int*)CTF_int::alloc(order*order*sizeof(int));
-    memset(this->sym_table, 0, order*order*sizeof(int));
+    sym_table = (int*)CTF_int::alloc(order*order*sizeof(int));
+    this->set_sym (sym_);
     this->edge_map  = new mapping[order];
 
     /* initialize map array and symmetry table */
@@ -238,17 +282,17 @@ namespace CTF_int {
       this->edge_map[i].type       = NOT_MAPPED;
       this->edge_map[i].has_child  = 0;
       this->edge_map[i].np         = 1;
-      if (this->sym[i] != NS) {
+      /*if (this->sym[i] != NS) {
         //FIXME: keep track of capabilities of algberaic structure and add more robust property checking
-/*        if (this->sym[i] == AS && !sr->is_ring){ 
+        if (this->sym[i] == AS && !sr->is_ring){ 
           if (wrld->rank == 0){
             printf("CTF ERROR: It is illegal to define antisymmetric tensor must be defined on a ring, yet no additive inverse was provided for this algstrct (see algstrct constructor), aborting.\n");
           }
           ABORT;
-        }*/
+        }
         this->sym_table[(i+1)+i*order] = 1;
         this->sym_table[(i+1)*order+i] = 1;
-      }
+      }*/
     }
     /* Set tensor data to zero. */
     if (alloc_data){
@@ -451,9 +495,11 @@ namespace CTF_int {
   void tensor::print_map(FILE * stream, bool allcall) const {
     if (!allcall || wrld->rank == 0){
       printf("printing mapping of %s\n",name);
-      printf("mapped to order %d topology with dims:",topo->order);
-      for (int dim=0; dim<topo->order; dim++){
-        printf(" %d ",topo->lens[dim]);
+      if (topo != NULL){
+        printf("mapped to order %d topology with dims:",topo->order);
+        for (int dim=0; dim<topo->order; dim++){
+          printf(" %d ",topo->lens[dim]);
+        }
       }
       printf("\nCTF: sym  len  tphs  pphs  vphs\n");
       for (int dim=0; dim<order; dim++){
@@ -1783,46 +1829,20 @@ namespace CTF_int {
     return SUCCESS;
 
   }
-  
-  void tensor::repack(int const * new_sym){
-    bool has_chng=false, less_sym=false, more_sym=false;
+ 
+  void tensor::set_sym(int const * sym_){
+    if (sym_ == NULL)
+      std::fill(this->sym, this->sym+order, NS);
+    else
+      memcpy(this->sym, sym_, order*sizeof(int));
+ 
+    memset(sym_table, 0, order*order*sizeof(int));
     for (int i=0; i<order; i++){
-      if (sym[i] != new_sym[i]){
-        has_chng = true;
-        if (sym[i] == NS){
-          assert(!less_sym); 
-          more_sym = true;
-        }
-        if (new_sym[i] == NS){
-          assert(!more_sym); 
-          less_sym = true;
-        }
-      }
-    }  
-    if (!has_chng) return;
-    if (less_sym){
-      tensor tcpy(this, true, true);
-      memcpy(sym, new_sym, order*sizeof(int));
-      int idx[order];
-      for (int j=0; j<order; j++){
-        idx[j] = j;
-      }
-      summation ts(&tcpy, idx, sr->mulid(), this, idx, sr->addid());
-      ts.execute();
-    } else if (more_sym){
-      tensor tcpy(this, true, true);
-      memcpy(sym, new_sym, order*sizeof(int));
-      int idx[order];
-      for (int j=0; j<order; j++){
-        idx[j] = j;
+      if (this->sym[i] != NS) {
+        sym_table[(i+1)+i*order] = 1;
+        sym_table[(i+1)*order+i] = 1;
       }
-      summation ts(&tcpy, idx, sr->mulid(), this, idx, sr->addid());
-      ts.sum_tensors(true);
-    } else {
-      memcpy(sym, new_sym, order*sizeof(int));
-      zero_out_padding();
     }
   }
-
 }
 

src/tensor/untyped_tensor.h

@@ -135,15 +135,6 @@ namespace CTF_int {
 
       /** \brief destructor */
       void free_self();
-
-      /**
-       * \brief creates tensor copy, unfolds other if other is folded
-       * \param[in] other tensor to copy
-       * \param[in] copy whether to copy mapping and data
-       * \param[in] alloc_data whether th allocate data
-       */
-      tensor(tensor const * other, bool copy = 1, bool alloc_data = 1);
-
       /**
        * \brief defines a tensor object with some mapping (if alloc_data)
        * \param[in] sr defines the tensor arithmetic for this tensor
@@ -163,6 +154,23 @@ namespace CTF_int {
              char const *     name=NULL,
              bool             profile=1);
 
+      /**
+       * \brief creates tensor copy, unfolds other if other is folded
+       * \param[in] other tensor to copy
+       * \param[in] copy whether to copy mapping and data
+       * \param[in] alloc_data whether th allocate data
+       */
+      tensor(tensor const * other, bool copy = 1, bool alloc_data = 1);
+
+      /**
+       * \brief repacks the tensor other to a different symmetry 
+       *        (assumes existing data contains the symmetry and keeps only values with indices in increasing order)
+       * WARN: LIMITATION: new_sym must cause unidirectional structural changes, i.e. {NS,NS}->{SY,NS} OK, {SY,NS}->{NS,NS} OK, {NS,NS,SY,NS}->{SY,NS,NS,NS} NOT OK!
+       * \param[in] other tensor to copy
+       * \param[in] new_sym new symmetry array (replaces this->sym)
+       */
+      tensor(tensor * other, int const * new_sym);
+
       /**
        * \brief compute the cyclic phase of each tensor dimension
        * \return int * of cyclic phases
@@ -538,14 +546,10 @@ namespace CTF_int {
                        tensor *&   new_tsr,
                        int **      idx_map_new);
 
-      /**
-       * \brief repacks the tensor to a different symmetry 
-       *        (assumes existing data contains the symmetry and keeps only values with indices in increasing order)
-       * WARN: LIMITATION: new_sym must cause unidirectional structural changes, i.e. {NS,NS}->{SY,NS} OK, {SY,NS}->{NS,NS} OK, {NS,NS,SY,NS}->{SY,NS,NS,NS} NOT OK!
-       * \param[in] new_sym new symmetry array (replaces this->sym)
-       */
-      void repack(int const * new_sym);
-
+      /** \brief sets symmetry, WARNING: for internal use only !!!!
+        * \param[in] sym
+        */
+      void set_sym(int const * sym); 
   };
 }
 

test/Makefile

 include ../config.mk
 
-TESTS = test_suite diag_sym sym3 readwrite_test \
+TESTS = test_suite diag_sym sym3 readwrite_test sy_times_ns repack \
         ccsdt_t3_to_t2 ccsdt_map_test multi_tsr_sym diag_ctr readall_test
 
 SCALA_TESTS = pgemm_test nonsq_pgemm_test 

test/repack.cxx

+/*Copyright (c) 2011, Edgar Solomonik, all rights reserved.*/
+
+/** \addtogroup examples 
+  * @{ 
+  * \defgroup repack repack 
+  * @{ 
+  * \brief Tests contraction of a symmetric index group with a nonsymmetric one
+  */
+
+#include <ctf.hpp>
+
+using namespace CTF;
+
+int repack(int     n,
+           World & dw){
+  int rank, i, num_pes, pass;
+  int64_t np;
+  double * pairs;
+  int64_t * indices;
+  
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &num_pes);
+
+
+  int shapeN4[] = {NS,NS,NS,NS};
+  int shapeS4[] = {NS,NS,SY,NS};
+  int sizeN4[] = {n,n,n,n};
+
+  //* Creates distributed tensors initialized with zeros
+  Tensor<> An(4, sizeN4, shapeN4, dw);
+  Tensor<> As(4, sizeN4, shapeS4, dw);
+
+  As.read_local(&np, &indices, &pairs);
+  for (i=0; i<np; i++ ) pairs[i] = drand48()-.5; //(1.E-3)*sin(indices[i]);
+  As.write(np, indices, pairs);
+  An.write(np, indices, pairs);
+
+  Tensor<> Anr(An, shapeS4);
+ 
+  Anr["ijkl"] -= As["ijkl"];
+
+  double norm = Anr.norm2();
+
+  if (norm < 1.E-6)
+    pass = 1;
+  else
+    pass = 0;
+  
+  if (!pass)
+    printf("{NS -> SY repack} failed \n");
+  else {
+    Tensor<> Anur(As, shapeN4);
+    Tensor<> Asur(As, shapeN4);
+    Asur["ijkl"] = 0.0;
+    Asur.write(np, indices, pairs);
+    Anur["ijkl"] -= Asur["ijkl"];
+
+    norm = Anur.norm2();
+
+    if (norm < 1.E-6){
+      pass = 1;
+      if (rank == 0)
+        printf("{NS -> SY -> NS repack} passed \n");
+    } else {
+      pass = 0;
+      if (rank == 0)
+        printf("{SY -> NS repack} failed \n");
+    }
+
+  }
+  free(pairs);
+  free(indices);
+  return pass;
+} 
+
+
+#ifndef TEST_SUITE
+char* getCmdOption(char ** begin,
+                   char ** end,
+                   const   std::string & option){
+  char ** itr = std::find(begin, end, option);
+  if (itr != end && ++itr != end){
+    return *itr;
+  }
+  return 0;
+}
+
+
+int main(int argc, char ** argv){
+  int rank, np, n;
+  int in_num = argc;
+  char ** input_str = argv;
+
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &np);
+
+  if (getCmdOption(input_str, input_str+in_num, "-n")){
+    n = atoi(getCmdOption(input_str, input_str+in_num, "-n"));
+    if (n < 0) n = 7;
+  } else n = 7;
+
+
+  {
+    World dw(argc, argv);
+    repack(n, dw);
+  }
+
+  MPI_Finalize();
+  return 0;
+}
+/**
+ * @} 
+ * @}
+ */
+
+#endif

test/sy_times_ns.cxx

+/*Copyright (c) 2011, Edgar Solomonik, all rights reserved.*/
+
+/** \addtogroup examples 
+  * @{ 
+  * \defgroup sy_times_ns sy_times_ns 
+  * @{ 
+  * \brief Tests contraction of a symmetric index group with a nonsymmetric one
+  */
+
+#include <ctf.hpp>
+
+using namespace CTF;
+
+int sy_times_ns(int     n,
+                World & dw){
+  int rank, i, num_pes, pass;
+  int64_t np;
+  double * pairs;
+  int64_t * indices;
+  
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &num_pes);
+
+
+  int shapeN4[] = {NS,NS,NS,NS};
+  int sizeN4[] = {n,n,n,n};
+
+  //* Creates distributed tensors initialized with zeros
+  Tensor<> B(4, sizeN4, shapeN4, dw);
+
+  Matrix<> A(n, n, SY, dw);
+  Matrix<> An(n, n, NS, dw);
+  Matrix<> C(n, n, SY, dw, "C");
+  Matrix<> Cn(n, n, NS, dw, "Cn");
+
+  srand48(13*rank);
+
+
+  A.read_local(&np, &indices, &pairs);
+  for (i=0; i<np; i++ ) pairs[i] = drand48()-.5; //(1.E-3)*sin(indices[i]);
+  A.write(np, indices, pairs);
+  free(pairs);
+  free(indices);
+  B.read_local(&np, &indices, &pairs);
+  for (i=0; i<np; i++ ) pairs[i] = drand48()-.5; //(1.E-3)*sin(indices[i]);
+  B.write(np, indices, pairs);
+  free(pairs);
+  free(indices);
+
+  An["ij"] = A["ij"];
+
+  C["ij"] = A["ij"]*B["ijkl"];
+  Cn["ij"] = An["ij"]*B["ijkl"];
+  Cn["ji"] += An["ij"]*B["ijkl"];
+
+
+  Cn["ij"] -= C["ij"];
+
+  double norm = Cn.norm2();
+  
+  if (norm < 1.E-10){
+    pass = 1;
+    if (rank == 0)
+      printf("{C[\(ij)\"]=A[\"(ij)\"]*B[\"ijkl\"]} passed \n");
+  } else {
+    pass = 0;
+    if (rank == 0)
+      printf("{C[\(ij)\"]=A[\"(ij)\"]*B[\"ijkl\"]} failed \n");
+  }
+  return pass;
+} 
+
+
+#ifndef TEST_SUITE
+char* getCmdOption(char ** begin,
+                   char ** end,
+                   const   std::string & option){
+  char ** itr = std::find(begin, end, option);
+  if (itr != end && ++itr != end){
+    return *itr;
+  }
+  return 0;
+}
+
+
+int main(int argc, char ** argv){
+  int rank, np, n;
+  int in_num = argc;
+  char ** input_str = argv;
+
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &np);
+
+  if (getCmdOption(input_str, input_str+in_num, "-n")){
+    n = atoi(getCmdOption(input_str, input_str+in_num, "-n"));
+    if (n < 0) n = 7;
+  } else n = 7;
+
+
+  {
+    World dw(argc, argv);
+    sy_times_ns(n, dw);
+  }
+
+  MPI_Finalize();
+  return 0;
+}
+/**
+ * @} 
+ * @}
+ */
+
+#endif

test/test_suite.cxx

@@ -22,6 +22,8 @@
 #include "readall_test.cxx"
 #include "../examples/subworld_gemm.cxx"
 #include "multi_tsr_sym.cxx"
+#include "repack.cxx"
+#include "sy_times_ns.cxx"
 
 using namespace CTF;
 
@@ -152,9 +154,17 @@ int main(int argc, char ** argv){
       printf("Testing diagonal write with n = %d:\n",n);
     pass.push_back(readwrite_test(n, dw));
     
-  if (rank == 0)
+    if (rank == 0)
       printf("Testing readall test with n = %d m = %d:\n",n,n*n);
     pass.push_back(readall_test(n, n*n, dw));
+    
+    if (rank == 0)
+      printf("Testing repack with n = %d:\n",n);
+    pass.push_back(repack(n,dw));
+    
+    if (rank == 0)
+      printf("Testing SY times NS with n = %d:\n",n);
+    pass.push_back(sy_times_ns(n,dw));
 #if 0
     if (rank == 0)
       printf("Testing skew-symmetric Strassen's algorithm with n = %d:\n",n*n);