Cyclops Tensor Framework
parallel arithmetic on multidimensional arrays
trace.cxx
Go to the documentation of this file.
1 /*Copyright (c) 2011, Edgar Solomonik, all rights reserved.*/
9 #include <ctf.hpp>
10 using namespace CTF;
11 
12 int trace(int const n,
13  World &dw){
14  int rank, i, num_pes;
15  int64_t np;
16  double * pairs;
17  double tr1, tr2, tr3, tr4;
18  int64_t * indices;
19 
20  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
21  MPI_Comm_size(MPI_COMM_WORLD, &num_pes);
22 
23  Matrix<> A(n, n, NS, dw);
24  Matrix<> B(n, n, NS, dw);
25  Matrix<> C(n, n, NS, dw);
26  Matrix<> D(n, n, NS, dw);
27  Matrix<> C1(n, n, NS, dw);
28  Matrix<> C2(n, n, NS, dw);
29  Matrix<> C3(n, n, NS, dw);
30  Matrix<> C4(n, n, NS, dw);
31  Vector<> DIAG(n, dw);
32 
33  srand48(13*rank);
34 
35  A.get_local_data(&np, &indices, &pairs);
36  for (i=0; i<np; i++ ) pairs[i] = drand48();;
37  A.write(np, indices, pairs);
38  delete [] pairs;
39  free(indices);
40  B.get_local_data(&np, &indices, &pairs);
41  for (i=0; i<np; i++ ) pairs[i] = drand48();
42  B.write(np, indices, pairs);
43  delete [] pairs;
44  free(indices);
45  C.get_local_data(&np, &indices, &pairs);
46  for (i=0; i<np; i++ ) pairs[i] = drand48();
47  C.write(np, indices, pairs);
48  delete [] pairs;
49  free(indices);
50  D.get_local_data(&np, &indices, &pairs);
51  for (i=0; i<np; i++ ) pairs[i] = drand48();
52  D.write(np, indices, pairs);
53  delete [] pairs;
54  free(indices);
55 
56 
57  C1["ij"] = A["ia"]*B["ab"]*C["bc"]*D["cj"];
58  C2["ij"] = D["ia"]*A["ab"]*B["bc"]*C["cj"];
59  C3["ij"] = C["ia"]*D["ab"]*A["bc"]*B["cj"];
60  C4["ij"] = B["ia"]*C["ab"]*D["bc"]*A["cj"];
61 
62 
63  DIAG["i"] = C1["ii"];
64  tr1 = DIAG.reduce(CTF::OP_SUM);
65 
66  DIAG["i"] = C2["ii"];
67  tr2 = DIAG.reduce(CTF::OP_SUM);
68  DIAG["i"] = C3["ii"];
69  tr3 = DIAG.reduce(CTF::OP_SUM);
70  DIAG["i"] = C4["ii"];
71  tr4 = DIAG.reduce(CTF::OP_SUM);
72 
73  int pass = 1;
74  if (rank == 0){
75  MPI_Reduce(MPI_IN_PLACE, &pass, 1, MPI_INT, MPI_MIN, 0, MPI_COMM_WORLD);
76  /*printf("tr(ABCD)=%lf, tr(DABC)=%lf, tr(CDAB)=%lf, tr(BCDA)=%lf\n",
77  tr1, tr2, tr3, tr4);*/
78  if (fabs(tr1-tr2)/tr1>1.E-10 || fabs(tr2-tr3)/tr2>1.E-10 || fabs(tr3-tr4)/tr3>1.E-10){
79  pass = 0;
80  }
81  if (!pass){
82  printf("{ tr(ABCD) = tr(DABC) = tr(CDAB) = tr(BCDA) } failed\n");
83  } else {
84  printf("{ tr(ABCD) = tr(DABC) = tr(CDAB) = tr(BCDA) } passed\n");
85  }
86  } else
87  MPI_Reduce(&pass, MPI_IN_PLACE, 1, MPI_INT, MPI_MIN, 0, MPI_COMM_WORLD);
88  return pass;
89 
90 }
91 
92 
93 #ifndef TEST_SUITE
94 char* getCmdOption(char ** begin,
95  char ** end,
96  const std::string & option){
97  char ** itr = std::find(begin, end, option);
98  if (itr != end && ++itr != end){
99  return *itr;
100  }
101  return 0;
102 }
103 
104 
105 int main(int argc, char ** argv){
106  int rank, np, n;
107  int const in_num = argc;
108  char ** input_str = argv;
109 
110  MPI_Init(&argc, &argv);
111  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
112  MPI_Comm_size(MPI_COMM_WORLD, &np);
113 
114  if (getCmdOption(input_str, input_str+in_num, "-n")){
115  n = atoi(getCmdOption(input_str, input_str+in_num, "-n"));
116  if (n < 0) n = 7;
117  } else n = 7;
118 
119 
120  {
121  World dw(MPI_COMM_WORLD, argc, argv);
122  if (rank == 0){
123  printf("Checking trace calculation n = %d, p = %d:\n",n,np);
124  }
125  int pass = trace(n,dw);
126  assert(pass);
127  }
128 
129  MPI_Finalize();
130  return 0;
131 }
137 #endif
Matrix class which encapsulates a 2D tensor.
Definition: matrix.h:18
def rank(self)
Definition: core.pyx:312
Vector class which encapsulates a 1D tensor.
Definition: vector.h:14
dtype reduce(OP op)
performs a reduction on the tensor
Definition: tensor.cxx:731
int trace(int const n, World &dw)
Definition: trace.cxx:12
Definition: common.h:37
an instance of the CTF library (world) on a MPI communicator
Definition: world.h:19
string
Definition: core.pyx:456
void get_local_data(int64_t *npair, int64_t **global_idx, dtype **data, bool nonzeros_only=false, bool unpack_sym=false) const
Gives the global indices and values associated with the local data.
Definition: tensor.cxx:159
Definition: apsp.cxx:17
char * getCmdOption(char **begin, char **end, const std::string &option)
Definition: trace.cxx:94
int main(int argc, char **argv)
Definition: trace.cxx:105
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
def np(self)
Definition: core.pyx:315