1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
8 //
9 // File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
10 //////////////////////////////////////////////////////////////////////////////////////
11 // -*- C++ -*-
12 /** @file restart.cpp
13 * @brief developing restart IO
14 */
15
16 #include <Configuration.h>
17 #include "Message/CommOperators.h"
18 #include "Particle/MCWalkerConfiguration.h"
19 #include "Particle/HDFWalkerOutput.h"
20 #include "hdf/HDFVersion.h"
21 #include "Particle/HDFWalkerInput_0_4.h"
22 #include "OhmmsApp/RandomNumberControl.h"
23 #include "Sandbox/input.hpp"
24 #include "Sandbox/pseudo.hpp"
25 #include "Utilities/FairDivide.h"
26 #include "Utilities/Timer.h"
27 #include "Sandbox/common.hpp"
28 #include <getopt.h>
29 #include "mpi/collectives.h"
30 #include "ParticleBase/ParticleAttribOps.h"
31
32 using namespace std;
33 using namespace qmcplusplus;
34
setWalkerOffsets(MCWalkerConfiguration & W,Communicate * myComm)35 void setWalkerOffsets(MCWalkerConfiguration& W, Communicate* myComm)
36 {
37 std::vector<int> nw(myComm->size(),0),nwoff(myComm->size()+1,0);
38 nw[myComm->rank()]=W.getActiveWalkers();
39 myComm->allreduce(nw);
40 for(int ip=0; ip<myComm->size(); ip++)
41 nwoff[ip+1]=nwoff[ip]+nw[ip];
42 W.setGlobalNumWalkers(nwoff[myComm->size()]);
43 W.setWalkerOffsets(nwoff);
44 }
45
main(int argc,char ** argv)46 int main(int argc, char** argv)
47 {
48 #ifdef HAVE_MPI
49 mpi3::environment env(0, NULL);
50 OHMMS::Controller->initialize(env);
51 #endif
52 Communicate* myComm = OHMMS::Controller;
53 myComm->setName("restart");
54 myComm->barrier();
55
56 typedef QMCTraits::RealType RealType;
57 typedef ParticleSet::ParticlePos_t ParticlePos_t;
58 typedef ParticleSet::ParticleLayout_t LatticeType;
59 typedef ParticleSet::TensorType TensorType;
60 typedef ParticleSet::PosType PosType;
61 typedef RandomGenerator_t::uint_type uint_type;
62 typedef MCWalkerConfiguration::Walker_t Walker_t;
63
64 //use the global generator
65
66 int na=4;
67 int nb=4;
68 int nc=1;
69 int nsteps=100;
70 int iseed=11;
71 int AverageWalkersPerNode=0;
72 int nwtot;
73 std::vector<int> wPerNode;
74 RealType Rmax(1.7);
75
76 const int NumThreads=omp_get_max_threads();
77
78 char *g_opt_arg;
79 int opt;
80 while((opt = getopt(argc, argv, "hg:i:r:")) != -1)
81 {
82 switch(opt)
83 {
84 case 'h':
85 printf("[-g \"n0 n1 n2\"]\n");
86 return 1;
87 case 'g': //tiling1 tiling2 tiling3
88 sscanf(optarg,"%d %d %d",&na,&nb,&nc);
89 break;
90 case 'i': //number of MC steps
91 nsteps=atoi(optarg);
92 break;
93 case 's'://random seed
94 iseed=atoi(optarg);
95 break;
96 case 'w'://the number of walkers
97 AverageWalkersPerNode=atoi(optarg);
98 break;
99 case 'r'://rmax
100 Rmax=atof(optarg);
101 break;
102 }
103 }
104
105 // set the number of walkers equal to the threads.
106 if(!AverageWalkersPerNode) AverageWalkersPerNode=NumThreads;
107 //set nwtot, to be random
108 nwtot=std::abs(AverageWalkersPerNode+myComm->rank()%5-2);
109 FairDivideLow(nwtot,NumThreads,wPerNode);
110
111 //Random.init(0,1,iseed);
112 Tensor<int,3> tmat(na,0,0,0,nb,0,0,0,nc);
113
114 //turn off output
115 if(myComm->rank())
116 {
117 outputManager.shutOff();
118 }
119
120 int nptcl=0;
121 double t0=0.0,t1=0.0;
122
123 RandomNumberControl::make_seeds();
124 std::vector<RandomGenerator_t> myRNG(NumThreads);
125 std::vector<uint_type> mt(Random.state_size(),0);
126 std::vector<MCWalkerConfiguration> elecs(NumThreads);
127
128 ParticleSet ions;
129 OHMMS_PRECISION scale=1.0;
130 tile_cell(ions,tmat,scale);
131
132 #pragma omp parallel reduction(+:t0)
133 {
134 int ip=omp_get_thread_num();
135
136 MCWalkerConfiguration& els=elecs[ip];
137
138 //create generator within the thread
139 myRNG[ip]=*RandomNumberControl::Children[ip];
140 RandomGenerator_t& random_th=myRNG[ip];
141
142 const int nions=ions.getTotalNum();
143 const int nels=count_electrons(ions);
144 const int nels3=3*nels;
145
146 #pragma omp master
147 nptcl=nels;
148
149 {//create up/down electrons
150 els.Lattice.BoxBConds = 1;
151 els.Lattice = ions.Lattice;
152 vector<int> ud(2); ud[0]=nels/2; ud[1]=nels-ud[0];
153 els.create(ud);
154 els.R.InUnit = PosUnit::Lattice;
155 random_th.generate_uniform(&els.R[0][0],nels3);
156 els.convert2Cart(els.R); // convert to Cartiesian
157 els.setCoordinates(els.R);
158 }
159
160 if(!ip) elecs[0].createWalkers(nwtot);
161 #pragma omp barrier
162
163 for(MCWalkerConfiguration::iterator wi=elecs[0].begin()+wPerNode[ip]; wi!=elecs[0].begin()+wPerNode[ip+1]; wi++)
164 els.saveWalker(**wi);
165
166 // save random seeds and electron configurations.
167 *RandomNumberControl::Children[ip]=myRNG[ip];
168 //MCWalkerConfiguration els_save(els);
169
170 } //end of omp parallel
171 Random.save(mt);
172
173 setWalkerOffsets(elecs[0], myComm);
174
175 //storage variables for timers
176 double h5write = 0.0, h5read = 0.0; //random seed R/W speeds
177 double walkerWrite = 0.0, walkerRead =0.0; //walker R/W speeds
178 Timer h5clock; //timer for the program
179
180 // dump random seeds
181 myComm->barrier();
182 h5clock.restart(); //start timer
183 RandomNumberControl::write("restart",myComm);
184 myComm->barrier();
185 h5write += h5clock.elapsed(); //store timer
186
187 // flush random seeds to zero
188 #pragma omp parallel
189 {
190 int ip=omp_get_thread_num();
191 RandomGenerator_t& random_th=*RandomNumberControl::Children[ip];
192 std::vector<uint_type> vt(random_th.state_size(),0);
193 random_th.load(vt);
194 }
195 std::vector<uint_type> mt_temp(Random.state_size(),0);
196 Random.load(mt_temp);
197
198 // load random seeds
199 myComm->barrier();
200 h5clock.restart(); //start timer
201 RandomNumberControl::read("restart",myComm);
202 myComm->barrier();
203 h5read += h5clock.elapsed(); //store timer
204
205 // validate random seeds
206 int mismatch_count=0;
207 #pragma omp parallel reduction(+:mismatch_count)
208 {
209 int ip=omp_get_thread_num();
210 RandomGenerator_t& random_th=myRNG[ip];
211 std::vector<uint_type> vt_orig(random_th.state_size());
212 std::vector<uint_type> vt_load(random_th.state_size());
213 random_th.save(vt_orig);
214 RandomNumberControl::Children[ip]->save(vt_load);
215 for(int i=0; i<random_th.state_size(); i++)
216 if(vt_orig[i]!=vt_load[i]) mismatch_count++;
217 }
218 Random.save(mt_temp);
219 for(int i=0; i<Random.state_size(); i++)
220 if(mt_temp[i]!=mt[i]) mismatch_count++;
221
222 myComm->allreduce(mismatch_count);
223
224 if(!myComm->rank())
225 {
226 if(mismatch_count!=0)
227 std::cout << "Fail: random seeds mismatch between write and read!\n"
228 << " state_size= " << myRNG[0].state_size() << " mismatch_cout=" << mismatch_count << std::endl;
229 else
230 std::cout << "Pass: random seeds match exactly between write and read!\n";
231 }
232
233 // dump electron coordinates.
234 HDFWalkerOutput wOut(elecs[0],"restart",myComm);
235 myComm->barrier();
236 h5clock.restart(); //start timer
237 wOut.dump(elecs[0],1);
238 myComm->barrier();
239 walkerWrite += h5clock.elapsed(); //store timer
240 if(!myComm->rank()) std::cout << "Walkers are dumped!\n";
241
242 // save walkers before destroying them
243 std::vector<Walker_t> saved_walkers;
244 for(int wi=0; wi<elecs[0].getActiveWalkers(); wi++)
245 saved_walkers.push_back(*elecs[0][wi]);
246 elecs[0].destroyWalkers(elecs[0].begin(),elecs[0].end());
247
248 // load walkers
249 const char *restart_input = \
250 "<tmp> \
251 <mcwalkerset fileroot=\"restart\" node=\"-1\" version=\"3 0\" collected=\"yes\"/> \
252 </tmp> \
253 ";
254
255 Libxml2Document doc;
256 bool okay = doc.parseFromString(restart_input);
257 xmlNodePtr root = doc.getRoot();
258 xmlNodePtr restart_leaf = xmlFirstElementChild(root);
259
260 HDFVersion in_version(0,4);
261 HDFWalkerInput_0_4 wIn(elecs[0],myComm,in_version);
262 myComm->barrier();
263 h5clock.restart(); //start timer
264 wIn.put(restart_leaf);
265 myComm->barrier();
266 walkerRead += h5clock.elapsed(); //store time spent
267
268 if(saved_walkers.size()!=elecs[0].getActiveWalkers())
269 std::cout << "Fail: Rank " << myComm->rank() << " had " << saved_walkers.size()
270 << " walkers but loaded only " << elecs[0].getActiveWalkers() << " walkers from file!" << std::endl;
271
272 mismatch_count=0;
273 for(int wi=0; wi<saved_walkers.size(); wi++)
274 {
275 saved_walkers[wi].R=saved_walkers[wi].R-elecs[0][wi]->R;
276 if(Dot(saved_walkers[wi].R,saved_walkers[wi].R)>std::numeric_limits<RealType>::epsilon()) mismatch_count++;
277 }
278 myComm->allreduce(mismatch_count);
279
280 if(!myComm->rank())
281 {
282 if(mismatch_count!=0)
283 std::cout << "Fail: electron coordinates mismatch between write and read!\n"
284 << " mismatch_cout=" << mismatch_count << std::endl;
285 else
286 std::cout << "Pass: electron coordinates match exactly between write and read!\n";
287 }
288
289 //print out hdf5 R/W times
290 TinyVector<double,4> timers(h5read, h5write, walkerRead, walkerWrite);
291 mpi::reduce(*myComm, timers);
292 h5read = timers[0]/myComm->size();
293 h5write = timers[1]/myComm->size();
294 walkerRead = timers[2]/myComm->size();
295 walkerWrite = timers[3]/myComm->size();
296 if(myComm->rank() == 0)
297 {
298 cout << "\nTotal time of writing random seeds to HDF5 file: " << setprecision(6) << h5write << "\n";
299 cout << "\nTotal time of reading random seeds in HDF5 file: " << setprecision(6) << h5read << "\n";
300 cout << "\nTotal time of writing walkers to HDF5 file: " << setprecision(6) << walkerWrite << "\n";
301 cout << "\nTotal time of reading walkers in HDF5 file: " << setprecision(6) << walkerRead << "\n";
302 }
303
304 if(myComm->size()>1)
305 {
306 Communicate* subComm = new Communicate(*myComm, 2);
307 subComm->setName("restart2");
308
309 if(subComm->getGroupID() == 0)
310 {
311 elecs[0].destroyWalkers(elecs[0].begin(),elecs[0].end());
312 HDFWalkerInput_0_4 subwIn(elecs[0],subComm,in_version);
313 subwIn.put(restart_leaf);
314 subComm->barrier();
315 if(!subComm->rank()) std::cout << "Walkers are loaded again by the subgroup!\n";
316 setWalkerOffsets(elecs[0], subComm);
317 HDFWalkerOutput subwOut(elecs[0],"XXXX",subComm);
318 subwOut.dump(elecs[0],1);
319 if(!subComm->rank()) std::cout << "Walkers are dumped again by the subgroup!\n";
320 }
321 delete subComm;
322 }
323
324 OHMMS::Controller->finalize();
325
326 return 0;
327 }
328