tensor-0.1.0
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Enumerator Groups Pages
tensor_permute.cc
1 // -*- mode: c++; fill-column: 80; c-basic-offset: 2; indent-tabs-mode: nil -*-
2 /*
3  Copyright (c) 2010 Juan Jose Garcia Ripoll
4 
5  Tensor is free software; you can redistribute it and/or modify it
6  under the terms of the GNU Library General Public License as published
7  by the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Library General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with this program; if not, write to the Free Software Foundation, Inc.,
17  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #define TENSOR_LOAD_IMPL
21 #include <tensor/tensor.h>
22 
23 namespace tensor {
24 
25  template<typename n>
26  void permute_12(Tensor<n> &b, const Tensor<n> &a, index a1, index a2, index a3)
27  {
28  // In an abstract sense b(j,i,k) = a(i,j,k)
29  // Both tensors are stored in row-major order. So for A we have
30  // (0,0,0), (1,0,0), ... (2,1,0), ... (a1-1,a2-1,a3-1)
31  //
32  typename Tensor<n>::const_iterator ijk_a = a.begin();
33  typename Tensor<n>::iterator k_b = b.begin();
34  for (; a3; a3--, k_b += (a1*a2)) {
35  typename Tensor<n>::iterator jk_b = k_b;
36  for (index j = a2; j--; jk_b++) {
37  typename Tensor<n>::iterator ijk_b = jk_b;
38  for (index i = a1; i--; ijk_a++, ijk_b += a2) {
39  //assert(ijk_a >= a.begin() && ijk_a < a.end());
40  //assert(ijk_b >= b.begin() && ijk_b < b.end());
41  *ijk_b = *ijk_a;
42  }
43  }
44  }
45  }
46 
47  template<typename n>
48  void permute_23(Tensor<n> &b, const Tensor<n> &a, index a1, index a2, index a3, index a4)
49  {
50  // In an abstract sense
51  // b(i,k,j,l) = a(i,j,k,l)
52  // b(a1,a3,a2,a4) = a(a1,a2,a3,a4)
53  // Both tensors are stored in row-major order. So for A we have
54  // (0,0,0), (1,0,0), ... (2,1,0), ... (a1-1,a2-1,a3-1)
55  //
56  typename Tensor<n>::const_iterator ijkl_a = a.begin();
57  typename Tensor<n>::iterator l_b = b.begin();
58  index a13 = a1*a3;
59  index a132 = a13*a2;
60  for (; a4--; l_b += a132) {
61  typename Tensor<n>::iterator kl_b = l_b;
62  for (index k = a3; k--; kl_b += a1) {
63  typename Tensor<n>::iterator kjl_b = kl_b;
64  for (index j = a2; j--; kjl_b += a13) {
65  typename Tensor<n>::iterator ikjl_b = kjl_b;
66  for (index i = a1; i--; ijkl_a++, ikjl_b++) {
67  //assert(ijkl_a >= a.begin() && ijkl_a < a.end());
68  //assert(ikjl_b >= b.begin() && ikjl_b < b.end());
69  *ikjl_b = *ijkl_a;
70  }
71  }
72  }
73  }
74  }
75 
76  template<typename n>
77  void permute_24(Tensor<n> &b, const Tensor<n> &a, index a1, index a2, index a3,
78  index a4, index a5)
79  {
80  // In an abstract sense
81  // b(i,l,k,j,m) = a(i,j,k,l,m)
82  // b(a1,a4,a3,a2,a5) = a(a1,a2,a3,a4,a5)
83  // Both tensors are stored in row-major order. So for A we have
84  // (0,0,0), (1,0,0), ... (2,1,0), ... (a1-1,a2-1,a3-1)
85  //
86  typename Tensor<n>::const_iterator ijklm_a = a.begin();
87  typename Tensor<n>::iterator m_b = b.begin();
88  index a14 = a1*a4;
89  index a143 = a14*a3;
90  index a1432 = a143*a2;
91  for (; a5--; m_b += a1432) {
92  typename Tensor<n>::iterator lm_b = m_b;
93  for (index l = a4; l--; lm_b += a1) {
94  typename Tensor<n>::iterator lkm_b = lm_b;
95  for (index k = a3; k--; lkm_b += a14) {
96  typename Tensor<n>::iterator lkjm_b = lkm_b;
97  for (index j = a2; j--; lkjm_b += a143) {
98  typename Tensor<n>::iterator ilkjm_b = lkjm_b;
99  for (index i = a1; i--; ijklm_a++, ilkjm_b++) {
100  //assert(ilkjm_b < b.begin() || b.end() <= ilkjm_b);
101  //assert(ijklm_a < a.begin() || a.end() <= ijklm_a);
102  *ilkjm_b = *ijklm_a;
103  }
104  }
105  }
106  }
107  }
108  }
109 
110  template<typename n>
111  void permute_13(Tensor<n> &b, const Tensor<n> &a, index a1, index a2, index a3,
112  index a4)
113  {
114  // In an abstract sense
115  // b(k,j,i,l) = a(i,j,k,l)
116  // b(a3,a2,a1,a4) = a(a1,a2,a3,a4)
117  // Both tensors are stored in row-major order. So for A we have
118  // (0,0,0), (1,0,0), ... (2,1,0), ... (a1-1,a2-1,a3-1)
119  //
120  typename Tensor<n>::const_iterator ijkl_a = a.begin();
121  typename Tensor<n>::iterator l_b = b.begin();
122  index a32 = a3*a2;
123  index a321 = a32*a1;
124  for (; a4--; l_b += a321) {
125  typename Tensor<n>::iterator kl_b = l_b;
126  for (index k = a3; k--; kl_b++) {
127  typename Tensor<n>::iterator kjl_b = kl_b;
128  for (index j = a2; j--; kjl_b += a3) {
129  typename Tensor<n>::iterator kjil_b = kjl_b;
130  for (index i = a1; i--; ijkl_a++, kjil_b += a32) {
131  //assert(kjil_b < b.begin() || b.end() <= kjil_b);
132  //assert(ijkl_a < a.begin() || a.end() <= ijkl_a);
133  *kjil_b = *ijkl_a;
134  }
135  }
136  }
137  }
138  }
139 
140  template<typename n>
141  const Tensor<n> do_permute(const Tensor<n> &a, index ndx1, index ndx2)
142  {
143  index n1 = normalize_index(ndx1, a.rank());
144  index n2 = normalize_index(ndx2, a.rank());
145  if (n2 < n1) {
146  std::swap(n1,n2);
147  } else if (n2 == n1) {
148  return a;
149  }
150  Indices new_dims = a.dimensions();
151  index i,a1,a2,a3,a4,a5;
152  for (i = 0, a1 = 1; i < n1; )
153  a1 *= new_dims[i++];
154  a2 = new_dims[i++];
155  for (a3 = 1; i < n2;)
156  a3 *= new_dims[i++];
157  a4 = new_dims[i++];
158  for (a5 = 1; i < a.rank(); )
159  a5 *= new_dims[i++];
160 
161  std::swap(new_dims.at(n1), new_dims.at(n2));
162  Tensor<n> output(new_dims);
163 
164  if (output.size()) {
165  if (a3 > 1) {
166  if (a1 > 1) {
167  permute_24(output, a, a1,a2,a3,a4,a5);
168  } else {
169  permute_13(output, a, a2,a3,a4,a5);
170  }
171  } else if (a1 > 1) {
172  permute_23(output, a, a1,a2,a4,a5);
173  } else {
174  permute_12(output, a, a2,a4,a5);
175  }
176  }
177  return output;
178  }
179 
180 } // namespace tensor