Cyclops Tensor Framework
parallel arithmetic on multidimensional arrays
hosvd.cxx
Go to the documentation of this file.
1 #include <ctf.hpp>
2 #include <stdlib.h>
3 #include <iostream>
4 #include <stdio.h>
5 #include <complex>
6 
7 int icontxt;
8 
9 
10 extern "C" {
11 
12  void Cblacs_pinfo(int*, int*);
13 
14  void Cblacs_get(int, int, int*);
15 
16  void Cblacs_gridinit(int*, char*, int, int);
17 
18  void Cblacs_gridinfo(int, int*, int*, int*, int*);
19 
20  void Cblacs_gridmap(int*, int*, int, int, int);
21 
22  void Cblacs_barrier(int , char*);
23 
24  void Cblacs_gridexit(int);
25 
26 }
27 
28 
29 
30 using namespace CTF;
31 
32 
34 
35  int64_t * inds_X;
36  double * vals_X;
37  int64_t n_X;
38  //if global index ordering is preserved between the two tensors, we can fold simply
39  X.read_local(&n_X, &inds_X, &vals_X);
40  Y.write(n_X, inds_X, vals_X);
41 
42 }
43 
44 std::vector< Matrix <> > get_factor_matrices(Tensor<>& T, int ranks[], World& dw) {
45 
46  std::vector< Matrix <> > factor_matrices (T.order);
47 
48  char chars[] = {'i','j','k','l','m','n','o','p','\0'};
49  char arg[T.order+1];
50  int transformed_lens[T.order];
51  char transformed_arg[T.order+1];
52  transformed_arg[T.order] = '\0';
53  for (int i = 0; i < T.order; i++) {
54  arg[i] = chars[i];
55  transformed_arg[i] = chars[i];
56  transformed_lens[i] = T.lens[i];
57  }
58  arg[T.order] = '\0';
59 
60 
61  for (int i = 0; i < T.order; i++) {
62  for (int j = i; j > 0; j--) {
63  transformed_lens[j] = T.lens[j-1];
64  }
65 
66  transformed_lens[0] = T.lens[i];
67  for (int j = 0; j < i; j++) {
68  transformed_arg[j] = arg[j+1];
69  }
70  transformed_arg[i] = arg[0];
71 
72  int unfold_lens [2];
73  unfold_lens[0] = T.lens[i];
74  int ncol = 1;
75 
76  for (int j = 0; j < T.order; j++) {
77  if (j != i)
78  ncol *= T.lens[j];
79  }
80  unfold_lens[1] = ncol;
81 
82  Tensor<double> transformed_T(T.order, transformed_lens, dw);
83  transformed_T[arg] = T[transformed_arg];
84 
85  Tensor<double> cur_unfold(2, unfold_lens, dw);
86  fold_unfold(transformed_T, cur_unfold);
87 
88  Matrix<double> M(cur_unfold);
89  Matrix<> U;
90  Matrix<> VT;
91  Vector<> S;
92  M.svd(U, S, VT, ranks[i]);
93 /*
94  printf("%d-mode unfolding of T:\n", i+1);
95  M.print_matrix();
96  printf("SVD of %d-mode unfolding of T\n", i+1);
97  printf("Left singular vectors (U): \n");
98  U.print_matrix();
99  printf("Singular values (S):\n");
100  S.print();
101  printf("Right singular vectors (VT):\n");
102  VT.print_matrix();
103 */
104  factor_matrices[i] = U;
105 
106  }
107 
108  return factor_matrices;
109 }
110 
111 Tensor<> get_core_tensor(Tensor<>& T, std::vector< Matrix <> > factor_matrices, int ranks[], World& dw) {
112 
113  std::vector< Tensor <> > core_tensors(T.order+1);
114  core_tensors[0] = T;
115  int lens[T.order];
116  for (int i = 0; i < T.order; i++) {
117  lens[i] = T.lens[i];
118  }
119  for (int i = 1; i < T.order+1; i++) {
120  lens[i-1] = ranks[i-1];
121  Tensor<double> core(T.order, lens, dw);
122  core_tensors[i] = core;
123  }
124 
125  //calculate core tensor
126  char chars[] = {'i','j','k','l','m','n','o','p','\0'};
127  char arg[T.order+1];
128  char core_arg[T.order+1];
129  for (int i = 0; i < T.order; i++) {
130  arg[i] = chars[i];
131  core_arg[i] = chars[i];
132  }
133  arg[T.order] = '\0';
134  core_arg[T.order] = '\0';
135  char matrix_arg[3];
136  matrix_arg[0] = 'a';
137  matrix_arg[2] = '\0';
138  for (int i = 0; i < T.order; i++) {
139  core_arg[i] = 'a';
140  matrix_arg[1] = arg[i];
141  Matrix<double> transpose(factor_matrices[i].ncol, factor_matrices[i].nrow, dw);
142  transpose["ij"] = factor_matrices[i]["ji"];
143  //printf("transpose of factor matrix %d is\n", i+1);
144  //transpose.print_matrix();
145  //printf("core_arg is %s \n", core_arg);
146  //printf("matrix_arg is %s \n", matrix_arg);
147  //printf("arg is %s \n", arg);
148  core_tensors[i+1][core_arg] = transpose[matrix_arg] * core_tensors[i][arg];
149  core_arg[i] = arg[i];
150  }
151  return core_tensors[T.order];
152 }
153 
154 void hosvd(Tensor<>& T, Tensor<>& core, std::vector< Matrix <> > factor_matrices, int * ranks, World& dw) {
155  factor_matrices = get_factor_matrices(T, ranks, dw);
156  core = Tensor<double>(get_core_tensor(T, factor_matrices, ranks, dw));
157 }
158 
159 Tensor<> get_core_tensor_hooi(Tensor<>& T, std::vector< Matrix <> > factor_matrices, int ranks[], World& dw, int j = -1) {
160  printf("j is %d \n", j);
161 
162  std::vector< Tensor <> > core_tensors(T.order);
163  core_tensors[0] = T;
164  int lens[T.order];
165  for (int i = 0; i < T.order; i++) {
166  lens[i] = T.lens[i];
167  }
168  for (int i = 1; i < T.order+1; i++) {
169  if (i < j+1) {
170  lens[i-1] = ranks[i-1];
171  Tensor<double> core(T.order, lens, dw);
172  core_tensors[i] = core;
173  }
174  if (i < j+1) {
175  lens[i-1] = ranks[i-1];
176  Tensor<double> core(T.order, lens, dw);
177  core_tensors[i-1] = core;
178  }
179  }
180 
181  //calculate core tensor
182  char chars[] = {'i','j','k','l','m','n','o','p','\0'};
183  char arg[T.order+1];
184  char core_arg[T.order+1];
185  for (int i = 0; i < T.order; i++) {
186  arg[i] = chars[i];
187  core_arg[i] = chars[i];
188  }
189  arg[T.order] = '\0';
190  core_arg[T.order] = '\0';
191  char matrix_arg[3];
192  matrix_arg[0] = 'a';
193  matrix_arg[2] = '\0';
194  for (int i = 0; i < T.order; i++) {
195  if (i != j) {
196  core_arg[i] = 'a';
197  matrix_arg[1] = arg[i];
198  Matrix<double> transpose(factor_matrices[i].ncol, factor_matrices[i].nrow, dw);
199  transpose["ij"] = factor_matrices[i]["ji"];
200  //printf("transpose of factor matrix %d is\n", i+1);
201  //transpose.print_matrix();
202  //printf("core_arg is %s \n", core_arg);
203  //printf("matrix_arg is %s \n", matrix_arg);
204  //printf("arg is %s \n", arg);
205  if (i < j) {
206  core_tensors[i+1][core_arg] = transpose[matrix_arg] * core_tensors[i][arg];
207  core_arg[i] = arg[i];
208  }
209  if (i > j) {
210  core_tensors[i+1][core_arg] = transpose[matrix_arg] * core_tensors[i-1][arg];
211  core_arg[i] = arg[i];
212  }
213  }
214  }
215  return core_tensors[T.order-1];
216 }
217 
218 void hooi(Tensor<>& T, Tensor<>& core, std::vector< Matrix <> > factor_matrices, int * ranks, World& dw) {
219  factor_matrices = get_factor_matrices(T, ranks, dw);
220 
221  for (int i = 0; i < T.order; i++) {
222  Tensor<double> temp_core = get_core_tensor_hooi(T, factor_matrices, ranks, dw, i);
223  printf("here");
224  std::vector< Matrix<double> > temp_factor_matrices = get_factor_matrices(temp_core, ranks, dw);
225  factor_matrices[i] = temp_factor_matrices[i];
226  }
227 }
228 
229 int main(int argc, char ** argv) {
230 
231  int rank, np;
232  char cC = 'C';
233 
234  MPI_Init(&argc, &argv);
235  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
236  MPI_Comm_size(MPI_COMM_WORLD, &np);
237  World dw(argc, argv);
238 
239  Cblacs_get(-1, 0, &icontxt);
240  Cblacs_gridinit(&icontxt, &cC, np, 1);
241 
242 /*
243  int T_lens[] = {2 ,2 ,3};
244  Tensor<double> T(3, T_lens, dw);
245  T.fill_random(0,10);
246  printf("Tensor T \n");
247  T.print();
248 
249  int ranks[] = {2,1,2};
250 
251  std::vector< Matrix<> > hosvd_factor_matrices;
252  Tensor<> hosvd_core;
253  hosvd(T, hosvd_core, hosvd_factor_matrices, ranks, dw);
254  hosvd_core.print();
255 
256  std::vector< Matrix<> > hooi_factor_matrices;
257  Tensor<> hooi_core;
258  printf("here");
259  hooi(T, hooi_core, hooi_factor_matrices, ranks, dw);
260  hooi_core.print();
261 */
262 
263  int lens[] = {4, 4};
264  Tensor<double> T2(2, lens, dw);
265  T2.fill_random(0,10);
266  Matrix<double> M(T2);
267  M.print_matrix();
268  Matrix<> U;
269  Matrix<> VT;
270  Vector<> S;
271  M.svd(U, S, VT);
272  U.print_matrix();
273 
274  T2.print();
275 
276 
277 
278  MPI_Finalize();
279 
280  return 0;
281 }
282 
283 
void Cblacs_gridinit(int *, char *, int, int)
Matrix class which encapsulates a 2D tensor.
Definition: matrix.h:18
def rank(self)
Definition: core.pyx:312
void Cblacs_get(int, int, int *)
int icontxt
Definition: hosvd.cxx:7
Vector class which encapsulates a 1D tensor.
Definition: vector.h:14
an instance of the CTF library (world) on a MPI communicator
Definition: world.h:19
int order
number of tensor dimensions
void read_local(int64_t *npair, int64_t **global_idx, dtype **data, bool unpack_sym=false) const
Using get_local_data(), which returns an array that must be freed with delete [], is more efficient...
Definition: tensor.cxx:182
void fill_random(dtype rmin, dtype rmax)
fills local unique tensor elements to random values in the range [min,max] works only for dtype in {f...
Definition: tensor.cxx:928
void Cblacs_gridinfo(int, int *, int *, int *, int *)
void Cblacs_gridexit(int)
int * lens
unpadded tensor edge lengths
void svd(Matrix< dtype > &U, Vector< dtype > &S, Matrix< dtype > &VT, int rank=0)
Definition: matrix.cxx:481
int main(int argc, char **argv)
Definition: hosvd.cxx:229
void hosvd(Tensor<> &T, Tensor<> &core, std::vector< Matrix<> > factor_matrices, int *ranks, World &dw)
Definition: hosvd.cxx:154
void print(FILE *fp, dtype cutoff) const
prints tensor data to file using process 0 (modify print(...) overload in set.h if you would like a d...
Definition: tensor.cxx:407
Tensor< dtype > get_core_tensor(Tensor< dtype > &T, std::vector< Matrix< dtype > > factor_matrices, int *ranks)
void Cblacs_gridmap(int *, int *, int, int, int)
def transpose(init_A, axes=None)
Definition: core.pyx:4850
void Cblacs_barrier(int, char *)
std::vector< Matrix< dtype > > get_factor_matrices(Tensor< dtype > &T, int *ranks)
Definition: apsp.cxx:17
Tensor get_core_tensor_hooi(Tensor<> &T, std::vector< Matrix<> > factor_matrices, int ranks[], World &dw, int j=-1)
Definition: hosvd.cxx:159
an instance of a tensor within a CTF world
Definition: tensor.h:74
void hooi(Tensor<> &T, Tensor<> &core, std::vector< Matrix<> > factor_matrices, int *ranks, World &dw)
Definition: hosvd.cxx:218
void fold_unfold(Tensor< dtype > &X, Tensor< dtype > &Y)
void print_matrix()
Definition: matrix.cxx:144
void write(int64_t npair, int64_t const *global_idx, dtype const *data)
writes in values associated with any set of indices The sparse data is defined in coordinate format...
Definition: tensor.cxx:264
void Cblacs_pinfo(int *, int *)
def np(self)
Definition: core.pyx:315