Cyclops Tensor Framework parallel arithmetic on multidimensional arrays
sparse_mp3.cxx
Go to the documentation of this file.
1
8 #include <ctf.hpp>
9 #include <float.h>
10 using namespace CTF;
11
12 struct dpair {
13  double a, b;
14  dpair(){ a=0.0; b=0.0; }
15  dpair(double a_, double b_){ a=a_, b=b_;}
16  dpair operator+(dpair const & p) const { return dpair(a+p.a, b+p.b); }
17 };
18
19 namespace CTF {
20  template <>
21  inline void Set<dpair>::print(char const * p, FILE * fp) const {
22  fprintf(fp,"(a=%lf b=%lf)",((dpair*)p)->a, ((dpair*)p)->b);
23  }
24 }
25
27  Tensor<> & Ei,
28  Tensor<> & T,
29  bool sparse_T){
30  if (!sparse_T){
31  Tensor<> D(4,T.lens,*T.wrld);
32  D["abij"] += Ei["i"];
33  D["abij"] += Ei["j"];
34  D["abij"] -= Ea["a"];
35  D["abij"] -= Ea["b"];
36
37  Transform<> div([](double & b){ b=1./b; });
38  div(D["abij"]);
39  T["abij"] = T["abij"]*D["abij"];
40  } else {
41  Tensor<dpair> TD(4,sparse_T,T.lens,*T.wrld,Monoid<dpair,false>(dpair(0.0,0.0)));
42  TD["abij"] = Function<double,dpair>(
43  [](double d){
44  return dpair(d,0.0);
45  })(T["abij"]);
47  [](double d, dpair & p){
48  return p.b += d;
49  });
53  [](double d, dpair & p){
54  return p.b -= d;
55  });
56  bsub(Ea["a"],TD["abij"]);
57  bsub(Ea["b"],TD["abij"]);
58  T["abij"] = Function<dpair,double>(
59  [](dpair p){
60  return p.a/p.b;
61  })(TD["abij"]);
62
63  }
64
65 }
66
67 double mp3(Tensor<> & Ea,
68  Tensor<> & Ei,
69  Tensor<> & Fab,
70  Tensor<> & Fij,
71  Tensor<> & Vabij,
72  Tensor<> & Vijab,
73  Tensor<> & Vabcd,
74  Tensor<> & Vijkl,
75  Tensor<> & Vaibj,
76  bool sparse_T){
77  Tensor<> T(4,sparse_T,Vabij.lens,*Vabij.wrld);
78  T["abij"] = Vabij["abij"];
79
80  divide_EaEi(Ea, Ei, T, sparse_T);
81
82  Tensor<> Z(4,Vabij.lens,*Vabij.wrld);
83  Z["abij"] = Vijab["ijab"];
84  Z["abij"] += Fab["af"]*T["fbij"];
85  Z["abij"] -= Fij["ni"]*T["abnj"];
86  Z["abij"] += 0.5*Vabcd["abef"]*T["efij"];
87  Z["abij"] += 0.5*Vijkl["mnij"]*T["abmn"];
88  Z["abij"] += Vaibj["amei"]*T["ebmj"];
89
90  divide_EaEi(Ea, Ei, Z, 0);
91
92  double MP3_energy = Z["abij"]*Vabij["abij"];
93  return MP3_energy;
94 }
95
96 int sparse_mp3(int nv, int no, World & dw, double sp=.8, bool test=1, int niter=0, bool bnd=1, bool bns=1, bool sparse_T=1){
97  int vvvv[] = {nv,nv,nv,nv};
98  int vovo[] = {nv,no,nv,no};
99  int vvoo[] = {nv,nv,no,no};
100  int oovv[] = {no,no,nv,nv};
101  int oooo[] = {no,no,no,no};
102
103  srand48(dw.rank);
104
105  Vector<> Ea(nv,dw);
106  Vector<> Ei(no,dw);
107
108  Ea.fill_random(1.0*nv*nv,2.0*nv*nv);
109  Ei.fill_random(-2.0*no*no,-1.0*no*no);
110
111  Matrix<> Fab(nv,nv,dw);
112  Matrix<> Fij(no,no,dw);
113
114  Fab.fill_random(-1.0,1.0);
115  Fij.fill_random(-1.0,1.0);
116
117  Tensor<> Vabij(4,vvoo,dw);
118  Tensor<> Vijab(4,oovv,dw);
119  Tensor<> Vabcd(4,vvvv,dw);
120  Tensor<> Vijkl(4,oooo,dw);
121  Tensor<> Vaibj(4,vovo,dw);
122
123  Vabij.fill_random(-1.0,1.0);
124  Vijab.fill_random(-1.0,1.0);
125  Vabcd.fill_random(-1.0,1.0);
126  Vijkl.fill_random(-1.0,1.0);
127  Vaibj.fill_random(-1.0,1.0);
128
129  Transform<> fltr([=](double & d){ if (fabs(d)<sp) d=0.0; });
130  fltr(Vabij["abij"]);
131  fltr(Vijab["ijab"]);
132  fltr(Vabcd["abcd"]);
133  fltr(Vijkl["ijkl"]);
134  fltr(Vaibj["aibj"]);
135
136  double dense_energy, sparse_energy;
137
138  if (test){
139  dense_energy = mp3(Ea, Ei, Fab, Fij, Vabij, Vijab, Vabcd, Vijkl, Vaibj, 0);
140 #ifndef TEST_SUITE
141  if (dw.rank == 0)
142  printf("Calculated MP3 energy %lf with dense integral tensors.\n",dense_energy);
143 #endif
144  } else
145  dense_energy = 0.0;
146
147 #ifndef TEST_SUITE
148  double min_time = DBL_MAX;
149  double max_time = 0.0;
150  double tot_time = 0.0;
151  double times[niter];
152  if (bnd){
153  if (dw.rank == 0){
154  printf("Starting %d benchmarking iterations of dense MP3...\n", niter);
155  }
156  Timer_epoch dmp3("dense MP3");
157  dmp3.begin();
158  for (int i=0; i<niter; i++){
159  double start_time = MPI_Wtime();
160  mp3(Ea, Ei, Fab, Fij, Vabij, Vijab, Vabcd, Vijkl, Vaibj, 0);
161  double end_time = MPI_Wtime();
162  double iter_time = end_time-start_time;
163  times[i] = iter_time;
164  tot_time += iter_time;
165  if (iter_time < min_time) min_time = iter_time;
166  if (iter_time > max_time) max_time = iter_time;
167  }
168  dmp3.end();
169
170  if (dw.rank == 0){
171  printf("Completed %d benchmarking iterations of dense MP3 (no=%d nv=%d sp=%lf).\n", niter, no, nv, sp);
172  printf("All iterations times: ");
173  for (int i=0; i<niter; i++){
174  printf("%lf ", times[i]);
175  }
176  printf("\n");
177  std::sort(times,times+niter);
178  printf("Dense MP3 (no=%d nv=%d sp=%lf p=%d) Min time = %lf, Avg time = %lf, Med time = %lf, Max time = %lf\n",no,nv,sp,dw.np,min_time,tot_time/niter, times[niter/2], max_time);
179  }
180  }
181 #endif
182
183  Vabcd.sparsify();
184  Vabij.sparsify();
185  Vabcd.sparsify();
186  Vijkl.sparsify();
187  Vaibj.sparsify();
188
189  if (test)
190  sparse_energy = mp3(Ea, Ei, Fab, Fij, Vabij, Vijab, Vabcd, Vijkl, Vaibj, sparse_T);
191  else
192  sparse_energy = 0.0;
193
194  bool pass;
195  if (test){
196  pass = fabs((dense_energy-sparse_energy)/dense_energy)<1.E-6;
197
198  if (Ea.wrld->rank == 0){
199  if (!sparse_T){
200  if (pass)
201  printf("{ third-order Moller-Plesset petrubation theory (MP3) using sparse*dense } passed \n");
202  else
203  printf("{ third-order Moller-Plesset petrubation theory (MP3) using sparse*dense } failed \n");
204  } else {
205  if (pass)
206  printf("{ third-order Moller-Plesset petrubation theory (MP3) using sparse*sparse } passed \n");
207  else
208  printf("{ third-order Moller-Plesset petrubation theory (MP3) using sparse*sparse } failed \n");
209  }
210  }
211 #ifndef TEST_SUITE
212  if (dw.rank == 0)
213  printf("Calcluated MP3 energy %lf with sparse integral tensors.\n",sparse_energy);
214 #endif
215  } else pass = 1;
216
217 #ifndef TEST_SUITE
218  if (bns){
219  if (dw.rank == 0){
220  printf("Starting %d benchmarking iterations of sparse MP3...\n", niter);
221  }
222  min_time = DBL_MAX;
223  max_time = 0.0;
224  tot_time = 0.0;
225  Timer_epoch smp3("sparse MP3");
226  smp3.begin();
227  for (int i=0; i<niter; i++){
228  double start_time = MPI_Wtime();
229  mp3(Ea, Ei, Fab, Fij, Vabij, Vijab, Vabcd, Vijkl, Vaibj, sparse_T);
230  double end_time = MPI_Wtime();
231  double iter_time = end_time-start_time;
232 #ifdef TUNE
234 #endif
235  times[i] = iter_time;
236  tot_time += iter_time;
237  if (iter_time < min_time) min_time = iter_time;
238  if (iter_time > max_time) max_time = iter_time;
239  }
240  smp3.end();
241
242  if (dw.rank == 0){
243  printf("Completed %d benchmarking iterations of sparse MP3 (no=%d nv=%d).\n", niter, no, nv);
244  printf("All iterations times: ");
245  for (int i=0; i<niter; i++){
246  printf("%lf ", times[i]);
247  }
248  printf("\n");
249  std::sort(times,times+niter);
250  printf("Sparse MP3 (no=%d nv=%d sp=%lf p=%d spT=%d) Min time = %lf, Avg time = %lf, Med time = %lf, Max time = %lf\n",no,nv,sp,dw.np,sparse_T,min_time,tot_time/niter, times[niter/2], max_time);
251  }
252  }
253 #endif
254  return pass;
255 }
256
257
258 #ifndef TEST_SUITE
259 char* getCmdOption(char ** begin,
260  char ** end,
261  const std::string & option){
262  char ** itr = std::find(begin, end, option);
263  if (itr != end && ++itr != end){
264  return *itr;
265  }
266  return 0;
267 }
268
269
270 int main(int argc, char ** argv){
271  int rank, np, nv, no, pass, niter, bnd, bns, test;
272  bool sparse_T;
273  double sp;
274  int const in_num = argc;
275  char ** input_str = argv;
276
277  MPI_Init(&argc, &argv);
278  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
279  MPI_Comm_size(MPI_COMM_WORLD, &np);
280
281  if (getCmdOption(input_str, input_str+in_num, "-nv")){
282  nv = atoi(getCmdOption(input_str, input_str+in_num, "-nv"));
283  if (nv < 0) nv = 7;
284  } else nv = 7;
285
286  if (getCmdOption(input_str, input_str+in_num, "-no")){
287  no = atoi(getCmdOption(input_str, input_str+in_num, "-no"));
288  if (no < 0) no = 7;
289  } else no = 7;
290
291  if (getCmdOption(input_str, input_str+in_num, "-sp")){
292  sp = atof(getCmdOption(input_str, input_str+in_num, "-sp"));
293  if (sp < 0.0 || sp > 1.0) sp = .8;
294  } else sp = .8;
295
296  if (getCmdOption(input_str, input_str+in_num, "-niter")){
297  niter = atof(getCmdOption(input_str, input_str+in_num, "-niter"));
298  if (niter < 0) niter = 10;
299  } else niter = 10;
300  if (getCmdOption(input_str, input_str+in_num, "-sparse_T")){
301  sparse_T = (bool)atoi(getCmdOption(input_str, input_str+in_num, "-sparse_T"));
302  } else sparse_T = 1;
303
304  if (getCmdOption(input_str, input_str+in_num, "-bnd")){
305  bnd = atoi(getCmdOption(input_str, input_str+in_num, "-bnd"));
306  if (bnd != 0 && bnd != 1) bnd = 0;
307  } else bnd = 0;
308
309  if (getCmdOption(input_str, input_str+in_num, "-bns")){
310  bns = atoi(getCmdOption(input_str, input_str+in_num, "-bns"));
311  if (bns != 0 && bns != 1) bns = 0;
312  } else bns = 0;
313
314  if (getCmdOption(input_str, input_str+in_num, "-test")){
315  test = atoi(getCmdOption(input_str, input_str+in_num, "-test"));
316  if (test != 0 && test != 1) test = 1;
317  } else test = 1;
318
319  if (rank == 0){
320  printf("Running sparse (%lf zeros) third-order Moller-Plesset petrubation theory (MP3) method on %d virtual and %d occupied orbitals and T sparsity turned ot %d\n",sp,nv,no,sparse_T);
321  }
322
323  {
324  World dw;
325  pass = sparse_mp3(nv, no, dw, sp, test, niter, bnd, bns, sparse_T);
326  assert(pass);
327  }
328
329  MPI_Finalize();
330  return 0;
331 }
337 #endif
void update_all_models(MPI_Comm comm)
Definition: model.cxx:15
CTF_int::CommData cdt
communicator data for MPI comm defining this world
Definition: world.h:32
Matrix class which encapsulates a 2D tensor.
Definition: matrix.h:18
def rank(self)
Definition: core.pyx:312
int sparse_mp3(int nv, int no, World &dw, double sp=.8, bool test=1, int niter=0, bool bnd=1, bool bns=1, bool sparse_T=1)
Definition: sparse_mp3.cxx:96
double a
Definition: sparse_mp3.cxx:13
Vector class which encapsulates a 1D tensor.
Definition: vector.h:14
dpair operator+(dpair const &p) const
Definition: sparse_mp3.cxx:16
an instance of the CTF library (world) on a MPI communicator
Definition: world.h:19
string
Definition: core.pyx:456
int main(int argc, char **argv)
Definition: sparse_mp3.cxx:270
void divide_EaEi(Tensor<> &Ea, Tensor<> &Ei, Tensor<> &T, bool sparse_T)
Definition: sparse_mp3.cxx:26
CTF::World * wrld
distributed processor context on which tensor is defined
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
int rank
rank of local processor
Definition: world.h:24
int * lens
void print(char const *a, FILE *fp=stdout) const
prints the value
Definition: set.h:386
void sparsify()
reduce tensor to sparse format, storing only nonzero data, or data above a specified threshold...
Definition: tensor.cxx:449
epoch during which to measure timers
Definition: timer.h:69
dpair(double a_, double b_)
Definition: sparse_mp3.cxx:15
MPI_Comm cm
Definition: common.h:129
double b
Definition: sparse_mp3.cxx:13
char * getCmdOption(char **begin, char **end, const std::string &option)
Definition: sparse_mp3.cxx:259
Definition: apsp.cxx:17
A Monoid is a Set equipped with a binary addition operator &#39;+&#39; or a custom function addition must hav...
Definition: monoid.h:69
an instance of a tensor within a CTF world
Definition: tensor.h:74
double mp3(Tensor<> &Ea, Tensor<> &Ei, Tensor<> &Fab, Tensor<> &Fij, Tensor<> &Vabij, Tensor<> &Vijab, Tensor<> &Vabcd, Tensor<> &Vijkl, Tensor<> &Vaibj, bool sparse_T)
Definition: sparse_mp3.cxx:67
int np
number of processors
Definition: world.h:26
def np(self)
Definition: core.pyx:315