diff --git a/dimension_t/test/is_good.c b/dimension_t/test/is_good.c
index 19adc8c..6bd7f5e 100644
--- a/dimension_t/test/is_good.c
+++ b/dimension_t/test/is_good.c
@@ -51,6 +51,23 @@ TEST(SubDim){
 
 }
 
+TEST(SubDim){
+  dimension *D=create_dim(4);
+  D->perm[0]=2;
+  D->perm[1]=3;
+  D->perm[2]=5;
+  D->perm[3]=6;
+
+
+  dimension *d_head2 = sub_dim_head(D,2);
+  
+  
+  EXPECT_EQ(d_head2->rank, 2*3);
+
+  dimension *d_tail2 = sub_dim_tail(D,2);
+  EXPECT_EQ(d_tail2->rank, 5*6);
+
+}
 
 TEST(Coord_linear){
   dimension *D=create_dim(4);
diff --git a/tensor_t/src/tensor_t/tensor_t.c b/tensor_t/src/tensor_t/tensor_t.c
index c25f8e6..b5b8cb9 100644
--- a/tensor_t/src/tensor_t/tensor_t.c
+++ b/tensor_t/src/tensor_t/tensor_t.c
@@ -33,6 +33,53 @@ void printArraySzt(size_t *a, size_t sz,char *msg){
     return r_tens;\
   }\
 \
+  tensor_##type * sub_tensor_head_##type(tensor_##type *rootens, size_t minuSubdim, size_t rankInDim){\
+    dimension *rdim= rootens->dim;\
+    dimension *dS_t = sub_dim_tail(rdim,rdim->size - minuSubdim);\
+    if(rankInDim < dS_t->rank){\
+      dimension *dS_h = sub_dim_head(rdim,minuSubdim);\
+      tensor_##type *ret_ens = malloc(sizeof(tensor_##type));\
+      ret_ens->dim = dS_h;\
+      /*ret_ens->x = (rootens->x)+rankInDim*dS_t->rank;*/\
+      if(endian){\
+        ret_ens->x = malloc(sizeof(type)*dS_h->rank);\
+        for(size_t i=0; i<dS_h->rank; ++i){\
+          ret_ens->x[i]=rootens->x[i*dS_t->rank + rankInDim];\
+          /*printf("%ld: [i:%ld] | %ld : [%ld ]\n",dS_t->rank, i,dS_h->rank,i*dS_h->rank + rankInDim);*/\
+          \
+        }\
+      }else{\
+        ret_ens->x = (rootens->x)+rankInDim*dS_h->rank;\
+      \
+      }\
+      return ret_ens;\
+    }\
+    return NULL;\
+  }\
+  \
+  tensor_##type * sub_tensor_tail_##type(tensor_##type *rootens, size_t minuSubdim, size_t rankInDim){\
+    dimension *rdim= rootens->dim;\
+    dimension *dS_h = sub_dim_head(rdim,rdim->size - minuSubdim);\
+    if(rankInDim < dS_h->rank){\
+      dimension *dS_t = sub_dim_tail(rdim,minuSubdim);\
+      tensor_##type *ret_ens = malloc(sizeof(tensor_##type));\
+      ret_ens->dim = dS_t;\
+      if(endian==false){\
+        ret_ens->x = malloc(sizeof(type)*dS_t->rank);\
+        for(size_t i=0; i<dS_t->rank; ++i){\
+          ret_ens->x[i]=rootens->x[i*dS_h->rank + rankInDim];\
+          /*printf("%ld: [i:%ld] | %ld : [%ld ]\n",dS_t->rank, i,dS_h->rank,i*dS_h->rank + rankInDim);*/\
+          \
+        }\
+      }else{\
+        ret_ens->x = (rootens->x)+rankInDim*dS_t->rank;\
+      \
+      }\
+      return ret_ens;\
+    }\
+    return NULL;\
+  }\
+  \
 void tensorProdNotOpt_##type(tensor_##type **MM, tensor_##type *M0, tensor_##type *M1) {  \
   dimension *dd;  \
     add_dimension(&dd, M0->dim, M1->dim); \
