Cyclops Tensor Framework parallel arithmetic on multidimensional arrays
ccsd.cxx
Go to the documentation of this file.
2
10 #include <ctf.hpp>
11 using namespace CTF;
12
13 double divide(double a, double b){
14  return a/b;
15 }
16
17
18 class Integrals {
19  public:
36
37  Integrals(int no, int nv, World &dw_){
38  int shapeASAS[] = {AS,NS,AS,NS};
39  int shapeASNS[] = {AS,NS,NS,NS};
40  int shapeNSNS[] = {NS,NS,NS,NS};
41  int shapeNSAS[] = {NS,NS,AS,NS};
42  int vvvv[] = {nv,nv,nv,nv};
43  int vvvo[] = {nv,nv,nv,no};
44  int vovv[] = {nv,no,nv,nv};
45  int vovo[] = {nv,no,nv,no};
46  int vvoo[] = {nv,nv,no,no};
47  int oovv[] = {no,no,nv,nv};
48  int vooo[] = {nv,no,no,no};
49  int oovo[] = {no,no,nv,no};
50  int oooo[] = {no,no,no,no};
51
52  dw = &dw_;
53
54  aa = new CTF_Vector(nv,dw_);
55  ii = new CTF_Vector(no,dw_);
56
57  ab = new CTF_Matrix(nv,nv,AS,dw_,"Vab",1);
58  ai = new CTF_Matrix(nv,no,NS,dw_,"Vai",1);
59  ia = new CTF_Matrix(no,nv,NS,dw_,"Via",1);
60  ij = new CTF_Matrix(no,no,AS,dw_,"Vij",1);
61
62  abcd = new Tensor<>(4,vvvv,shapeASAS,dw_,"Vabcd",1);
63  abci = new Tensor<>(4,vvvo,shapeASNS,dw_,"Vabci",1);
64  aibc = new Tensor<>(4,vovv,shapeNSAS,dw_,"Vaibc",1);
65  aibj = new Tensor<>(4,vovo,shapeNSNS,dw_,"Vaibj",1);
66  abij = new Tensor<>(4,vvoo,shapeASAS,dw_,"Vabij",1);
67  ijab = new Tensor<>(4,oovv,shapeASAS,dw_,"Vijab",1);
68  aijk = new Tensor<>(4,vooo,shapeNSAS,dw_,"Vaijk",1);
69  ijak = new Tensor<>(4,oovo,shapeASNS,dw_,"Vijak",1);
70  ijkl = new Tensor<>(4,oooo,shapeASAS,dw_,"Vijkl",1);
71  }
72
74  delete aa;
75  delete ii;
76
77  delete ab;
78  delete ai;
79  delete ia;
80  delete ij;
81
82  delete abcd;
83  delete abci;
84  delete aibc;
85  delete aibj;
86  delete abij;
87  delete ijab;
88  delete aijk;
89  delete ijak;
90  delete ijkl;
91  }
92
93  void fill_rand(){
94  int i, rank;
95  int64_t j, sz, * indices;
96  double * values;
97
98  Tensor<> * tarr[] = {aa, ii, ab, ai, ia, ij,
99  abcd, abci, aibc, aibj,
100  abij, ijab, aijk, ijak, ijkl};
101  MPI_Comm comm = dw->comm;
102  MPI_Comm_rank(comm, &rank);
103
104  srand48(rank*13);
105
106  for (i=0; i<15; i++){
107  tarr[i]->get_local_data(&sz, &indices, &values);
108 // for (j=0; j<sz; j++) values[j] = drand48()-.5;
109  for (j=0; j<sz; j++) values[j] = ((indices[j]*16+i)%13077)/13077. -.5;
110  tarr[i]->write(sz, indices, values);
111  free(indices), delete [] values;
112  }
113  }
114
115  Idx_Tensor operator[](char const * idx_map_){
116  int i, lenm, no, nv;
117  lenm = strlen(idx_map_);
118  char new_idx_map[lenm+1];
119  new_idx_map[lenm]='\0';
120  no = 0;
121  nv = 0;
122  for (i=0; i<lenm; i++){
123  if (idx_map_[i] >= 'a' && idx_map_[i] <= 'h'){
124  new_idx_map[i] = 'a'+nv;
125  nv++;
126  } else if (idx_map_[i] >= 'i' && idx_map_[i] <= 'n'){
127  new_idx_map[i] = 'i'+no;
128  no++;
129  }
130  }
131 // printf("indices %s are %s\n",idx_map_,new_idx_map);
132  if (0 == strcmp("a",new_idx_map)) return (*aa)[idx_map_];
133  if (0 == strcmp("i",new_idx_map)) return (*ii)[idx_map_];
134  if (0 == strcmp("ab",new_idx_map)) return (*ab)[idx_map_];
135  if (0 == strcmp("ai",new_idx_map)) return (*ai)[idx_map_];
136  if (0 == strcmp("ia",new_idx_map)) return (*ia)[idx_map_];
137  if (0 == strcmp("ij",new_idx_map)) return (*ij)[idx_map_];
138  if (0 == strcmp("abcd",new_idx_map)) return (*abcd)[idx_map_];
139  if (0 == strcmp("abci",new_idx_map)) return (*abci)[idx_map_];
140  if (0 == strcmp("aibc",new_idx_map)) return (*aibc)[idx_map_];
141  if (0 == strcmp("aibj",new_idx_map)) return (*aibj)[idx_map_];
142  if (0 == strcmp("abij",new_idx_map)) return (*abij)[idx_map_];
143  if (0 == strcmp("ijab",new_idx_map)) return (*ijab)[idx_map_];
144  if (0 == strcmp("aijk",new_idx_map)) return (*aijk)[idx_map_];
145  if (0 == strcmp("ijak",new_idx_map)) return (*ijak)[idx_map_];
146  if (0 == strcmp("ijkl",new_idx_map)) return (*ijkl)[idx_map_];
147  printf("Invalid integral indices\n");
148  assert(0);
149 //shut up compiler
150  return (*aa)[idx_map_];
151  }
152 };
153
154 class Amplitudes {
155  public:
159
160  Amplitudes(int no, int nv, World &dw_){
161  dw = &dw_;
162  int shapeASAS[] = {AS,NS,AS,NS};
163  int vvoo[] = {nv,nv,no,no};
164
165  ai = new CTF_Matrix(nv,no,NS,dw_,"Tai",1);
166  abij = new Tensor<>(4,vvoo,shapeASAS,dw_,"Tabij",1);
167  }
168
170  delete ai;
171  delete abij;
172  }
173
174  Idx_Tensor operator[](char const * idx_map_){
175  if (strlen(idx_map_) == 4) return (*abij)[idx_map_];
176  else return (*ai)[idx_map_];
177  }
178
179  void fill_rand(){
180  int i, rank;
181  int64_t j, sz, * indices;
182  double * values;
183
184  Tensor<> * tarr[] = {ai, abij};
185  MPI_Comm comm = dw->comm;
186  MPI_Comm_rank(comm, &rank);
187
188  srand48(rank*25);
189
190  for (i=0; i<2; i++){
191  tarr[i]->get_local_data(&sz, &indices, &values);
192 // for (j=0; j<sz; j++) values[j] = drand48()-.5;
193  for (j=0; j<sz; j++) values[j] = ((indices[j]*13+i)%13077)/13077. -.5;
194  tarr[i]->write(sz, indices, values);
195  free(indices), delete [] values;
196  }
197  }
198 };
199
200 void ccsd(Integrals &V,
201  Amplitudes &T,
202  int sched_nparts = 0){
203  int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
204 #ifdef SCHEDULE_CCSD
205  double timer = MPI_Wtime();
206  tCTF_Schedule<double> sched(V.dw);
207  sched.set_max_partitions(sched_nparts);
208  sched.record();
209 #endif
210
211  Tensor<> T21 = Tensor<>(T.abij);
212  T21["abij"] += .5*T["ai"]*T["bj"];
213
214  Tensor<> tFme(*V["me"].parent);
215  Idx_Tensor Fme(&tFme,"me");
216  Fme += V["me"];
217  Fme += V["mnef"]*T["fn"];
218
219  Tensor<> tFae(*V["ae"].parent);
220  Idx_Tensor Fae(&tFae,"ae");
221  Fae += V["ae"];
222  Fae -= Fme*T["am"];
223  Fae -=.5*V["mnef"]*T["afmn"];
224  Fae += V["anef"]*T["fn"];
225
226  Tensor<> tFmi(*V["mi"].parent);
227  Idx_Tensor Fmi(&tFmi,"mi");
228  Fmi += V["mi"];
229  Fmi += Fme*T["ei"];
230  Fmi += .5*V["mnef"]*T["efin"];
231  Fmi += V["mnfi"]*T["fn"];
232
233  Tensor<> tWmnei(*V["mnei"].parent);
234  Idx_Tensor Wmnei(&tWmnei,"mnei");
235  Wmnei += V["mnei"];
236  Wmnei += V["mnei"];
237  Wmnei += V["mnef"]*T["fi"];
238
239  Tensor<> tWmnij(*V["mnij"].parent);
240  Idx_Tensor Wmnij(&tWmnij,"mnij");
241  Wmnij += V["mnij"];
242  Wmnij -= V["mnei"]*T["ej"];
243  Wmnij += V["mnef"]*T21["efij"];
244
245  Tensor<> tWamei(*V["amei"].parent);
246  Idx_Tensor Wamei(&tWamei,"amei");
247  Wamei += V["amei"];
248  Wamei -= Wmnei*T["an"];
249  Wamei += V["amef"]*T["fi"];
250  Wamei += .5*V["mnef"]*T["afin"];
251
252  Tensor<> tWamij(*V["amij"].parent);
253  Idx_Tensor Wamij(&tWamij,"amij");
254  Wamij += V["amij"];
255  Wamij += V["amei"]*T["ej"];
256  Wamij += V["amef"]*T["efij"];
257
258  Tensor<> tZai(*V["ai"].parent);
259  Idx_Tensor Zai(&tZai,"ai");
260  Zai += V["ai"];
261  Zai -= Fmi*T["am"];
262  Zai += V["ae"]*T["ei"];
263  Zai += V["amei"]*T["em"];
264  Zai += V["aeim"]*Fme;
265  Zai += .5*V["amef"]*T21["efim"];
266  Zai -= .5*Wmnei*T21["eamn"];
267
268  Tensor<> tZabij(*V["abij"].parent);
269  Idx_Tensor Zabij(&tZabij,"abij");
270  Zabij += V["abij"];
271  Zabij += V["abei"]*T["ej"];
272  Zabij += Wamei*T["ebmj"];
273  Zabij -= Wamij*T["bm"];
274  Zabij += Fae*T["ebij"];
275  Zabij -= Fmi*T["abmj"];
276  Zabij += .5*V["abef"]*T21["efij"];
277  Zabij += .5*Wmnij*T21["abmn"];
278
279 #ifdef SCHEDULE_CCSD
280  if (rank == 0) {
281  printf("Record: %lf\n",
282  MPI_Wtime()-timer);
283  }
284
285  timer = MPI_Wtime();
286  tCTF_ScheduleTimer schedule_time = sched.execute();
287 #endif
288
289
290  Tensor<> Dai(2, V.ai->lens, V.ai->sym, *V.dw);
291  int sh_sym[4] = {SH, NS, SH, NS};
292  Tensor<> Dabij(4, V.abij->lens, sh_sym, *V.dw);
293  Dai["ai"] += V["i"];
294  Dai["ai"] -= V["a"];
295
296  Dabij["abij"] += V["i"];
297  Dabij["abij"] += V["j"];
298  Dabij["abij"] -= V["a"];
299  Dabij["abij"] -= V["b"];
300
301
302  Function<> fctr(&divide);
303
304  T.ai->contract(1.0, *(Zai.parent), "ai", Dai, "ai", 0.0, "ai", fctr);
305  T.abij->contract(1.0, *(Zabij.parent), "abij", Dabij, "abij", 0.0, "abij", fctr);
306 #ifdef SCHEDULE_CCSD
307  if (rank == 0) {
308  printf("Schedule comm down: %lf\n", schedule_time.comm_down_time);
309  printf("Schedule execute: %lf\n", schedule_time.exec_time);
310  printf("Schedule imbalance, wall: %lf\n", schedule_time.imbalance_wall_time);
311  printf("Schedule imbalance, accum: %lf\n", schedule_time.imbalance_acuum_time);
312  printf("Schedule comm up: %lf\n", schedule_time.comm_up_time);
313  printf("Schedule total: %lf\n", schedule_time.total_time);
314  printf("All execute: %lf\n",
315  MPI_Wtime()-timer);
316  }
317 #endif
318 }
319
320 #ifndef TEST_SUITE
321
322 char* getCmdOption(char ** begin,
323  char ** end,
324  const std::string & option){
325  char ** itr = std::find(begin, end, option);
326  if (itr != end && ++itr != end){
327  return *itr;
328  }
329  return 0;
330 }
331
332
333 int main(int argc, char ** argv){
334  int rank, np, niter, no, nv, sched_nparts, i;
335  int const in_num = argc;
336  char ** input_str = argv;
337
338  MPI_Init(&argc, &argv);
339  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
340  MPI_Comm_size(MPI_COMM_WORLD, &np);
341
342  if (getCmdOption(input_str, input_str+in_num, "-no")){
343  no = atoi(getCmdOption(input_str, input_str+in_num, "-no"));
344  if (no < 0) no= 4;
345  } else no = 4;
346  if (getCmdOption(input_str, input_str+in_num, "-nv")){
347  nv = atoi(getCmdOption(input_str, input_str+in_num, "-nv"));
348  if (nv < 0) nv = 6;
349  } else nv = 6;
350  if (getCmdOption(input_str, input_str+in_num, "-niter")){
351  niter = atoi(getCmdOption(input_str, input_str+in_num, "-niter"));
352  if (niter < 0) niter = 1;
353  } else niter = 1;
354  if (getCmdOption(input_str, input_str+in_num, "-nparts")){
355  sched_nparts = atoi(getCmdOption(input_str, input_str+in_num, "-nparts"));
356  if (sched_nparts < 0) sched_nparts = 0;
357  } else sched_nparts = 0;
358
359  {
360  World dw(argc, argv);
361  {
362  Integrals V(no, nv, dw);
363  V.fill_rand();
364  Amplitudes T(no, nv, dw);
365  Timer_epoch tccsd("CCSD");
366  tccsd.begin();
367  for (i=0; i<niter; i++){
368  T.fill_rand();
369  double d = MPI_Wtime();
370  ccsd(V,T,sched_nparts);
371  if (rank == 0)
372  printf("(%d nodes) Completed %dth CCSD iteration in time = %lf, |T| is %lf\n",
373  np, i, MPI_Wtime()-d, T.ai->norm2()+T.abij->norm2());
374  else {
375  T.ai->norm2();
376  T.abij->norm2();
377  }
378  T["ai"] = (1./T.ai->norm2())*T["ai"];
379  T["abij"] = (1./T.abij->norm2())*T["abij"];
380  }
381  tccsd.end();
382  }
383  }
384
385  MPI_Finalize();
386  return 0;
387 }
393 #endif
394
Tensor * ijkl
Definition: ccsd.cxx:35
Idx_Tensor operator[](char const *idx_map_)
Definition: ccsd.cxx:174
Tensor * abij
Definition: ccsd.cxx:31
void fill_rand()
Definition: ccsd.cxx:179
Tensor * aibj
Definition: ccsd.cxx:30
MPI_Comm comm
Definition: int_timer.cxx:22
Tensor * ai
Definition: ccsd.cxx:24
int * sym
symmetries among tensor dimensions
CTF::Matrix CTF_Matrix
Definition: back_comp.h:39
Tensor * ijak
Definition: ccsd.cxx:34
Idx_Tensor operator[](char const *idx_map_)
Definition: ccsd.cxx:115
def rank(self)
Definition: core.pyx:312
Tensor * ai
Definition: ccsd.cxx:156
double divide(double a, double b)
Definition: ccsd.cxx:13
void fill_rand()
Definition: ccsd.cxx:93
Definition: common.h:37
Tensor * aa
Definition: ccsd.cxx:21
an instance of the CTF library (world) on a MPI communicator
Definition: world.h:19
int main(int argc, char **argv)
Definition: ccsd.cxx:333
Tensor * ijab
Definition: ccsd.cxx:32
string
Definition: core.pyx:456
dtype norm2()
computes the frobenius norm of the tensor (needs sqrt()!)
Definition: tensor.h:811
Tensor * ab
Definition: ccsd.cxx:23
int * lens
Tensor * ij
Definition: ccsd.cxx:26
Amplitudes(int no, int nv, World &dw_)
Definition: ccsd.cxx:160
Tensor * abci
Definition: ccsd.cxx:28
Tensor * abij
Definition: ccsd.cxx:157
epoch during which to measure timers
Definition: timer.h:69
Tensor * ii
Definition: ccsd.cxx:22
~Amplitudes()
Definition: ccsd.cxx:169
Integrals(int no, int nv, World &dw_)
Definition: ccsd.cxx:37
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
~Integrals()
Definition: ccsd.cxx:73
char * getCmdOption(char **begin, char **end, const std::string &option)
Definition: ccsd.cxx:322
Definition: apsp.cxx:17
an instance of a tensor within a CTF world
Definition: tensor.h:74
Tensor * aibc
Definition: ccsd.cxx:29
World * dw
Definition: ccsd.cxx:158
CTF::Vector CTF_Vector
Definition: back_comp.h:40
Tensor * abcd
Definition: ccsd.cxx:27
a tensor with an index map associated with it (necessary for overloaded operators) ...
Definition: idx_tensor.h:15
Definition: common.h:37
Tensor * aijk
Definition: ccsd.cxx:33
World * dw
Definition: ccsd.cxx:20
void ccsd(Integrals &V, Amplitudes &T, int sched_nparts=0)
Definition: ccsd.cxx:200
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
Tensor * ia
Definition: ccsd.cxx:25
Definition: common.h:37
MPI_Comm comm
set of processors making up this world
Definition: world.h:22
def np(self)
Definition: core.pyx:315
CTF_int::tensor * parent
Definition: idx_tensor.h:17