Cyclops Tensor Framework parallel arithmetic on multidimensional arrays
decomposition.cxx
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include "decomposition.h"
4
5 namespace CTF {
6
7  template<typename dtype>
9
10  int64_t * inds_X;
11  double * vals_X;
12  int64_t n_X;
13  //if global index ordering is preserved between the two tensors, we can fold simply
15  Y.write(n_X, inds_X, vals_X);
16
17  }
18
19
20  template<typename dtype>
21  Contract_Term HoSVD::operator[](char const * idx_map){
22  char int_inds[core_tensor.order];
23  for (int i=0; i<core_tensor.order; i++){
24  int_inds = "["+i;
25  //FIXME: how to make this robust?
26  }
27  char factor_inds[2];
28  factor_inds[0] = idx_map[0];
29  factor_inds[1] = int_inds[0];
30  Contract_Term t(core_tensor[ind_inds],factor_matrices[0][factor_inds]);
31  for (int i=1; i<core_tensor.order; i++){
32  factor_inds[0] = idx_map[i];
33  factor_inds[1] = int_inds[i];
34  Contract_Term t(t,factor_matrices[0][factor_inds]);
35  }
36  return t;
37  }
38
39
44  template<typename dtype>
45  std::vector< Matrix <dtype> > get_factor_matrices(Tensor<dtype>& T, int * ranks) {
46
47  std::vector< Matrix <dtype> > factor_matrices (T.order);
48
49  char arg[T.order];
50  int transformed_lens[T.order];
51  char transformed_arg[T.order];
52  for (int i = 0; i < T.order; i++) {
53  arg[i] = 'i' + i;
54  transformed_arg[i] = 'i' + i;
55  transformed_lens[i] = T.lens[i];
56  }
57
58  for (int i = 0; i < T.order; i++) {
59  for (int j = i; j > 0; j--) {
60  transformed_lens[j] = T.lens[j-1];
61  }
62
63  transformed_lens[0] = T.lens[i];
64  for (int j = 0; j < i; j++) {
65  transformed_arg[j] = arg[j+1];
66  }
67  transformed_arg[i] = arg[0];
68
69  int unfold_lens [2];
70  unfold_lens[0] = T.lens[i];
71  int ncol = 1;
72
73  for (int j = 0; j < T.order; j++) {
74  if (j != i)
75  ncol *= T.lens[j];
76  }
77  unfold_lens[1] = ncol;
78
79  Tensor<dtype> transformed_T(T.order, transformed_lens, *(T.wrld));
80  transformed_T[arg] = T[transformed_arg];
81
82  Tensor<dtype> cur_unfold(2, unfold_lens, *(T.wrld));
83  fold_unfold(transformed_T, cur_unfold);
84
85  Matrix<dtype> M(cur_unfold);
86  Matrix<dtype> U;
87  Matrix<dtype> VT;
88  Vector<dtype> S;
89  M.svd(U, S, VT, ranks[i]);
90  factor_matrices[i] = U;
91
92  }
93
94  return factor_matrices;
95  }
96  template<typename dtype>
98
99  std::vector< Tensor <> > core_tensors(T.order+1);
100  core_tensors[0] = T;
101  int lens[T.order];
102  for (int i = 0; i < T.order; i++) {
103  lens[i] = T.lens[i];
104  }
105  for (int i = 1; i < T.order+1; i++) {
106  lens[i-1] = ranks[i-1];
107  Tensor<dtype> core(T.order, lens, *(T.wrld));
108  core_tensors[i] = core;
109  }
110
111  //calculate core tensor
112  char arg[T.order];
113  char core_arg[T.order];
114  for (int i = 0; i < T.order; i++) {
115  arg[i] = 'i' + i;
116  core_arg[i] = 'i' + i;
117  }
118  char matrix_arg[2];
119  matrix_arg[0] = 'a';
120  for (int i = 0; i < T.order; i++) {
121  core_arg[i] = 'a';
122  matrix_arg[1] = arg[i];
124  transpose["ij"] = factor_matrices[i]["ji"];
125  core_tensors[i+1][core_arg] = transpose[matrix_arg] * core_tensors[i][arg];
126  core_arg[i] = arg[i];
127  }
128  return core_tensors[T.order];
129  }
130
131  template<typename dtype>
132  void HoSVD::HoSVD(Tensor<dtype>& T, int * ranks) {
135  }
136
137  template<typename dtype>
138  void HoSVD::HoSVD(int * lens, int * ranks) {
139  //initialize to zero
140  assert(0);
141  /*factor_matrices = get_factor_matrices(T, ranks);
142  core_tensor = get_core_tensor(T, factor_matrices, ranks); */
143  }
144
145 }
146
147
148
Tensor< dtype > core_tensor
Definition: decomposition.h:23
Contract_Term operator[](char const *idx_map)
associated an index map with the tensor decomposition for algebra
Matrix class which encapsulates a 2D tensor.
Definition: matrix.h:18
std::vector< Matrix< dtype > > factor_matrices
Definition: decomposition.h:24
Vector class which encapsulates a 1D tensor.
Definition: vector.h:14
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
CTF::World * wrld
distributed processor context on which tensor is defined
int * lens
void svd(Matrix< dtype > &U, Vector< dtype > &S, Matrix< dtype > &VT, int rank=0)
Definition: matrix.cxx:481
Tensor< dtype > get_core_tensor(Tensor< dtype > &T, std::vector< Matrix< dtype > > factor_matrices, int *ranks)
def transpose(init_A, axes=None)
Definition: core.pyx:4850
std::vector< Matrix< dtype > > get_factor_matrices(Tensor< dtype > &T, int *ranks)
HoSVD(Tensor< dtype > T, int *ranks)
Definition: apsp.cxx:17
an instance of a tensor within a CTF world
Definition: tensor.h:74
void fold_unfold(Tensor< dtype > &X, Tensor< dtype > &Y)
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