@@ -65,7 +112,10 @@ void tensorProd_##type(tensor_##type **MM, tensor_##type *M0, tensor_##type *M1)
     size_t m_idx;\
     for(size_t i=0; i<M0->dim->rank; ++i){\
       for(size_t j=0; j<M1->dim->rank; ++j){\
-          m_idx= i*M1->dim->rank + j ;\
+          if(endian)\
+            m_idx= i*M1->dim->rank + j ;\
+          else\
+            m_idx= i+M0->dim->rank * j ;\
           M->x[m_idx]=M0->x[i]*M1->x[j];\
           /*printf("[%ld|%ld:(%ld,%ld)]",x_idx++,m_idx,i,j);*/\
       }\
@@ -133,8 +183,14 @@ void tensorContractnProd_##type(tensor_##type** MM, tensor_##type *M0, tensor_##
 \
     size_t a0_id, a1_id, n0_id, n1_id;\
     for (size_t i = 0; i < M->dim->rank; i++) {\
-        a0_id=i/dSub1->rank;\
-        a1_id=i%dSub1->rank;\
+        if(endian){\
+          a0_id=i/dSub1->rank;\
+          a1_id=i%dSub1->rank;\
+        }\
+        else{\
+          a0_id=i%dSub0->rank;\
+          a1_id=i/dSub0->rank;\
+        }\
         /*vCoordFromLin(coord, i, M->dim);\
         subArray(coord0, coord, 0, len0, 0);\
         subArray(coord1, coord, 0, len1, len0);\
@@ -142,8 +198,14 @@ void tensorContractnProd_##type(tensor_##type** MM, tensor_##type *M0, tensor_##
         printf("i:%ld=>  c1: %ld vs %ld \n",i,LineFromCoord(coord1,dSub1),a1_id);*/\
         M->x[i] = 0;\
         for (size_t k = 0; k < dM->rank; k++) {\
+          if(endian){\
             n0_id= a0_id*dM->rank + k;\
             n1_id= a1_id + dSub1->rank * k;\
+          }\
+          else{\
+            n0_id= a0_id + dSub0->rank * k;\
+            n1_id= a1_id*dM->rank + k;\
+          }\
             M->x[i] += M0->x[n0_id] * M1->x[n1_id];\
             /*vCoordFromLin(Koord, k, dM);\
             concatArray(coordM0, coord0, Koord, 0, 0, len0, 0, contractionNumber);\
diff --git a/tensor_t/src/tensor_t/tensor_t.h b/tensor_t/src/tensor_t/tensor_t.h
index f55c2e7..5a3dabf 100644
--- a/tensor_t/src/tensor_t/tensor_t.h
+++ b/tensor_t/src/tensor_t/tensor_t.h
@@ -12,6 +12,8 @@ struct tensor_##type{\
 };\
 typedef struct tensor_##type tensor_##type;\
 tensor_##type * CREATE_TENSOR_##type(dimension *dim); \
+tensor_##type * sub_tensor_head_##type(tensor_##type *rootens, size_t minuSubdim, size_t rankInDim); \
+tensor_##type * sub_tensor_tail_##type(tensor_##type *rootens, size_t minuSubdim, size_t rankInDim); \
 void tensorProdNotOpt_##type(tensor_##type **MM, tensor_##type *M0, tensor_##type *M1); \
 void tensorProd_##type(tensor_##type **MM, tensor_##type *M0, tensor_##type *M1); \
 void tensorContractnProd_##type(tensor_##type **MM, tensor_##type *M0, tensor_##type *M1, size_t contractionNumber); \
diff --git a/tensor_t/test/is_good.c b/tensor_t/test/is_good.c
index b560424..5005507 100644
--- a/tensor_t/test/is_good.c
+++ b/tensor_t/test/is_good.c
@@ -31,6 +31,10 @@ TEST(rank){
 }
 
 void print_tensor_float(tensor_TYPE_FLOAT *M, char *msg){
+  if(M==NULL) {
+    LOG("empty tensor | %s ===============\n",msg);
+    return;
+  }
   LOG("================= %s ===============\n",msg);
   for(size_t i=0; i<M->dim->rank;++i)
       LOG("[%ld]: %f ",i,M->x[i]);
@@ -47,6 +51,100 @@ void print_tensor_double(tensor_TYPE_DOUBLE *M, char *msg){
     LOG("%s","\n");
 }
 
+TEST(tensorSubhead_ ){
+  dimension *d0=create_dim(3);
+
+  d0->perm[0]=4;
+  d0->perm[1]=3;
+  d0->perm[2]=5;
+  
+
+  updateRankDim(d0);
+
+
+  tensor_TYPE_FLOAT *M0 = CREATE_TENSOR_TYPE_FLOAT(d0);
+
+  LOG("M0->dim->rank = %ld\n",M0->dim->rank);
+  for(size_t i=0; i<M0->dim->rank;++i) M0->x[i]=i*0.1 +1;
+
+  endian=false;
+  size_t rnkId = 1;
+  tensor_TYPE_FLOAT *s2h = sub_tensor_head_TYPE_FLOAT(M0,1,rnkId);
+
+  print_tensor_float(M0, "M0 for sub");
+  print_tensor_float(s2h, "sub s2h of M0");
+
+  printDebug_dimension(M0->dim," M0 dimension  ");
+  printDebug_dimension(s2h->dim," s2h dimension  ");
+ //EXPECT_ARRAY_EQ_TYPE_FLOAT(M0->x,s2h->dim->rank,s2h->x,s2h->dim->rank);
+ 
+  for(size_t i1=0; i1<d0->perm[1];++i1){
+    for(size_t i0=0; i0<d0->perm[0];++i0){
+      printf("EXPECT_EQ_TYPE_FLOAT c0[%ld](%ld,%ld,%ld) : s[%ld](%ld,%ld) \n",LineFromCoord((size_t[]){i0,i1,rnkId},M0->dim),rnkId,i0,i1,LineFromCoord((size_t[]){i0,i1},s2h->dim),i0,i1);
+      EXPECT_EQ_TYPE_FLOAT(M0->x[LineFromCoord((size_t[]){i0,i1,rnkId},M0->dim)],s2h->x[LineFromCoord((size_t[]){i0,i1},s2h->dim)]);
+    }  
+  }
+}
+
+
+TEST(tensorSubtail ){
+  dimension *d0=create_dim(3);
+
+  d0->perm[0]=4;
+  d0->perm[1]=3;
+  d0->perm[2]=5;
+  
+
+  updateRankDim(d0);
+
+
+  tensor_TYPE_FLOAT *M0 = CREATE_TENSOR_TYPE_FLOAT(d0);
+
+  LOG("M0->dim->rank = %ld\n",M0->dim->rank);
+  for(size_t i=0; i<M0->dim->rank;++i) M0->x[i]=i*0.1 +1;
+
+
+  tensor_TYPE_FLOAT *s1t = sub_tensor_tail_TYPE_FLOAT(M0,1,5);
+
+  print_tensor_float(M0, "M0 for sub");
+  print_tensor_float(s1t, "sub s2t of M0 from 5");
+
+}
+
+
+TEST(tensorSubtail ){
+  dimension *d0=create_dim(3);
+
+  d0->perm[0]=4;
+  d0->perm[1]=3;
+  d0->perm[2]=5;
+  
+
+  updateRankDim(d0);
+
+
+  tensor_TYPE_FLOAT *M0 = CREATE_TENSOR_TYPE_FLOAT(d0);
+
+  LOG("M0->dim->rank = %ld\n",M0->dim->rank);
+  for(size_t i=0; i<M0->dim->rank;++i) M0->x[i]=i*0.1 +1;
+
+  //endian=false;
+  size_t rnkId=3;
+  tensor_TYPE_FLOAT *s2t = sub_tensor_tail_TYPE_FLOAT(M0,1,rnkId);
+
+
+
+  print_tensor_float(M0, "M0 for sub");
+  print_tensor_float(s2t, "sub s2t of M0 from 3");
+
+  for(size_t i1=0; i1<d0->perm[1];++i1){
+    for(size_t i2=0; i2<d0->perm[2];++i2){
+      printf("EXPECT_EQ_TYPE_FLOAT c0[%ld](%ld,%ld,%ld) : s[%ld](%ld,%ld) \n",LineFromCoord((size_t[]){rnkId,i1,i2},M0->dim),rnkId,i1,i2,LineFromCoord((size_t[]){i1,i2},s2t->dim),i1,i2);
+      EXPECT_EQ_TYPE_FLOAT(M0->x[LineFromCoord((size_t[]){rnkId,i1,i2},M0->dim)],s2t->x[LineFromCoord((size_t[]){i1,i2},s2t->dim)]);
+    }  
+  }
+
+}
 
 TEST(tensorProd ){
   dimension *d0=create_dim(3);