QDP++
qdp_hdf5.cc
Go to the documentation of this file.
1//
2// QDP data parallel interface
7
8#include "qdp.h"
9
10#ifdef QDP_USE_HDF5
11
12// optional hdf5 code goes here
13namespace QDP {
14 using std::string;
15
16 // HDF5 classes
17 //--------------------------------------------------------------------------------
18 //--------------------------------------------------------------------------------
19 // Base class
20 //--------------------------------------------------------------------------------
21 //--------------------------------------------------------------------------------
22 //constructor, never called by user:
23
24 HDF5::HDF5(const long int& stripesizee, const long int& maxalignn) :
25 error_stack(H5E_DEFAULT), par_init(false), file_comm(H5P_DEFAULT),
26 file_id(-1), current_group(-1),
27 stripesize(stripesizee),
28 maxalign(maxalignn)
29 {
30 //no profiling
31 profile=false;
32 isprefetched=false;
33
34 //do not do cleanup after finishing:
35 herr_t err=H5dont_atexit();
36
37 //set error handler to custom error handler:
38 H5E_auto2_t oldfunc;
39 void *client_data;
40
41 //get old error handler
42 H5Eget_auto(error_stack, &oldfunc, &client_data);
43
44 //get new error handler
45 //H5Eset_auto(error_stack,errorHandler,NULL);
46 //H5Eset_auto(error_stack,NULL,NULL);
47 }
48
49 //checks whether an HDF5 file already exists:
50 bool HDF5::check_exists(const std::string& filename){
51 //on node 0, test if file exists:
52 bool exists=false;
53 if(Layout::nodeNumber()==0){
54 std::ifstream input(filename.c_str(),std::ios_base::binary);
55 if(input.good()){
56 exists=true;
57 }
58 input.close();
59
60 if(exists){
61 htri_t ex=H5Fis_hdf5(filename.c_str());
62 if(ex<=0) exists=false;
63 }
64 }
65 QDPInternal::broadcast(exists);
66
67 return exists;
68 }
69
70 //error handler:
71 hid_t HDF5::errorHandler(hid_t errstack, void* unused){
72 QDPIO::cout << "Some error occured, but we do not care at the moment!" << std::endl;
73 return errstack;
74 }
75
76 //lookup routines:
77 std::string HDF5::pwd()const{
78 return getNameById(current_group);
79 }
80
81 std::string HDF5::parentDir()const{
82 size_t pos;
83 std::string dir=pwd();
84 if((pos=dir.find_last_of("/"))==0) return std::string("/");
85 return dir.substr(0,pos);
86 }
87
88 std::string HDF5::getNameById(hid_t id)const{
89 ssize_t size=H5Iget_name(id,NULL,0);
90 char* val=new char[size+1];
91 H5Iget_name(id,val,size+1);
92 std::string name(val);
93 delete [] val;
94 return name;
95 }
96
97 void HDF5::tokenize(const ::std::string& str, ::std::vector< ::std::string >& tokens, const ::std::string& delimiters){
98 // Skip delimiters at beginning.
99 ::std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
100 // Find first "non-delimiter".
101 ::std::string::size_type pos = str.find_first_of(delimiters, lastPos);
102
103 while (::std::string::npos != pos || ::std::string::npos != lastPos){
104 // Found a token, add it to the vector.
105 tokens.push_back(str.substr(lastPos, pos - lastPos));
106 // Skip delimiters. Note the "not_of"
107 lastPos = str.find_first_not_of(delimiters, pos);
108 // Find next "non-delimiter"
109 pos = str.find_first_of(delimiters, lastPos);
110 }
111 }
112
113 std::vector<std::string> HDF5::splitPathname(const std::string& name){
114 //first split the string according to the file separators, such that the linking is correct:
115 ::std::vector<std::string> dirlist;
116 tokenize(name, dirlist, "/");
117
118 //remove zero strings, they can occur when name contained expressions such as //:
119 ::std::vector<std::string>::iterator it;
120 for(it=dirlist.begin(); it!=dirlist.end(); it++){
121 if(*it=="") dirlist.erase(it);
122 }
123
124 return dirlist;
125 }
126
127 //check if object exists:
128 bool HDF5::objectExists(const ::std::string& name){
129 return objectExists(file_id,name);
130 }
131
132 bool HDF5::objectExists(hid_t loc_id, const std::string& name){
133 std::vector<std::string> dirlist=splitPathname(name);
134
135 hid_t start_group;
136 if(name.find_first_of("/")==0){
137 start_group=file_id;
138 }
139 else start_group=loc_id;
140
141 //iterate through the tree and check whether everything on the way exists:
142 htri_t exists=H5Lexists(start_group,dirlist[0].c_str(),H5P_DEFAULT);
143 if(exists!=1) return false;
144 exists=H5Oexists_by_name(start_group,dirlist[0].c_str(),H5P_DEFAULT);
145 if(exists!=1) return false;
146 std::string tmpstring=dirlist[0];
147 for(unsigned int i=1; i<dirlist.size(); i++){
148 tmpstring+="/"+dirlist[i];
149 exists=H5Lexists(start_group,tmpstring.c_str(),H5P_DEFAULT);
150 if(exists!=1) return false;
151 exists=H5Oexists_by_name(start_group,tmpstring.c_str(),H5P_DEFAULT);
152 if(exists!=1) return false;
153 }
154 return true;
155 }
156
157 ::std::string HDF5::objectType(const ::std::string& name){
158 return objectType(file_id, name);
159 }
160
161 ::std::string HDF5::objectType(hid_t loc_id, const ::std::string& name){
162 std::string result="Null";
163 if(objectExists(loc_id, name)){
164 H5O_info_t objinfo;
165 herr_t errhandle=H5Oget_info_by_name(loc_id,name.c_str(),&objinfo,H5P_DEFAULT);
166 if(errhandle<0){
167 ::std::cerr << "HDF5::objectType: error, something went wrong with looking up " << name << "!" << std::endl;
168 }
169 else{
170 switch(objinfo.type){
171 case H5O_TYPE_GROUP:
172 result="Group";
173 break;
174
175 case H5O_TYPE_DATASET:
176 result="Dataset";
177 break;
178
179 case H5O_TYPE_NAMED_DATATYPE:
180 result="NamedDatatype";
181 break;
182
183 case H5O_TYPE_NTYPES:
184 result="VariousTypes";
185 break;
186
187 default:
188 result="Unknown";
189 break;
190 }
191 }
192 }
193 return result;
194 }
195
196
197 //navigation routines:
198 void HDF5::pop(){
199 //only do something if not at root-level:
200 if(pwd().compare("/")!=0){
201 hid_t tmp_group;
202
203 //lookup what current and parent directory is:
204 ::std::string cdir=pwd();
205 ::std::string pdir=parentDir();
206
207 //close group to prevent data loss and then go to new group:
208 H5Gclose(current_group);
209 tmp_group=H5Gopen(file_id,pdir.c_str(),H5P_DEFAULT);
210
211 current_group=tmp_group;
212 }
213 }
214
215 void HDF5::cd(const std::string& dirname){
216 std::string tmpdirname(dirname);
217 if(tmpdirname.compare("..")==0){
218 pop();
219 }
220 else if(tmpdirname.compare(pwd())!=0){
221 hid_t tmp_group;
222
223 //open HDF5 group relative to current_group:
224 tmp_group=H5Gopen(current_group,tmpdirname.c_str(),H5P_DEFAULT);
225
226 if(tmp_group<0){
227 QDPIO::cerr << "HDF5::cd: error, the group " << tmpdirname << " does not exist!" << std::endl;
228 return;
229 }
230 //only close current_group if tmp_group is not further below in the tree or completely outside:
231 //get the names of the new and current dir:
232 std::string ndir=getNameById(tmp_group);
233 std::string cdir=pwd();
234
235 //perform different tests:
236 if(ndir.find(cdir)==::std::string::npos || ndir.compare("/")==0){
237 //ndir is not inside cdir tree or ndir is root directory, so close cdir tree:
238 while(pwd().compare("/")!=0){
239 pop();
240 }
241 }
242 else if(cdir.find(ndir)==0){
243 //ndir is parent of cdir: iterate up from cdir and close all dirs on the way:
244 while(current_group!=tmp_group){
245 pop();
246 }
247 }
248 current_group=tmp_group;
249 }
250 }
251
252 //close-routine
253 int HDF5::close(){
254 if(file_id>0){
255 herr_t err=1;
256
257 unsigned int types=H5F_OBJ_DATASET | H5F_OBJ_GROUP |
258 H5F_OBJ_DATATYPE | H5F_OBJ_ATTR;
259
260 //get number of still open objects:
261 ssize_t num_open = H5Fget_obj_count(file_id,types);
262 if (num_open > 0) {
263 std::vector<hid_t> open_object_ids(num_open, 0);
264 H5Fget_obj_ids(file_id, types, num_open, &(open_object_ids.front()) );
265 for(unsigned int i=0; i<num_open; i++){
266 err=H5Oclose(open_object_ids[i]);
267 }
268 }
269 err=H5Fclose(file_id);
270 }
271
272 //reset file_id:
273 file_id=-1;
274
275 return EXIT_SUCCESS;
276 }
277
278 //prefetch mapping for CB->lexicographical:
279 int HDF5::prefetchLatticeCoordinates(){
280 const int mynode = Layout::nodeNumber();
281
282 //measure
283 reordermap.resize(Layout::sitesOnNode());
284
285 unsigned int run=0;
286 for(int site=0; site < Layout::vol(); ++site){
287 multi1d<int> coord = crtesn(site, Layout::lattSize());
288 int node = Layout::nodeNumber(coord);
289
290 if(node==mynode){
291 reordermap[run]=Layout::linearSiteIndex(coord);
292 run++;
293 }
294 }
295 isprefetched=true;
296
297 return EXIT_SUCCESS;
298 }
299
300 //***********************************************************************************************************************************
301 //***********************************************************************************************************************************
302 //DATATYPE HELPERS
303 //***********************************************************************************************************************************
304 //***********************************************************************************************************************************
305 //create complex
306 hid_t HDF5::createComplexType(const unsigned int& float_size){
307 hid_t complex_id = H5Tcreate(H5T_COMPOUND, 2*float_size);
308 hid_t base_type;
309
310 if(float_size==4) base_type=H5Tcopy(H5T_NATIVE_FLOAT);
311 else base_type=H5Tcopy(H5T_NATIVE_DOUBLE);
312
313 H5Tinsert(complex_id, "r", 0, base_type);
314 H5Tinsert(complex_id, "i", float_size, base_type);
315 H5Tclose(base_type);
316
317 return complex_id;
318 }
319
320 //check complex:
321 bool HDF5::checkComplexType(const hid_t& type_id, hid_t& base_type_id){
322 //check if type_id is complex class:
323 if(H5Tget_class(type_id)!=H5T_COMPOUND){
324 return false;
325 }
326
327 if(H5Tget_nmembers(type_id)!=2) return false;
328 base_type_id=H5Tget_member_type(type_id,0);
329 htri_t err=H5Tequal(base_type_id,H5Tget_member_type(type_id,1));
330 if(err<=0){
331 base_type_id=-1;
332 return false;
333 }
334
335 //check if data members are named r and i:
336 bool error=false;
337 char* name=H5Tget_member_name(type_id,0);
338 if( strcmp(name,"r") != 0 ) error=true;
339 free(name);
340 name=H5Tget_member_name(type_id,1);
341 if( strcmp(name,"i") != 0 ) error=true;
342 free(name);
343 if(error){
344 base_type_id=-1;
345 return false;
346 }
347
348 return true;
349 }
350
351 //create colormat
352 hid_t HDF5::createColorMatrixType(const hid_t& complex_id, const unsigned int& rank){
353 if(rank>Nc || rank==0){
354 HDF5_error_exit("HDF::createColorMatrixType: error, the rank should be between 1 and Nc");
355 }
356 hsize_t dims[]={rank,rank};
357 hid_t colmat_id = H5Tarray_create(complex_id,2,dims);
358
359 return colmat_id;
360 }
361
362 //check colormat:
363 bool HDF5::checkColorMatrixType(const hid_t& type_id, const unsigned int& rank, hid_t& base_type_id){
365 if( (H5Tget_class(type_id)!=H5T_ARRAY) || (H5Tget_array_ndims(type_id)!=2) ){
366 return false;
367 }
368
369 //each dimension should have rank rank:
370 hsize_t dims[2];
371 H5Tget_array_dims(type_id,dims);
372 if( (dims[0]!=dims[1]) || (dims[0]!=rank) ){
373 return false;
374 }
375
376 //get base datatype and see whether it is the same for both dimensions:
377 hid_t member_type=H5Tget_super(type_id);
378 if(!checkComplexType(member_type,base_type_id)){
379 return false;
380 }
381 H5Tclose(member_type);
382
383 return true;
384 }
385
386 //create propagator:
387 hid_t HDF5::createPropagatorType(const hid_t& colmat_id, const unsigned int& spinrank){
388 if(spinrank>Ns || spinrank==0){
389 HDF5_error_exit("HDF::createPropagatorType: error, the spinrank should be between 1 and Ns");
390 }
391 hsize_t dims[]={spinrank,spinrank};
392 hid_t prop_id = H5Tarray_create(colmat_id,2,dims);
393
394 return prop_id;
395 }
396
397 //check propagator:
398 bool HDF5::checkDiracPropagatorType(const hid_t& type_id, const unsigned int& spinrank, const unsigned int& colorrank, hid_t& base_type_id){
399 //datatype specification is: [4,4][3,3] {re,im}
400
402 if( (H5Tget_class(type_id)!=H5T_ARRAY) || (H5Tget_array_ndims(type_id)!=2) ){
403 return false;
404 }
405
406 //each dimension should have rank Ns:
407 hsize_t dims[2];
408 H5Tget_array_dims(type_id,dims);
409 if( (dims[0]!=dims[1]) || (dims[0]!=spinrank) ){
410 return false;
411 }
412
413 //check inner type
414 hid_t member_type=H5Tget_super(type_id);
415 if(!checkColorMatrixType(member_type,colorrank,base_type_id)){
416 return false;
417 }
418 H5Tclose(member_type);
419
420 return true;
421 }
422
423
424 //***********************************************************************************************************************************
425 //***********************************************************************************************************************************
426 //READING ATTRIBUTES
427 //***********************************************************************************************************************************
428 //***********************************************************************************************************************************
429 //single datum
430 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, short& datum){
431 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,true);
432 }
433
434 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, unsigned short& datum){
435 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,false);
436 }
437
438 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, int& datum){
439 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,true);
440 }
441
442 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, unsigned int& datum){
443 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,false);
444 }
445
446 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, unsigned long long& datum){
447 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,false);
448 }
449
450 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, float& datum){
451 rdAtt(obj_name,attr_name,datum,H5T_FLOAT,true);
452 }
453
454 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, double& datum){
455 rdAtt(obj_name,attr_name,datum,H5T_FLOAT,true);
456 }
457
458 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, std::string& datum){
459 std::string oname(obj_name), aname(attr_name);
460 bool exists=objectExists(current_group,oname);
461 if(!exists){
462 HDF5_error_exit("HDF5::readAttribute: error, object "+oname+" you try to read attribute from does not exists!");
463 }
464
465 //do sanity checks and get datatype
466 hid_t ex=H5Aexists_by_name(current_group,oname.c_str(),aname.c_str(),H5P_DEFAULT);
467 if(ex!=1){
468 HDF5_error_exit("HDF5::readAttribute: error, the attribute "+aname+" you try to read does not exists!");
469 }
470 hid_t attr_id=H5Aopen_by_name(current_group,oname.c_str(),aname.c_str(),H5P_DEFAULT,H5P_DEFAULT);
471 if(attr_id<0){
472 HDF5_error_exit("HDF5::readAttribute: error, cannot open attribute "+aname+" attached to "+oname+"!");
473 }
474 hid_t type_id=H5Aget_type(attr_id);
475 if(H5Tget_class(type_id)!=H5T_STRING){
476 HDF5_error_exit("HDF5::readAttribute: error, datatype mismatch in attribute "+aname+"!");
477 }
478
479 //memory size
480 hsize_t size=H5Tget_size(type_id);
481 char* datumcpy=new char[size];
482 hid_t nat_type_id=H5Tget_native_type(type_id,H5T_DIR_ASCEND);
483 H5Aread(attr_id,nat_type_id,reinterpret_cast<void*>(datumcpy));
484 datum=std::string(datumcpy);
485 delete [] datumcpy;
486 H5Aclose(attr_id);
487 H5Tclose(nat_type_id);
488 H5Tclose(type_id);
489 }
490
491 //array
492 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, multi1d<short>& datum){
493 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,true);
494 }
495
496 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, multi1d<unsigned short>& datum){
497 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,false);
498 }
499
500 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, multi1d<int>& datum){
501 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,true);
502 }
503
504 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, multi1d<unsigned int>& datum){
505 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,false);
506 }
507
508 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, multi1d<unsigned long long>& datum){
509 rdAtt(obj_name,attr_name,datum,H5T_INTEGER,false);
510 }
511
512 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, multi1d<float>& datum){
513 rdAtt(obj_name,attr_name,datum,H5T_FLOAT,true);
514 }
515
516 void HDF5::readAttribute(const std::string& obj_name, const std::string& attr_name, multi1d<double>& datum){
517 rdAtt(obj_name,attr_name,datum,H5T_FLOAT,true);
518 }
519
520 //***********************************************************************************************************************************
521 //***********************************************************************************************************************************
522 //READING DATASETS
523 //***********************************************************************************************************************************
524 //***********************************************************************************************************************************
525 //single datum
526 void HDF5::read(const std::string& obj_name, short& datum){
527 rd(obj_name,datum,H5T_INTEGER,true);
528 }
529
530 void HDF5::read(const std::string& obj_name, unsigned short& datum){
531 rd(obj_name,datum,H5T_INTEGER,false);
532 }
533
534 void HDF5::read(const std::string& obj_name, int& datum){
535 rd(obj_name,datum,H5T_INTEGER,true);
536 }
537
538 void HDF5::read(const std::string& obj_name, unsigned int& datum){
539 rd(obj_name,datum,H5T_INTEGER,false);
540 }
541
542 void HDF5::read(const std::string& obj_name, unsigned long long& datum){
543 rd(obj_name,datum,H5T_INTEGER,false);
544 }
545
546 void HDF5::read(const std::string& obj_name, float& datum){
547 rd(obj_name,datum,H5T_FLOAT,true);
548 }
549
550 void HDF5::read(const std::string& obj_name, double& datum){
551 rd(obj_name,datum,H5T_FLOAT,true);
552 }
553
554 void HDF5::read(const std::string& dataname, std::string& datum){
555 std::string dname(dataname);
556
557 bool exists=objectExists(current_group,dname);
558 if(!exists){
559 HDF5_error_exit("HDF5::read: error, dataset does not exists!");
560 }
561 H5O_info_t objinfo;
562 herr_t errhandle=H5Oget_info_by_name(current_group,dname.c_str(),&objinfo,H5P_DEFAULT);
563 if(objinfo.type!=H5O_TYPE_DATASET){
564 HDF5_error_exit("HDF5::read: error, "+dname+" exists but it is not a dataset!");
565 }
566 hid_t dset_id=H5Dopen(current_group,dname.c_str(),H5P_DEFAULT);
567 if(dset_id<0){
568 HDF5_error_exit("HDF5::read: error, cannot open dataset!");
569 }
570 hid_t type_id=H5Dget_type(dset_id);
571 if(H5Tget_class(type_id)!=H5T_STRING){
572 HDF5_error_exit("HDF5::read: error, datatype mismatch in dataset "+dataname+"!");
573 }
574
575 //memory size
576 hsize_t size=H5Tget_size(type_id);
577 char* datumcpy=new char[size];
578 hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
579 H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
580 hid_t nat_type_id=H5Tget_native_type(type_id,H5T_DIR_ASCEND);
581 H5Dread(dset_id,nat_type_id,H5S_ALL,H5S_ALL,plist_id,reinterpret_cast<void*>(datumcpy));
582 datum=std::string(datumcpy);
583 delete [] datumcpy;
584 H5Pclose(plist_id);
585 H5Dclose(dset_id);
586 H5Tclose(nat_type_id);
587 H5Tclose(type_id);
588 }
589
590 //array:
591 //1D
592 void HDF5::read(const std::string& obj_name, multi1d<short>& datum){
593 rd(obj_name,datum,H5T_INTEGER,true);
594 }
595
596 void HDF5::read(const std::string& obj_name, multi1d<unsigned short>& datum){
597 rd(obj_name,datum,H5T_INTEGER,false);
598 }
599
600 void HDF5::read(const std::string& obj_name, multi1d<int>& datum){
601 rd(obj_name,datum,H5T_INTEGER,true);
602 }
603
604 void HDF5::read(const std::string& obj_name, multi1d<unsigned int>& datum){
605 rd(obj_name,datum,H5T_INTEGER,false);
606 }
607
608 void HDF5::read(const std::string& obj_name, multi1d<unsigned long long>& datum){
609 rd(obj_name,datum,H5T_INTEGER,false);
610 }
611
612 void HDF5::read(const std::string& obj_name, multi1d<float>& datum){
613 rd(obj_name,datum,H5T_FLOAT,true);
614 }
615
616 void HDF5::read(const std::string& obj_name, multi1d<double>& datum){
617 rd(obj_name,datum,H5T_FLOAT,true);
618 }
619 //2D
620 void HDF5::read(const std::string& obj_name, multi2d<short>& datum){
621 rd(obj_name,datum,H5T_INTEGER,true);
622 }
623
624 void HDF5::read(const std::string& obj_name, multi2d<unsigned short>& datum){
625 rd(obj_name,datum,H5T_INTEGER,false);
626 }
627
628 void HDF5::read(const std::string& obj_name, multi2d<int>& datum){
629 rd(obj_name,datum,H5T_INTEGER,true);
630 }
631
632 void HDF5::read(const std::string& obj_name, multi2d<unsigned int>& datum){
633 rd(obj_name,datum,H5T_INTEGER,false);
634 }
635
636 void HDF5::read(const std::string& obj_name, multi2d<unsigned long long>& datum){
637 rd(obj_name,datum,H5T_INTEGER,false);
638 }
639
640 void HDF5::read(const std::string& obj_name, multi2d<float>& datum){
641 rd(obj_name,datum,H5T_FLOAT,true);
642 }
643
644 void HDF5::read(const std::string& obj_name, multi2d<double>& datum){
645 rd(obj_name,datum,H5T_FLOAT,true);
646 }
647
648 //***********************************************************************************************************************************
649 //***********************************************************************************************************************************
650 //READING Compound types:
651 //***********************************************************************************************************************************
652 //***********************************************************************************************************************************
653 void HDF5::readPrepare(const std::string& name, hid_t& type_id){
654 //determine whether there is a dataset with the specified name:
655 bool exists=objectExists(current_group,name);
656 if(!exists){
657 HDF5_error_exit("HDF5::read: error, dataset does not exists!");
658 }
659 H5O_info_t objinfo;
660 herr_t errhandle=H5Oget_info_by_name(current_group,name.c_str(),&objinfo,H5P_DEFAULT);
661 if(objinfo.type!=H5O_TYPE_DATASET){
662 HDF5_error_exit("HDF5::read: error, "+name+" exists but it is not a dataset!");
663 }
664
665 //open dataset
666 hid_t dset_id=H5Dopen(current_group,name.c_str(),H5P_DEFAULT);
667 if(dset_id<0){
668 HDF5_error_exit("HDF5::read: error, cannot open dataset!");
669 }
670 type_id=H5Dget_type(dset_id);
671 H5Dclose(dset_id);
672 }
673
674 //complex types:
675 //Single element:
676 template<>void HDF5::read< PScalar< PScalar< RComplex<float> > > >(const std::string& dataname, ComplexF& datum){
677 //get type:
678 hid_t type_id;
679 readPrepare(dataname,type_id);
680 if(type_id<0){
681 HDF5_error_exit("HDF5::read: error, cannot open datatype!");
682 }
683 hid_t base_type_id;
684 if(!checkComplexType(type_id,base_type_id)){
685 HDF5_error_exit("HDF5::read: error, datatype mismatch!");
686 }
687 unsigned int float_size=H5Tget_size(base_type_id);
688 if(float_size!=4){
689 HDF5_error_exit("HDF5:read: error, datatype size mismatch!");
690 }
691 H5Tclose(type_id);
692 H5Tclose(base_type_id);
693
694 //perform the actual read:
695 rd(dataname,datum,H5T_COMPOUND,true,false);
696 }
697
698 template<>void HDF5::read< PScalar< PScalar< RComplex<double> > > >(const std::string& dataname, ComplexD& datum){
699 //get type:
700 hid_t type_id;
701 readPrepare(dataname,type_id);
702 if(type_id<0){
703 HDF5_error_exit("HDF5::read: error, cannot open datatype!");
704 }
705 hid_t base_type_id;
706 if(!checkComplexType(type_id,base_type_id)){
707 HDF5_error_exit("HDF5::read: error, datatype mismatch!");
708 }
709 unsigned int float_size=H5Tget_size(base_type_id);
710 if(float_size!=8){
711 HDF5_error_exit("HDF5:read: error, datatype size mismatch!");
712 }
713 H5Tclose(type_id);
714 H5Tclose(base_type_id);
715
716 //perform the actual read:
717 rd(dataname,datum,H5T_COMPOUND,true,false);
718 }
719
720 //array types
721 template<>void HDF5::read< PScalar< PScalar< RComplex<float> > > >(const std::string& dataname, multi1d<ComplexF>& datum){
722 //get type:
723 hid_t type_id;
724 readPrepare(dataname,type_id);
725 if(type_id<0){
726 HDF5_error_exit("HDF5::read: error, cannot open datatype!");
727 }
728 hid_t base_type_id;
729 if(!checkComplexType(type_id,base_type_id)){
730 HDF5_error_exit("HDF5::read: error, datatype mismatch!");
731 }
732 unsigned int float_size=H5Tget_size(base_type_id);
733 if(float_size!=4){
734 HDF5_error_exit("HDF5:read: error, datatype size mismatch!");
735 }
736 H5Tclose(type_id);
737 H5Tclose(base_type_id);
738
739 //perform the actual read:
740 rd(dataname,datum,H5T_COMPOUND,true,false);
741 }
742
743 template<>void HDF5::read< PScalar< PScalar< RComplex<double> > > >(const std::string& dataname, multi1d<ComplexD>& datum){
744 //get type:
745 hid_t type_id;
746 readPrepare(dataname,type_id);
747 if(type_id<0){
748 HDF5_error_exit("HDF5::read: error, cannot open datatype!");
749 }
750 hid_t base_type_id;
751 if(!checkComplexType(type_id,base_type_id)){
752 HDF5_error_exit("HDF5::read: error, datatype mismatch!");
753 }
754 unsigned int float_size=H5Tget_size(base_type_id);
755 if(float_size!=8){
756 HDF5_error_exit("HDF5:read: error, datatype size mismatch!");
757 }
758 H5Tclose(type_id);
759 H5Tclose(base_type_id);
760
761 //perform the actual read:
762 rd(dataname,datum,H5T_COMPOUND,true,false);
763 }
764
765 //***********************************************************************************************************************************
766 //***********************************************************************************************************************************
767 //READING Lattice Types:
768 //***********************************************************************************************************************************
769 //***********************************************************************************************************************************
770 void HDF5::readPrepareLattice(const std::string& name, hid_t& type_id, multi1d<ullong>& sizes){
771 //determine whether there is a lattice with the specified name:
772 bool exists=objectExists(current_group,name);
773 if(!exists){
774 HDF5_error_exit("HDF5::read: error, dataset does not exists!");
775 }
776 H5O_info_t objinfo;
777 herr_t errhandle=H5Oget_info_by_name(current_group,name.c_str(),&objinfo,H5P_DEFAULT);
778 if(objinfo.type!=H5O_TYPE_DATASET){
779 HDF5_error_exit("HDF5::read: error, "+name+" exists but it is not a dataset!");
780 }
781
782 //determine dimension of dataset:
783 //open dataset and filespace:
784 hid_t dset_id=H5Dopen(current_group,name.c_str(),H5P_DEFAULT);
785 if(dset_id<0){
786 HDF5_error_exit("HDF5::read: error, cannot open dataset!");
787 }
788 type_id=H5Dget_type(dset_id);
789
790 //get extent of dataset:
791 hid_t filespace=H5Dget_space(dset_id);
792 int Ndims=H5Sget_simple_extent_ndims(filespace);
793 sizes.resize(Ndims);
794 hsize_t* inttoken=new hsize_t[Ndims];
795 errhandle=H5Sget_simple_extent_dims(filespace,inttoken, NULL);
796 for(unsigned int dd=0; dd<Ndims; dd++){
797 sizes[dd]=inttoken[dd];
798 }
799 delete [] inttoken;
800
801 H5Sclose(filespace);
802 H5Dclose(dset_id);
803
804 //prefetch for faster I/O:
805 if(!isprefetched && sizes.size()>1) prefetchLatticeCoordinates();
806 }
807
808 void HDF5::readLattice(const std::string& name, const hid_t& type_id,
809 const hid_t& base_type_id, const ullong& obj_size,
810 const ullong& tot_size, char* buf, bool invert_order) {
811 // determine local sizes
812 const int mynode = Layout::nodeNumber();
813 const int nodeSites = Layout::sitesOnNode();
814
815 // set up and create hyperslap:
816 unsigned int dimensions;
817 if(obj_size>1) dimensions=Nd+1;
818 else dimensions=Nd;
819 hsize_t* node_offset = new hsize_t[dimensions];
820 hsize_t* offset = new hsize_t[dimensions];
821 hsize_t* total_count = new hsize_t[dimensions];
822 hsize_t* dim_size = new hsize_t[dimensions];
823
824 //set up the chunks every node has to read. If the ordering is row-major or column-major
825 if(invert_order){
826 // reorder such that x is fastest: necessary for chroma HDF5-file formats
827 for(unsigned int i = 0; i < Nd; i ++) {
828 dim_size[i] = total_count[i] = Layout::subgridLattSize()[(Nd - 1) - i];
829 offset[i] = node_offset[i] = Layout::nodeCoord()[(Nd - 1) - i] * total_count[i];
830 } // for
831 }
832 else{
833 //do not reorder, read as-is. necessary for other formats, such as QLUA
834 for(unsigned int i = 0; i < Nd; i ++) {
835 dim_size[i] = total_count[i] = Layout::subgridLattSize()[i];
836 offset[i] = node_offset[i] = Layout::nodeCoord()[i] * total_count[i];
837 } // for
838 }
839 if(obj_size>1){
840 dim_size[Nd] = total_count[Nd] = obj_size;
841 offset[Nd] = node_offset[Nd] = 0;
842 }
843
844 //create property list for readin:
845 hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
846 H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
847
848 // open dataset and filespace and create hyperslab:
849 hid_t dset_id = H5Dopen(current_group, name.c_str(), H5P_DEFAULT);
850 hid_t filespace = H5Dget_space(dset_id);
851
852 unsigned int hdf5_float_size=H5Tget_size(base_type_id);
853 if( (hdf5_float_size!=4) && (hdf5_float_size!=8) ){
854 HDF5_error_exit("HDF5Reader::read: error, invalid base datatype while trying to read lattice!");
855 }
856 if(H5Tget_class(base_type_id)!=H5T_FLOAT){
857 HDF5_error_exit("HDF5Reader::read: error, invalid base datatype while trying to read lattice!");
858 }
859
860 //32 GB MPIO pointer protection
861 hid_t err;
862 size_t two_gb = (size_t) 2 * 1024 * 1024 * 1024;
863 size_t total_size = tot_size;
864 unsigned int blocks = 1;
865 while( (total_size * hdf5_float_size) > two_gb) {
866 dim_size[0] = dim_size[0] >> 1;
867 total_size = total_size >> 1;
868 blocks = blocks << 1;
869 } // while
870
871 //allocate buffers and do some error handling
872 char* buf_small = new (std::nothrow) char[total_size*hdf5_float_size];
873 if(buf_small == 0x0) {
874 HDF5_error_exit("Unable to allocate buf\n");
875 } // if
876
877 hsize_t rank = static_cast<hsize_t>(dimensions);
878 hid_t memspace = H5Screate_simple(rank, dim_size, NULL);
879 hid_t nat_type_id=H5Tget_native_type(type_id,H5T_DIR_ASCEND);
880
881 // read:
882 for(int i = 0; i < blocks; ++ i) {
883 offset[0] = node_offset[0] + i * dim_size[0];
884 H5Sselect_hyperslab(filespace, H5S_SELECT_SET, const_cast<const hsize_t*>(offset),
885 NULL, const_cast<const hsize_t*>(dim_size), NULL);
886 err = H5Dread(dset_id, nat_type_id, memspace, filespace, plist_id, static_cast<void*>(buf_small));
887 for(ullong j = 0; j < (total_size*hdf5_float_size); j++){
888 (buf + i * (total_size*hdf5_float_size) )[j] = buf_small[j];
889 }
890 } // for
891
892 //do local transposition if necessary:
893 if(!invert_order){
894 //get subgrid sizes
895 multi1d<int> locsizes=Layout::subgridLattSize();
896 int dstruct_size=tot_size/nodeSites;
897
898 //create temporary buffer
899 char* tmpbuf=new char[nodeSites*dstruct_size*hdf5_float_size];
900 for(unsigned int x=0; x<locsizes[0]; x++){
901 for(unsigned int y=0; y<locsizes[1]; y++){
902 for(unsigned int z=0; z<locsizes[2]; z++){
903 for(unsigned int t=0; t<locsizes[3]; t++){
904 //transpose from reversed input order to chroma input order
905 memcpy(&tmpbuf[dstruct_size*hdf5_float_size*(x+locsizes[0]*(y+locsizes[1]*(z+locsizes[2]*t)))],
906 &buf[dstruct_size*hdf5_float_size*(t+locsizes[3]*(z+locsizes[2]*(y+locsizes[1]*x)))],
907 dstruct_size*hdf5_float_size);
908 }
909 }
910 }
911 }
912 memcpy(buf,tmpbuf,nodeSites*dstruct_size*hdf5_float_size);
913 delete [] tmpbuf;
914 }
915
916 delete [] buf_small;
917 delete [] dim_size;
918 delete [] total_count;
919 delete [] offset;
920 delete [] node_offset;
921 H5Pclose(plist_id);
922 H5Sclose(memspace);
923 H5Sclose(filespace);
924 H5Tclose(nat_type_id);
925 H5Dclose(dset_id);
926 } // readLattice()
927
928
929 //read LatticeColorMatrix
930 template<>void HDF5::read< PScalar< PColorMatrix< RComplex<REAL64>, 3> > >(const std::string& name,
931 LatticeColorMatrixD3& field,
932 const HDF5Base::accessmode& accmode){
933 StopWatch swatch_prepare, swatch_datatypes, swatch_reorder, swatch_read;
934
935 bool invert_order;
936 switch(accmode){
937 case HDF5Base::transpose_order:
938 invert_order=true;
939 break;
940 case HDF5Base::maintain_order:
941 invert_order=false;
942 break;
943 }
944
945 //read dataset extents:
946 if(profile) swatch_prepare.start();
947 multi1d<ullong> sizes;
948 hid_t type_id;
949 readPrepareLattice(name,type_id,sizes);
950 if(profile) swatch_prepare.stop();
951
952 //do some datatype sanity checks and get the basis precision:
953 if(profile) swatch_datatypes.start();
954 hid_t base_type_id;
955 if( !checkColorMatrixType(type_id,Nc,base_type_id) ){
956 HDF5_error_exit("HDF5::read: object is not a color matrix!");
957 }
958
959 ullong hdf5_float_size=H5Tget_size(base_type_id);
960 ullong field_float_size=sizeof(REAL64);
961
962 if(sizes.size()!=Nd){
963 HDF5_error_exit("HDF5::read: error, wrong dimensionality!");
964 return;
965 }
966 if(invert_order){
967 for(unsigned int dd=0; dd<Nd; dd++){
968 if(sizes[Nd-dd-1]!=Layout::lattSize()[dd]){
969 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
970 }
971 }
972 }
973 else{
974 for(unsigned int dd=0; dd<Nd; dd++){
975 if(sizes[dd]!=Layout::lattSize()[dd]){
976 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
977 }
978 }
979 }
980 if(profile) swatch_datatypes.stop();
981
982 //determine local sizes, allocate memory and read
983 if(profile) swatch_read.start();
984 const int mynode=Layout::nodeNumber();
985 const int nodeSites = Layout::sitesOnNode();
986 size_t obj_size=sizeof(ColorMatrixD3)/field_float_size;
987 size_t tot_size = nodeSites*obj_size;
988 char* buf = new(std::nothrow) char[tot_size*hdf5_float_size];
989 if( buf == 0x0 ) {
990 HDF5_error_exit("Unable to allocate buf\n");
991 }
992
993 readLattice(name,type_id,base_type_id,1,tot_size,buf,invert_order);
994 H5Tclose(type_id);
995 H5Tclose(base_type_id);
996 if(profile) swatch_read.stop();
997
998 //put lattice into u-field and reconstruct as well as reorder them on the fly:
999 // Reconstruct the gauge field
1000 if(profile) swatch_reorder.start();
1001 /*#pragma omp parallel for firstprivate(nodeSites,obj_size,float_size) shared(buf,field)
1002 for(unsigned int run=0; run<nodeSites; run++){
1003 memcpy(&(field.elem(reordermap[run])),reinterpret_cast<char*>(buf+run*obj_size),float_size*obj_size);
1004 }*/
1005 if(hdf5_float_size==field_float_size){
1006 CvtToLayout(field,reinterpret_cast<void*>(buf),nodeSites,sizeof(ColorMatrixD3));
1007 }
1008 else{
1009 REAL64* tmpbuf=new REAL64[tot_size];
1010 REAL32 tmpfloat;
1011 for(unsigned int i=0; i<tot_size; i++){
1012 memcpy(&tmpfloat,&buf[i*hdf5_float_size],sizeof(REAL32));
1013 tmpbuf[i]=static_cast< REAL64 >(tmpfloat);
1014 }
1015 CvtToLayout(field,reinterpret_cast<void*>(tmpbuf),nodeSites,sizeof(ColorMatrixD3));
1016 delete [] tmpbuf;
1017 }
1018 delete [] buf;
1019 if(profile) swatch_reorder.stop();
1020
1021 if(profile){
1022 QDPIO::cout << "HDF5-I/O statistics. Read:" << std::endl;
1023 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
1024 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
1025 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
1026 QDPIO::cout << "\t read: " << swatch_read.getTimeInSeconds() << " s." << std::endl;
1027 QDPIO::cout << "\t MB read: " << Layout::vol()*obj_size*static_cast<int>(hdf5_float_size)/1024/1024 << std::endl;
1028 }
1029 }
1030
1031 //read LatticeDiracPropagatorF3
1032 template<>void HDF5::read< PSpinMatrix< PColorMatrix< RComplex<REAL32>, 3>, 4> >(const std::string& name,
1034 const HDF5Base::accessmode& accmode){
1035 StopWatch swatch_prepare, swatch_datatypes, swatch_reorder, swatch_read;
1036
1037 bool invert_order;
1038 switch(accmode){
1039 case HDF5Base::transpose_order:
1040 invert_order=true;
1041 break;
1042 case HDF5Base::maintain_order:
1043 invert_order=false;
1044 break;
1045 }
1046
1047 //read dataset extents:
1048 if(profile) swatch_prepare.start();
1049 multi1d<ullong> sizes;
1050 hid_t type_id;
1051 readPrepareLattice(name,type_id,sizes);
1052 if(profile) swatch_prepare.stop();
1053
1054 //do some datatype sanity checks and get the basis precision:
1055 if(profile) swatch_datatypes.start();
1056 hid_t base_type_id;
1057 if( !checkDiracPropagatorType(type_id,Ns,Nc,base_type_id) ){
1058 HDF5_error_exit("HDF5::read: object is not a dirac propagator!");
1059 }
1060
1061 ullong hdf5_float_size=H5Tget_size(base_type_id);
1062 ullong field_float_size=sizeof(REAL32);
1063 if(sizes.size()!=Nd){
1064 HDF5_error_exit("HDF5::read: error, wrong dimensionality!");
1065 return;
1066 }
1067 if(invert_order){
1068 for(unsigned int dd=0; dd<Nd; dd++){
1069 if(sizes[Nd-dd-1]!=Layout::lattSize()[dd]){
1070 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
1071 }
1072 }
1073 }
1074 else{
1075 for(unsigned int dd=0; dd<Nd; dd++){
1076 if(sizes[dd]!=Layout::lattSize()[dd]){
1077 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
1078 }
1079 }
1080 }
1081 if(profile) swatch_datatypes.stop();
1082
1083 //determine local sizes, allocate memory and read
1084 if(profile) swatch_read.start();
1085 const int mynode=Layout::nodeNumber();
1086 const int nodeSites = Layout::sitesOnNode();
1087 size_t obj_size=sizeof(DiracPropagatorF3)/field_float_size;
1088 size_t tot_size = nodeSites*obj_size;
1089 char* buf = new(std::nothrow) char[tot_size*hdf5_float_size];
1090 if( buf == 0x0 ) {
1091 HDF5_error_exit("Unable to allocate buf\n");
1092 }
1093
1094 readLattice(name,type_id,base_type_id,1,tot_size,buf,invert_order);
1095 H5Tclose(type_id);
1096 H5Tclose(base_type_id);
1097 if(profile) swatch_read.stop();
1098
1099 //put lattice into u-field and reconstruct as well as reorder them on the fly:
1100 // Reconstruct the gauge field
1101 if(profile) swatch_reorder.start();
1102 if(hdf5_float_size==field_float_size){
1103 //convert layout directly
1104 CvtToLayout(field,reinterpret_cast<void*>(buf),nodeSites,sizeof(DiracPropagatorF3));
1105 }
1106 else{
1107 //convert precision first
1108 REAL32* tmpbuf=new REAL32[tot_size];
1109 REAL64 tmpfloat;
1110 for(unsigned int i=0; i<tot_size; i++){
1111 memcpy(&tmpfloat,&buf[i*hdf5_float_size],sizeof(REAL64));
1112 tmpbuf[i]=static_cast< REAL32 >(tmpfloat);
1113 }
1114 CvtToLayout(field,reinterpret_cast<void*>(tmpbuf),nodeSites,sizeof(DiracPropagatorF3));
1115 delete [] tmpbuf;
1116 }
1117 delete [] buf;
1118 if(profile) swatch_reorder.stop();
1119
1120 if(profile){
1121 QDPIO::cout << "HDF5-I/O statistics. Read:" << std::endl;
1122 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
1123 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
1124 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
1125 QDPIO::cout << "\t read: " << swatch_read.getTimeInSeconds() << " s." << std::endl;
1126 QDPIO::cout << "\t MB read: " << Layout::vol()*sizeof(DiracPropagatorF3)/1024/1024 << std::endl;
1127 }
1128 }
1129
1130 //read LatticeDiracPropagatorD3
1131 template<>void HDF5::read< PSpinMatrix< PColorMatrix< RComplex<REAL64>, 3>, 4> >(const std::string& name,
1133 const HDF5Base::accessmode& accmode){
1134 StopWatch swatch_prepare, swatch_datatypes, swatch_reorder, swatch_read;
1135
1136 bool invert_order;
1137 switch(accmode){
1138 case HDF5Base::transpose_order:
1139 invert_order=true;
1140 break;
1141 case HDF5Base::maintain_order:
1142 invert_order=false;
1143 break;
1144 }
1145
1146 //read dataset extents:
1147 if(profile) swatch_prepare.start();
1148 multi1d<ullong> sizes;
1149 hid_t type_id;
1150 readPrepareLattice(name,type_id,sizes);
1151 if(profile) swatch_prepare.stop();
1152
1153 //do some datatype sanity checks and get the basis precision:
1154 if(profile) swatch_datatypes.start();
1155 hid_t base_type_id;
1156 if( !checkDiracPropagatorType(type_id,Ns,Nc,base_type_id) ){
1157 HDF5_error_exit("HDF5::read: object is not a dirac propagator!");
1158 }
1159
1160 ullong hdf5_float_size=H5Tget_size(base_type_id);
1161 ullong field_float_size=sizeof(REAL64);
1162 if(sizes.size()!=Nd){
1163 HDF5_error_exit("HDF5::read: error, wrong dimensionality!");
1164 return;
1165 }
1166 if(invert_order){
1167 for(unsigned int dd=0; dd<Nd; dd++){
1168 if(sizes[Nd-dd-1]!=Layout::lattSize()[dd]){
1169 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
1170 }
1171 }
1172 }
1173 else{
1174 for(unsigned int dd=0; dd<Nd; dd++){
1175 if(sizes[dd]!=Layout::lattSize()[dd]){
1176 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
1177 }
1178 }
1179 }
1180 if(profile) swatch_datatypes.stop();
1181
1182 //determine local sizes, allocate memory and read
1183 if(profile) swatch_read.start();
1184 const int mynode=Layout::nodeNumber();
1185 const int nodeSites = Layout::sitesOnNode();
1186 size_t obj_size=sizeof(DiracPropagatorD3)/field_float_size;
1187 size_t tot_size = nodeSites*obj_size;
1188 char* buf = new(std::nothrow) char[tot_size*hdf5_float_size];
1189 if( buf == 0x0 ) {
1190 HDF5_error_exit("Unable to allocate buf\n");
1191 }
1192
1193 readLattice(name,type_id,base_type_id,1,tot_size,buf,invert_order);
1194 H5Tclose(type_id);
1195 H5Tclose(base_type_id);
1196 if(profile) swatch_read.stop();
1197
1198 //put lattice into u-field and reconstruct as well as reorder them on the fly:
1199 // Reconstruct the gauge field
1200 if(profile) swatch_reorder.start();
1201 if(hdf5_float_size==field_float_size){
1202 //convert layout directly
1203 CvtToLayout(field,reinterpret_cast<void*>(buf),nodeSites,sizeof(DiracPropagatorD3));
1204 }
1205 else{
1206 //convert precision first
1207 REAL64* tmpbuf=new REAL64[tot_size];
1208 REAL32 tmpfloat;
1209 for(unsigned int i=0; i<tot_size; i++){
1210 memcpy(&tmpfloat,&buf[i*hdf5_float_size],sizeof(REAL32));
1211 tmpbuf[i]=static_cast< REAL64 >(tmpfloat);
1212 }
1213 CvtToLayout(field,reinterpret_cast<void*>(tmpbuf),nodeSites,sizeof(DiracPropagatorD3));
1214 delete [] tmpbuf;
1215 }
1216 delete [] buf;
1217 if(profile) swatch_reorder.stop();
1218
1219 if(profile){
1220 QDPIO::cout << "HDF5-I/O statistics. Read:" << std::endl;
1221 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
1222 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
1223 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
1224 QDPIO::cout << "\t read: " << swatch_read.getTimeInSeconds() << " s." << std::endl;
1225 QDPIO::cout << "\t MB read: " << Layout::vol()*sizeof(DiracPropagatorD3)/1024/1024 << std::endl;
1226 }
1227 }
1228
1229 //QDP Lattice IO:
1230 template<>void HDF5::read< PScalar< PColorMatrix< RComplex<REAL64>, 3> > >(const std::string& name,
1231 multi1d<LatticeColorMatrixD3>& fieldarray,
1232 const HDF5Base::accessmode& accmode){
1233 StopWatch swatch_prepare, swatch_datatypes, swatch_reorder, swatch_read;
1234
1235 bool invert_order;
1236 switch(accmode){
1237 case HDF5Base::transpose_order:
1238 invert_order=true;
1239 break;
1240 case HDF5Base::maintain_order:
1241 invert_order=false;
1242 break;
1243 }
1244
1245 //read dataset extents:
1246 if(profile) swatch_prepare.start();
1247
1248 multi1d<ullong> sizes;
1249 hid_t type_id;
1250 readPrepareLattice(name,type_id,sizes);
1251 unsigned int arr_size=sizes[Nd];
1252 if(profile) swatch_prepare.stop();
1253
1254 //do some datatype sanity checks and get the basis precision:
1255 if(profile) swatch_datatypes.start();
1256 hid_t base_type_id;
1257 if( !checkColorMatrixType(type_id,Nc,base_type_id) ){
1258 HDF5_error_exit("HDF5::read: object is not a color matrix!");
1259 }
1260
1261 ullong hdf5_float_size=H5Tget_size(base_type_id);
1262 ullong field_float_size=sizeof(REAL64);
1263 if(sizes.size()!=(Nd+1)){
1264 HDF5_error_exit("HDF5::read: error, wrong dimensionality!");
1265 return;
1266 }
1267 if(invert_order){
1268 for(unsigned int dd=0; dd<Nd; dd++){
1269 if(sizes[Nd-dd-1]!=Layout::lattSize()[dd]){
1270 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
1271 }
1272 }
1273 }
1274 else{
1275 for(unsigned int dd=0; dd<Nd; dd++){
1276 if(sizes[dd]!=Layout::lattSize()[dd]){
1277 HDF5_error_exit("HDF5::read: mismatching lattice extents.");
1278 }
1279 }
1280 }
1281 fieldarray.resize(sizes[Nd]);
1282 if(profile) swatch_datatypes.stop();
1283
1284 //determine local sizes, allocate memory and read
1285 if(profile) swatch_read.start();
1286 const int mynode=Layout::nodeNumber();
1287 const int nodeSites = Layout::sitesOnNode();
1288 size_t obj_size=sizeof(ColorMatrixD3)/field_float_size;
1289 size_t tot_size = nodeSites*arr_size*obj_size;
1290 char* buf = new(std::nothrow) char[tot_size*hdf5_float_size];
1291 if( buf == 0x0 ) {
1292 HDF5_error_exit("Unable to allocate buf\n");
1293 }
1294
1295 readLattice(name,type_id,base_type_id,sizes[Nd],tot_size,buf,invert_order);
1296 H5Tclose(type_id);
1297 H5Tclose(base_type_id);
1298 if(profile) swatch_read.stop();
1299
1300 //put lattice into u-field and reconstruct as well as reorder them on the fly:
1301 // Reconstruct the gauge field
1302 if(profile) swatch_reorder.start();
1303 /*#pragma omp parallel for firstprivate(nodeSites,arr_size,obj_size,float_size) shared(buf,field)
1304 for(unsigned int run=0; run<nodeSites; run++){
1305 for(unsigned int dd=0; dd<arr_size; dd++){
1306 memcpy(&(field[dd].elem(reordermap[run])),reinterpret_cast<char*>(buf+(dd+arr_size*run)*obj_size),float_size*obj_size);
1307 }
1308 }*/
1309 if(hdf5_float_size==field_float_size){
1310 //convert layout directly
1311 CvtToLayout(fieldarray,reinterpret_cast<void*>(buf),nodeSites,arr_size,sizeof(ColorMatrixD3));
1312 }
1313 else{
1314 //convert precision first
1315 REAL64* tmpbuf=new REAL64[tot_size];
1316 REAL32 tmpfloat;
1317 for(unsigned int i=0; i<tot_size; i++){
1318 memcpy(&tmpfloat,&buf[i*hdf5_float_size],sizeof(REAL32));
1319 tmpbuf[i]=static_cast< REAL64 >(tmpfloat);
1320 }
1321 CvtToLayout(fieldarray,reinterpret_cast<void*>(tmpbuf),nodeSites,arr_size,sizeof(ColorMatrixD3));
1322 delete [] tmpbuf;
1323 }
1324 delete [] buf;
1325 if(profile) swatch_reorder.stop();
1326
1327 if(profile){
1328 QDPIO::cout << "HDF5-I/O statistics. Read:" << std::endl;
1329 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
1330 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
1331 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
1332 QDPIO::cout << "\t read: " << swatch_read.getTimeInSeconds() << " s." << std::endl;
1333 QDPIO::cout << "\t MB read: " << Layout::vol()*fieldarray.size()*obj_size*static_cast<int>(hdf5_float_size)/1024/1024 << std::endl;
1334 }
1335 }
1336
1337 //***********************************************************************************************************************************
1338 //***********************************************************************************************************************************
1339 //READING Qlua Lattice Types:
1340 //***********************************************************************************************************************************
1341 //***********************************************************************************************************************************
1342 //LatticeColorMatrixD3
1343 void HDF5::readQlua(const std::string& name, multi1d<LatticeColorMatrixD3>& field){
1344 StopWatch swatch_complete;
1345
1346 QDPIO::cout << "Reading Qlua config ..." << std::flush;
1347
1348 //do some checks
1349 if(profile) swatch_complete.start();
1350 if(!objectExists(current_group,name)){
1351 HDF5_error_exit("HDF5::readQlua: error, configuration "+name+"does not exist!");
1352 }
1353 //check if the object is a group and if it contains the datasets 0 to Nd-1 but not more:
1354 H5O_info_t objinfo;
1355 herr_t errhandle=H5Oget_info_by_name(current_group,name.c_str(),&objinfo,H5P_DEFAULT);
1356 if(objinfo.type!=H5O_TYPE_GROUP){
1357 HDF5_error_exit("HDF5::readQlua: error, "+name+" exists but it is not a Qlua config!");
1358 }
1359 //check if datasets from 0 to Nd-1 exist:
1360 for(unsigned int dd=0; dd<Nd; dd++){
1361 std::stringstream stream;
1362 stream << name << "/" << dd;
1363 std::string dname=stream.str();
1364 if(!objectExists(current_group,dname)){
1365 HDF5_error_exit("HDF5::readQlua: error, "+name+" exists but it is not a Qlua config! Dataset "+dname+" was not found!");
1366 }
1367 }
1368 //now, check if there is another entry and whether all these entries are LatticeColorMatrices:
1369 std::stringstream stream;
1370 stream << name << "/" << Nd;
1371 std::string dname=stream.str();
1372 if(objectExists(current_group,dname)){
1373 HDF5_error_exit("HDF5::readQlua: error, "+name+" exists but it is not a Qlua config!");
1374 }
1375
1376 //read LatticeColorMatrices: set invert_ordering to false, as Qlua stores the fields with x fastest
1377 field.resize(Nd);
1378 for(unsigned int dd=0; dd<Nd; dd++){
1379 std::stringstream stream;
1380 stream << name << "/" << dd;
1381 std::string dname=stream.str();
1382 read(dname,field[dd],HDF5Base::maintain_order);
1383 }
1384 if(profile) swatch_complete.stop();
1385 QDPIO::cout << "done!" << std::endl;
1386
1387 if(profile){
1388 QDPIO::cout << "HDF5-I/O statistics for Qlua-read: " << std::endl;
1389 QDPIO::cout << "\t total: " << swatch_complete.getTimeInSeconds() << " s." << std::endl;
1390 QDPIO::cout << "\t MB read: " << Layout::vol()*field.size()*sizeof(ColorMatrixD3)/1024/1024 << std::endl;
1391 }
1392 }
1393
1394 //LatticePropagatorD3
1395 void HDF5::readQlua(const std::string& name, LatticeDiracPropagatorD3& prop){
1396 StopWatch swatch_complete;
1397
1398 QDPIO::cout << "Reading Qlua propagator ..." << std::flush;
1399
1400 //do some checks
1401 if(profile) swatch_complete.start();
1402 if(!objectExists(current_group,name)){
1403 HDF5_error_exit("HDF5::readQlua: error, dataset "+name+"does not exist!");
1404 }
1405 //check if the object is a group:
1406 H5O_info_t objinfo;
1407 herr_t errhandle=H5Oget_info_by_name(current_group,name.c_str(),&objinfo,H5P_DEFAULT);
1408 if(!objectExists(current_group,name)){
1409 HDF5_error_exit("HDF5::readQlua: error, "+name+" exists but it is not a Qlua propagator!");
1410 }
1411
1412 /*
1413 //check if the Datatype named .DiracPropagatorDouble3 exists and open it:
1414 if(!objectExists(file_id,".DiracPropagatorDouble3")){
1415 HDF5_error_exit("HDF5::readQlua: error, the required committed type .DiracPropagatorDouble3 does not exist!");
1416 }
1417 errhandle=H5Oget_info_by_name(file_id,".DiracPropagatorDouble3",&objinfo,H5P_DEFAULT);
1418 if(objinfo.type!=H5O_TYPE_NAMED_DATATYPE){
1419 HDF5_error_exit("HDF5::readQlua: error, object .DiracPropagatorDouble3 exists but it is not a committed datatype!");
1420 }
1421 hid_t type_id=H5Topen(file_id,".DiracPropagatorDouble3",H5P_DEFAULT);
1422 hid_t base_type_id;
1423 if(!checkDiracPropagatorType(type_id, 4, 3, base_type_id)){
1424 HDF5_error_exit("HDF5::readQlua: error, the type .DiracPropagatorDouble3 is not formatted as expected!");
1425 }
1426 H5Tclose(type_id);
1427 //finally, check the base type id:
1428 unsigned int float_size=H5Tget_size(base_type_id);
1429 if(float_size!=8){
1430 HDF5_error_exit("HDF5:readQlua: error, datatype size mismatch!");
1431 }
1432 H5Tclose(base_type_id);
1433 */
1434
1435 //read LatticePropagator: do not invert the order
1436 read(name,prop,HDF5Base::maintain_order);
1437 if(profile) swatch_complete.stop();
1438 QDPIO::cout << "done!" << std::endl;
1439
1440 if(profile){
1441 QDPIO::cout << "HDF5-I/O statistics for Qlua-read: " << std::endl;
1442 QDPIO::cout << "\t total: " << swatch_complete.getTimeInSeconds() << " s." << std::endl;
1443 QDPIO::cout << "\t MB read: " << Layout::vol()*sizeof(DiracPropagatorD3)/1024/1024 << std::endl;
1444 }
1445 }
1446
1447 //--------------------------------------------------------------------------------
1448 //--------------------------------------------------------------------------------
1449 // Reader class
1450 //--------------------------------------------------------------------------------
1451 //--------------------------------------------------------------------------------
1452 //constructors:
1453 HDF5Reader::HDF5Reader() : HDF5() {};
1454
1455 HDF5Reader::HDF5Reader(const long int& blocksize, const long int& maxalign) : HDF5(blocksize, maxalign) {};
1456
1457 HDF5Reader::HDF5Reader(const std::string& filename): HDF5(){
1458 open(filename);
1459 };
1460
1461 //destructors:
1462 HDF5Reader::~HDF5Reader(){
1463 close();
1464 };
1465
1466 void HDF5Reader::open(const std::string& filename){
1467 //close actual file:
1468 close();
1469
1470 //check if file exists first:
1471 bool exists=check_exists(filename);
1472 if(!exists){
1473 HDF5_error_exit("HDF5Reader::open: error, file does not exist!");
1474 }
1475
1476 //create file access property control list:
1477 hid_t fapl_id = H5Pcreate(H5P_FILE_ACCESS);
1478
1479 //switch on LUSTRE optimizations:
1480 if(stripesize>0){
1481 //memory alignment:
1482 H5Pset_alignment(fapl_id,maxalign,stripesize);
1483 }
1484
1485 //create MPI_IO accessor:
1486 MPI_Info info = MPI_INFO_NULL;
1487 QMP_get_hidden_comm(QMP_comm_get_default(),reinterpret_cast<void**>(&mpicomm));
1488 H5Pset_fapl_mpio(fapl_id,*mpicomm, info);
1489
1490
1491 file_id=H5Fopen(filename.c_str(),H5F_ACC_RDONLY,fapl_id);
1492 H5Pclose(fapl_id);
1493
1494 if(file_id<0){
1495 HDF5_error_exit("HDF5Reader::open: could not open file "+filename+" for reading. Some error occured!");
1496 }
1497
1498 current_group=file_id;
1499 }
1500
1501
1502 //--------------------------------------------------------------------------------
1503 //--------------------------------------------------------------------------------
1504 // Writer class
1505 //--------------------------------------------------------------------------------
1506 //--------------------------------------------------------------------------------
1507 //constructors
1508 HDF5Writer::HDF5Writer() : HDF5() {};
1509
1510 HDF5Writer::HDF5Writer(const long int& stripesizee, const long int& maxalignn) : HDF5(stripesizee,maxalignn) {};
1511
1512 HDF5Writer::HDF5Writer(const std::string& filename, const HDF5Base::writemode& mode) : HDF5(){
1513 open(filename,mode);
1514 };
1515
1516 //destructors
1517 HDF5Writer::~HDF5Writer(){
1518 close();
1519 };
1520
1521 //member functions
1522 void HDF5Writer::open(const std::string& filename, const HDF5Base::writemode& mode){
1523 //close actual file if open:
1524 close();
1525
1526 //try to create file. First create file access control list:
1527 hid_t fapl_id = H5Pcreate(H5P_FILE_ACCESS);
1528 hid_t fcpl_id = H5Pcreate(H5P_FILE_CREATE);
1529
1530 //switch on LUSTRE optimizations:
1531 if(stripesize>0){
1532 //memory alignment:
1533 H5Pset_alignment(fapl_id,maxalign,stripesize);
1534 }
1535
1536 //activate parallel IO:
1537 //this is MPI stuff and breaks the whole USQCD software design paradigm, however, at the
1538 //moment it is not possible to get MPI objects from QMP routines. This will be implemented
1539 //as soon as this hack works:
1540 MPI_Info info = MPI_INFO_NULL;
1541 QMP_get_hidden_comm(QMP_comm_get_default(),reinterpret_cast<void**>(&mpicomm));
1542 H5Pset_fapl_mpio(fapl_id,*mpicomm,info);
1543
1544 //collective file creation/opening:
1545 //on node 0, test if file exists:
1546 bool exists=false;
1547 if(Layout::nodeNumber()==0){
1548 std::ifstream input(filename.c_str(),std::ios_base::binary);
1549 if(input.good()){
1550 exists=true;
1551 }
1552 input.close();
1553 }
1554 QDPInternal::broadcast(exists);
1555
1556 if(exists){
1557 //file exists, check if it is an HDF5 file:
1558 htri_t ex=H5Fis_hdf5(filename.c_str());
1559 if(ex>0){
1560 //file is existing HDF5 file:
1561 if(mode&HDF5Base::trunc){
1562 //it exists and should be overwritten, so truncate it:
1563 file_id=H5Fcreate(filename.c_str(),H5F_ACC_TRUNC,fcpl_id,fapl_id);
1564 }
1565 else{
1566 //it exists and should not be overwritten, so open it rw:
1567 file_id=H5Fopen(filename.c_str(),H5F_ACC_RDWR,fapl_id);
1568 }
1569 }
1570 else{
1571 //file exists but is not HDF5 file:
1572 if(mode&HDF5Base::trunc){
1573 //it should be overwritten, so truncate it:
1574 remove(filename.c_str());
1575 file_id=H5Fcreate(filename.c_str(),H5F_ACC_EXCL,fcpl_id,fapl_id);
1576 }
1577 else{
1578 HDF5_error_exit("HDF5Writer::open: error, file "+filename+" already exists and is not HDF5 file! Please use overwrite=true to create anew file!");
1579 }
1580 }
1581 }
1582 else{
1583 //it does not exists, so simply create a new file:
1584 file_id=H5Fcreate(filename.c_str(),H5F_ACC_EXCL,fcpl_id,fapl_id);
1585 }
1586 H5Pclose(fapl_id);
1587 H5Pclose(fcpl_id);
1588
1589 if(file_id<0){
1590 HDF5_error_exit("HDF5Writer::open: could not open file "+filename+" for writing.");
1591 }
1592
1593 current_group=file_id;
1594 }
1595
1596 //create a new group inside current one w/o steping into it:
1597 void HDF5Writer::mkdir(const ::std::string& name){
1598 if(name!="" && name!="/"){
1599 std::string cwd=pwd();
1600 push(name);
1601 cd(cwd);
1602 }
1603 }
1604
1605 //create a new group inside current one and step into it:
1606 void HDF5Writer::push(const std::string& name){
1607 std::vector<std::string> dirlist=splitPathname(name);
1608
1609 //if this is an absolute path, go to top dir first:
1610 if(name.find_first_of("/")==0){
1611 cd("/");
1612 }
1613
1614 //plist identifier
1615 hid_t plist_id = H5Pcreate (H5P_GROUP_CREATE);
1616
1617 //create groups iteratively:
1618 hid_t last_group;
1619 for(unsigned int i=0; i<static_cast<unsigned int>(dirlist.size()); i++){
1620 last_group=current_group;
1621
1622 //check if group exists:
1623 htri_t ex=H5Lexists(last_group,dirlist[i].c_str(),H5P_DEFAULT);
1624 if(ex==0){
1625 //does not exists, create a link:
1626 current_group=H5Gcreate(last_group,dirlist[i].c_str(),H5P_DEFAULT,plist_id,H5P_DEFAULT);
1627 }
1628 else{
1629 //the link exists, check if the group exists as well:
1630 ex=H5Oexists_by_name(last_group,dirlist[i].c_str(),H5P_DEFAULT);
1631 if(ex==0){
1632 QDPIO::cout << "HDF5Writer::push: path traversal error! You seem to have dangling links, I will try to clean up" << std::endl;
1633 //delete old link:
1634 H5Ldelete(last_group,dirlist[i].c_str(),H5P_DEFAULT);
1635 //create new group:
1636 current_group=H5Gcreate(last_group,dirlist[i].c_str(),H5P_DEFAULT,plist_id,H5P_DEFAULT);
1637 }
1638 else{
1639 //check if object is a group: if yes, step into it, if no, cast an error:
1640 H5O_info_t objinfo;
1641 herr_t errhandle=H5Oget_info_by_name(last_group,dirlist[i].c_str(),&objinfo,H5P_DEFAULT);
1642 if(errhandle<0){
1643 QDPIO::cout << "HDF5Writer::push: error, cannot get properties of " << name << "!" << std::endl;
1644 current_group=last_group;
1645 break;
1646 }
1647 if(objinfo.type!=H5O_TYPE_GROUP){
1648 QDPIO::cout << "HDF5Writer::push: error, the object " << dirlist[i] << " in " << getNameById(last_group) << " already exists and is not a group!" << std::endl;
1649 current_group=last_group;
1650 break;
1651 }
1652 else{
1653 current_group=H5Gopen(last_group,dirlist[i].c_str(),H5P_DEFAULT);
1654 }
1655 }
1656 }
1657
1658 if(current_group<0){
1659 HDF5_error_exit("HDF5Writer::push: something went wrong, aborting!");
1660 }
1661 }
1662
1663 //close plist-identifier
1664 H5Pclose(plist_id);
1665 }
1666
1667 //attribute handling:
1668 herr_t HDF5Writer::rmAtt(hid_t location_id, const char *attr_name, const H5A_info_t* attrinfo, void* opdata){
1669 return H5Adelete(location_id,attr_name);
1670 }
1671
1672 void HDF5Writer::deleteAllAttributes(const std::string& obj_name){
1673 H5Aiterate_by_name(current_group,obj_name.c_str(),H5_INDEX_NAME,H5_ITER_NATIVE,NULL,rmAtt,NULL,H5P_DEFAULT);
1674 }
1675
1676 void HDF5Writer::deleteAttribute(const std::string& obj_name, const std::string& attr_name){
1677 std::string oname(obj_name), aname(attr_name);
1678 htri_t exists=H5Aexists_by_name(current_group,oname.c_str(),aname.c_str(),H5P_DEFAULT);
1679 if(exists!=1){
1680 QDPIO::cout << "HDF5Writer::deleteAttribute: error, attribute does not exists!" << std::endl; \
1681 return;
1682 }
1683 herr_t errhandle=H5Adelete_by_name(current_group,oname.c_str(),aname.c_str(),H5P_DEFAULT);
1684 }
1685
1686 //***********************************************************************************************************************************
1687 //***********************************************************************************************************************************
1688 //WRITING ATTRIBUTES
1689 //***********************************************************************************************************************************
1690 //***********************************************************************************************************************************
1691
1692 template<typename ctype>
1693 bool get_global(ctype& global, const ctype& local)
1694 {
1695 global = local;
1696 QDPInternal::broadcast(global);
1697
1698 int chcksum = global != local ? 1 : 0;
1699 QDPInternal::globalSum(chcksum);
1700
1701 return chcksum==0;
1702 }
1703
1704 template<typename ctype>
1705 bool get_global(multi1d<ctype>& global, const multi1d<ctype>& local)
1706 {
1707 int s0;
1708 int chcksum = not get_global(s0, local.size());
1709
1710 if (Layout::nodeNumber()==0)
1711 global = local;
1712 else
1713 global.resize(s0);
1714
1715 QDPInternal::broadcast(&(global[0]), sizeof(ctype)*s0);
1716
1717 if (chcksum!=0) // this is and MUST be global, otherwise a possible deadlock
1718 return false;
1719
1720 for (int j=0; j<s0; ++j)
1721 chcksum += global[j] != local[j] ? 1 : 0;
1722
1723 QDPInternal::globalSum(chcksum);
1724
1725 return chcksum==0;
1726 }
1727
1728 template bool get_global<short>( multi1d<short>& global, const multi1d<short>& local);
1729 template bool get_global<int>( multi1d<int>& global, const multi1d<int>& local);
1730 template bool get_global<unsigned short>( multi1d<unsigned short>& global, const multi1d<unsigned short>& local);
1731 template bool get_global<unsigned int>( multi1d<unsigned int>& global, const multi1d<unsigned int>& local);
1732 template bool get_global<unsigned long>( multi1d<unsigned long>& global, const multi1d<unsigned long>& local);
1733 template bool get_global<unsigned long long>( multi1d<unsigned long long>& global, const multi1d<unsigned long long>& local);
1734 template bool get_global<float>( multi1d<float>& global, const multi1d<float>& local);
1735 template bool get_global<double>( multi1d<double>& global, const multi1d<double>& local);
1736
1737 //single datum
1738 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const short& datum, const HDF5Base::writemode& mode){
1739 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_SHORT,mode);
1740 }
1741
1742 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const unsigned short& datum, const HDF5Base::writemode& mode){
1743 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_USHORT,mode);
1744 }
1745
1746 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const int& datum, const HDF5Base::writemode& mode){
1747 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_INT,mode);
1748 }
1749
1750 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const unsigned int& datum, const HDF5Base::writemode& mode){
1751 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_UINT,mode);
1752 }
1753
1754 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const unsigned long long& datum, const HDF5Base::writemode& mode){
1755 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_ULLONG,mode);
1756 }
1757
1758 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const float& datum, const HDF5Base::writemode& mode){
1759 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_FLOAT,mode);
1760 }
1761
1762 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const double& datum, const HDF5Base::writemode& mode){
1763 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_DOUBLE,mode);
1764 }
1765
1766 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const std::string& datum, const HDF5Base::writemode& mode){
1767 std::string oname(obj_name), aname(attr_name);
1768
1769 std::string datum_0;
1770 if (not get_global(datum_0,datum)) {
1771 QDPIO::cerr << "HDF5Writer::writeAttribute() warning: " << obj_name
1772 << ".attrib(" << attr_name << ") was NOT global. Using node=0 value now." << std::endl;
1773 }
1774
1775 bool exists=objectExists(current_group,oname);
1776 if(!exists){
1777 HDF5_error_exit("HDF5Writer::writeAttribute: error, object "+oname+" you try to write attribute to does not exists!");
1778 }
1779
1780 hid_t ex=H5Aexists_by_name(current_group,oname.c_str(),aname.c_str(),H5P_DEFAULT);
1781 if(ex==1){
1782 if(!(mode&HDF5Base::trunc)){
1783 HDF5_error_exit("HDF5Writer::writeAttribute: error, attribute "+aname+" already exists!");
1784 }
1785 herr_t errhandle=H5Adelete_by_name(current_group,oname.c_str(),aname.c_str(),H5P_DEFAULT);
1786 }
1787
1788 if (datum_0.length()+1>64*1024) {
1789 QDPIO::cerr << "HDF5Writer::writeAttribute() error: " << obj_name
1790 << ".attrib(" << attr_name
1791 << ") exceeds the maximum hdf5 attrib size (64kB)." << std::endl;
1792
1793 HDF5_error_exit("bad string attrib write");
1794 }
1795
1796 //create string datatytpe and set encoding to UTF-8:
1797 hid_t typid=H5Tcreate(H5T_STRING,datum_0.length()+1);
1798 H5Tset_cset(typid,H5T_CSET_UTF8);
1799 //create space:
1800 hid_t attr_space_id=H5Screate(H5S_SCALAR);
1801 hid_t attr_id=H5Acreate_by_name(current_group,oname.c_str(),aname.c_str(),typid,attr_space_id,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
1802 H5Sclose(attr_space_id);
1803
1804 //write:
1805 H5Awrite(attr_id,typid,reinterpret_cast<void*>(const_cast<char*>(datum_0.c_str())));
1806 H5Aclose(attr_id);
1807 H5Tclose(typid);
1808 }
1809
1810 //array.
1811 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const multi1d<short>& datum, const HDF5Base::writemode& mode){
1812 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_SHORT,mode);
1813 }
1814
1815 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const multi1d<unsigned short>& datum, const HDF5Base::writemode& mode){
1816 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_USHORT,mode);
1817 }
1818
1819 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const multi1d<int>& datum, const HDF5Base::writemode& mode){
1820 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_INT,mode);
1821 }
1822
1823 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const multi1d<unsigned int>& datum, const HDF5Base::writemode& mode){
1824 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_UINT,mode);
1825 }
1826
1827 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const multi1d<unsigned long long>& datum, const HDF5Base::writemode& mode){
1828 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_ULLONG,mode);
1829 }
1830
1831 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const multi1d<float>& datum, const HDF5Base::writemode& mode){
1832 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_FLOAT,mode);
1833 }
1834
1835 void HDF5Writer::writeAttribute(const std::string& obj_name, const std::string& attr_name, const multi1d<double>& datum, const HDF5Base::writemode& mode){
1836 wtAtt(obj_name,attr_name,datum,H5T_NATIVE_DOUBLE,mode);
1837 }
1838
1839 //***********************************************************************************************************************************
1840 //***********************************************************************************************************************************
1841 //WRITING DATASETS
1842 //***********************************************************************************************************************************
1843 //***********************************************************************************************************************************
1844 //single datum
1845 void HDF5Writer::write(const std::string& obj_name, const short& datum, const HDF5Base::writemode& mode){
1846 wt(obj_name,datum,H5T_NATIVE_SHORT,mode);
1847 }
1848
1849 void HDF5Writer::write(const std::string& obj_name, const unsigned short& datum, const HDF5Base::writemode& mode){
1850 wt(obj_name,datum,H5T_NATIVE_USHORT,mode);
1851 }
1852
1853 void HDF5Writer::write(const std::string& obj_name, const int& datum, const HDF5Base::writemode& mode){
1854 wt(obj_name,datum,H5T_NATIVE_INT,mode);
1855 }
1856
1857 void HDF5Writer::write(const std::string& obj_name, const unsigned int& datum, const HDF5Base::writemode& mode){
1858 wt(obj_name,datum,H5T_NATIVE_UINT,mode);
1859 }
1860
1861 void HDF5Writer::write(const std::string& obj_name, const unsigned long long& datum, const HDF5Base::writemode& mode){
1862 wt(obj_name,datum,H5T_NATIVE_ULLONG,mode);
1863 }
1864
1865 void HDF5Writer::write(const std::string& obj_name, const float& datum, const HDF5Base::writemode& mode){
1866 wt(obj_name,datum,H5T_NATIVE_FLOAT,mode);
1867 }
1868
1869 void HDF5Writer::write(const std::string& obj_name, const double& datum, const HDF5Base::writemode& mode){
1870 wt(obj_name,datum,H5T_NATIVE_DOUBLE,mode);
1871 }
1872
1873 void HDF5Writer::write(const std::string& dataname, const std::string& datum, const HDF5Base::writemode& mode){
1874 std::string dname(dataname);
1875 std::string datum_0;
1876 if (not get_global(datum_0,datum)) {
1877 QDPIO::cerr << "HDF5Writer::write() warning: " << dataname << " was NOT global. Using node=0 value now." << std::endl;
1878 }
1879
1880 bool exists=objectExists(current_group,dname);
1881 if(exists){
1882 if(!(mode&HDF5Base::trunc)){
1883 QDPIO::cout << "HDF5Writer::write: error, object named " << dname << " already exists!" << std::endl;
1884 return;
1885 }
1886 H5O_info_t objinfo;
1887 hid_t errhandle=H5Oget_info_by_name(current_group,dname.c_str(),&objinfo,H5P_DEFAULT);
1888 if(objinfo.type!=H5O_TYPE_DATASET){
1889 QDPIO::cout << "HDF5Writer::write: error, object you try to write does already exist and is of different type!" << std::endl;
1890 return;
1891 }
1892 errhandle=H5Ldelete(current_group,dname.c_str(),H5P_DEFAULT);
1893 }
1894
1895 //create string datatytpe and set encoding to UTF-8:
1896 hid_t dataid, spaceid, typid=H5Tcreate(H5T_STRING,datum_0.length()+1);
1897 H5Tset_cset(typid,H5T_CSET_UTF8);
1898 //create space:
1899 spaceid=H5Screate(H5S_SCALAR);
1900 dataid=H5Dcreate(current_group,dname.c_str(),typid,spaceid,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
1901 H5Sclose(spaceid);
1902
1903 hid_t plist_id = H5Pcreate (H5P_DATASET_XFER);
1904 H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
1905 H5Dwrite(dataid,typid,H5S_ALL,H5S_ALL,plist_id,reinterpret_cast<void*>(const_cast<char*>(datum_0.c_str())));
1906 H5Pclose(plist_id);
1907 H5Tclose(typid);
1908 H5Dclose(dataid);
1909 }
1910
1911 //array:
1912 //1D
1913 void HDF5Writer::write(const std::string& obj_name, const multi1d<short>& datum, const HDF5Base::writemode& mode){
1914 wt(obj_name,datum,H5T_NATIVE_SHORT,mode);
1915 }
1916
1917 void HDF5Writer::write(const std::string& obj_name, const multi1d<unsigned short>& datum, const HDF5Base::writemode& mode){
1918 wt(obj_name,datum,H5T_NATIVE_USHORT,mode);
1919 }
1920
1921 void HDF5Writer::write(const std::string& obj_name, const multi1d<int>& datum, const HDF5Base::writemode& mode){
1922 wt(obj_name,datum,H5T_NATIVE_INT,mode);
1923 }
1924
1925 void HDF5Writer::write(const std::string& obj_name, const multi1d<unsigned int>& datum, const HDF5Base::writemode& mode){
1926 wt(obj_name,datum,H5T_NATIVE_UINT,mode);
1927 }
1928
1929 void HDF5Writer::write(const std::string& obj_name, const multi1d<unsigned long long>& datum, const HDF5Base::writemode& mode){
1930 wt(obj_name,datum,H5T_NATIVE_ULLONG,mode);
1931 }
1932
1933 void HDF5Writer::write(const std::string& obj_name, const multi1d<float>& datum, const HDF5Base::writemode& mode){
1934 wt(obj_name,datum,H5T_NATIVE_FLOAT,mode);
1935 }
1936
1937 void HDF5Writer::write(const std::string& obj_name, const multi1d<double>& datum, const HDF5Base::writemode& mode){
1938 wt(obj_name,datum,H5T_NATIVE_DOUBLE,mode);
1939 }
1940 //2D
1941 void HDF5Writer::write(const std::string& obj_name, const multi2d<short>& datum, const HDF5Base::writemode& mode){
1942 wt(obj_name,datum,H5T_NATIVE_SHORT,mode);
1943 }
1944
1945 void HDF5Writer::write(const std::string& obj_name, const multi2d<unsigned short>& datum, const HDF5Base::writemode& mode){
1946 wt(obj_name,datum,H5T_NATIVE_USHORT,mode);
1947 }
1948
1949 void HDF5Writer::write(const std::string& obj_name, const multi2d<int>& datum, const HDF5Base::writemode& mode){
1950 wt(obj_name,datum,H5T_NATIVE_INT,mode);
1951 }
1952
1953 void HDF5Writer::write(const std::string& obj_name, const multi2d<unsigned int>& datum, const HDF5Base::writemode& mode){
1954 wt(obj_name,datum,H5T_NATIVE_UINT,mode);
1955 }
1956
1957 void HDF5Writer::write(const std::string& obj_name, const multi2d<unsigned long long>& datum, const HDF5Base::writemode& mode){
1958 wt(obj_name,datum,H5T_NATIVE_ULLONG,mode);
1959 }
1960
1961 void HDF5Writer::write(const std::string& obj_name, const multi2d<float>& datum, const HDF5Base::writemode& mode){
1962 wt(obj_name,datum,H5T_NATIVE_FLOAT,mode);
1963 }
1964
1965 void HDF5Writer::write(const std::string& obj_name, const multi2d<double>& datum, const HDF5Base::writemode& mode){
1966 wt(obj_name,datum,H5T_NATIVE_DOUBLE,mode);
1967 }
1968
1969 //***********************************************************************************************************************************
1970 //***********************************************************************************************************************************
1971 //WRITING Compound types:
1972 //***********************************************************************************************************************************
1973 //***********************************************************************************************************************************
1974 //complex types:
1975 //Single element:
1976 template<>
1977 void HDF5Writer::write< PScalar< PScalar< RComplex<float> > > >(const std::string& dataname, const ComplexF& datum, const HDF5Base::writemode& mode){
1978 hid_t type_id;
1979 bool exists=objectExists(file_id,".ComplexFloat");
1980 if(!exists){
1981 type_id=createComplexType(sizeof(REAL32));
1982 commitType(".ComplexFloat",type_id);
1983 }
1984 else{
1985 type_id=H5Topen(file_id,".ComplexFloat",H5P_DEFAULT);
1986 if(type_id<0){
1987 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
1988 }
1989 }
1990
1991 //perform the actual write:
1992 wt(dataname,datum,type_id,mode);
1993 H5Tclose(type_id);
1994 }
1995
1996 template<>
1997 void HDF5Writer::write< PScalar< PScalar< RComplex<double> > > >(const std::string& dataname, const ComplexD& datum, const HDF5Base::writemode& mode){
1998 hid_t type_id;
1999 bool exists=objectExists(file_id,".ComplexDouble");
2000 if(!exists){
2001 type_id=createComplexType(sizeof(REAL64));
2002 commitType(".ComplexDouble",type_id);
2003 }
2004 else{
2005 type_id=H5Topen(file_id,".ComplexDouble",H5P_DEFAULT);
2006 if(type_id<0){
2007 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2008 }
2009 }
2010
2011 //perform the actual write:
2012 wt(dataname,datum,type_id,mode);
2013 H5Tclose(type_id);
2014 }
2015
2016 //array types
2017 template<>
2018 void HDF5Writer::write< PScalar< PScalar< RComplex<float> > > >(const std::string& dataname, const multi1d<ComplexF>& datum, const HDF5Base::writemode& mode){
2019 hid_t type_id;
2020 bool exists=objectExists(file_id,".ComplexFloat");
2021 if(!exists){
2022 type_id=createComplexType(sizeof(REAL32));
2023 commitType(".ComplexFloat",type_id);
2024 }
2025 else{
2026 type_id=H5Topen(file_id,".ComplexFloat",H5P_DEFAULT);
2027 if(type_id<0){
2028 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2029 }
2030 }
2031
2032 //perform the actual write:
2033 wt(dataname,datum,type_id,mode);
2034 H5Tclose(type_id);
2035 }
2036
2037 template<>
2038 void HDF5Writer::write< PScalar< PScalar< RComplex<double> > > >(const std::string& dataname, const multi1d<ComplexD>& datum, const HDF5Base::writemode& mode){
2039 hid_t type_id;
2040 bool exists=objectExists(file_id,".ComplexDouble");
2041 if(!exists){
2042 type_id=createComplexType(sizeof(REAL64));
2043 commitType(".ComplexDouble",type_id);
2044 }
2045 else{
2046 type_id=H5Topen(file_id,".ComplexDouble",H5P_DEFAULT);
2047 if(type_id<0){
2048 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2049 }
2050 }
2051
2052 //perform the actual write:
2053 wt(dataname,datum,type_id,mode);
2054 H5Tclose(type_id);
2055 }
2056
2057 //ColorMatrix:
2058 template<>
2059 void HDF5Writer::write< PScalar< PColorMatrix< RComplex<REAL32>, 3> > >(const std::string& dataname, const ColorMatrixF3& datum, const HDF5Base::writemode& mode){
2060 //first get complex type:
2061 hid_t complex_id, colmat_id;
2062 bool exists=objectExists(file_id,".ComplexFloat");
2063 if(!exists){
2064 complex_id=createComplexType(sizeof(REAL32));
2065 commitType(".ComplexFloat",complex_id);
2066 }
2067 else{
2068 complex_id=H5Topen(file_id,".ComplexFloat",H5P_DEFAULT);
2069 if(complex_id<0){
2070 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2071 }
2072 }
2073
2074 //get color matrix type:
2075 exists=objectExists(file_id,".ColorMatrixFloat3");
2076 if(!exists){
2077 colmat_id=createColorMatrixType(complex_id,Nc);
2078 commitType(".ColorMatrixFloat3",colmat_id);
2079 }
2080 else{
2081 colmat_id=H5Topen(file_id,".ColorMatrixFloat3",H5P_DEFAULT);
2082 if(colmat_id<0){
2083 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2084 }
2085 }
2086
2087 wt(dataname,datum,colmat_id,mode);
2088 H5Tclose(complex_id);
2089 H5Tclose(colmat_id);
2090 }
2091
2092 template<>
2093 void HDF5Writer::write< PScalar< PColorMatrix< RComplex<REAL64>, 3> > >(const std::string& dataname, const ColorMatrixD3& datum, const HDF5Base::writemode& mode){
2094 //first get complex type:
2095 hid_t complex_id, colmat_id;
2096 bool exists=objectExists(file_id,".ComplexDouble");
2097 if(!exists){
2098 complex_id=createComplexType(sizeof(REAL64));
2099 commitType(".ComplexDouble",complex_id);
2100 }
2101 else{
2102 complex_id=H5Topen(file_id,".ComplexDouble",H5P_DEFAULT);
2103 if(complex_id<0){
2104 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2105 }
2106 }
2107
2108 //get color matrix type:
2109 exists=objectExists(file_id,".ColorMatrixDouble3");
2110 if(!exists){
2111 colmat_id=createColorMatrixType(complex_id,Nc);
2112 commitType(".ColorMatrixDouble3",colmat_id);
2113 }
2114 else{
2115 colmat_id=H5Topen(file_id,".ColorMatrixDouble3",H5P_DEFAULT);
2116 if(colmat_id<0){
2117 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2118 }
2119 }
2120
2121 wt(dataname,datum,colmat_id,mode);
2122 H5Tclose(complex_id);
2123 H5Tclose(colmat_id);
2124 }
2125
2126 //***********************************************************************************************************************************
2127 //***********************************************************************************************************************************
2128 //Lattice IO
2129 //***********************************************************************************************************************************
2130 //***********************************************************************************************************************************
2131 //helper routines for Lattice field I/O:
2132 void HDF5Writer::writePrepare(const std::string& name, const HDF5Base::writemode& mode){
2133 //before writing is performed, check if dataset exists:
2134 herr_t errhandle;
2135
2136 bool exists=objectExists(current_group,name);
2137 if(exists){
2138 if(!(mode&HDF5Base::trunc)){
2139 HDF5_error_exit("HDF5Writer::write: error, dataset already exists and you specified not to overwrite!");
2140 }
2141 H5O_info_t objinfo;
2142 errhandle=H5Oget_info_by_name(current_group,name.c_str(),&objinfo,H5P_DEFAULT);
2143 if(objinfo.type!=H5O_TYPE_DATASET){
2144 HDF5_error_exit("HDF5Writer::write: error, the object you try to write does already exist and is of different type!");
2145 }
2146 //delete attributes (if present) and unlink storage:
2147 deleteAllAttributes(name);
2148 errhandle=H5Ldelete(current_group,name.c_str(),H5P_DEFAULT);
2149 }
2150
2151 //prefetch for faster I/O:
2152 if(!isprefetched) prefetchLatticeCoordinates();
2153 }
2154
2155 void HDF5Writer::writeLattice(const std::string& name, const hid_t& datatype, const ullong& obj_size, char* buf){
2156 //writing out reordered data array:
2157 //determine the dimension of the array: this is useful since not for all classes a write routine will be implemented. In that case,
2158 //it falls back to a floating point array. In other cases, more sophisticated datatypes will be written and stored
2159 unsigned int dimension;
2160 if(obj_size>1) dimension=Nd+1;
2161 else dimension=Nd;
2162
2163 //create dataspace and dataset:
2164 hsize_t rank = static_cast<hsize_t>(dimension);
2165 hsize_t* spacesize = new hsize_t[dimension];
2166 //reorder such that x is fastest:
2167 for(unsigned int i = 0; i < Nd; i ++) {
2168 spacesize[i] = Layout::lattSize()[(Nd - 1) - i];
2169 } // for
2170 if(obj_size>1) spacesize[Nd] = obj_size;
2171 hid_t filespace = H5Screate_simple(rank, const_cast<const hsize_t*>(spacesize), NULL);
2172 hid_t dcpl_id = H5Pcreate(H5P_DATASET_CREATE);
2173
2174 //hyperslab creation and selection:
2175 hsize_t* node_offset = new hsize_t[dimension];
2176 hsize_t* offset = new hsize_t[dimension];
2177 hsize_t* total_count = new hsize_t[dimension];
2178 hsize_t* dim_size = new hsize_t[dimension];
2179 //reorder such that x is fastest:
2180 for(unsigned int i = 0; i < Nd; ++ i){
2181 dim_size[i] = total_count[i] = Layout::subgridLattSize()[(Nd - 1) - i];
2182 offset[i] = node_offset[i] = Layout::nodeCoord()[(Nd - 1) - i] * total_count[i];
2183 } // for
2184 if(obj_size>1){
2185 dim_size[Nd] = total_count[Nd] = obj_size;
2186 offset[Nd] = node_offset[Nd] = 0;
2187 }
2188
2189 //LUSTRE optimization:
2190 if(stripesize > 0) H5Pset_chunk(dcpl_id,dimension, total_count);
2191
2192 //create dataset:
2193 hid_t dset_id = H5Dcreate(current_group, name.c_str(), datatype, filespace,
2194 H5P_DEFAULT, dcpl_id, H5P_DEFAULT);
2195 H5Sclose(filespace);
2196 H5Pclose(dcpl_id);
2197 delete [] spacesize;
2198
2199 //create property list for writeout:
2200 hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
2201 H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
2202
2203 size_t typesize=H5Tget_size(datatype);
2204 size_t total_size = typesize * obj_size * Layout::sitesOnNode();
2205 size_t two_gb = (size_t) 2 * 1024 * 1024 * 1024;
2206
2207 int power = 0;
2208 while(total_size * sizeof(REAL) > two_gb) {
2209 dim_size[0] = dim_size[0] >> 1; // assuming there is enough room in this dim
2210 total_size = total_size >> 1;
2211 ++ power;
2212 } // while
2213 int blocks = 1 << power; // number of blocks to write
2214
2215 // create memspace
2216 hid_t memspace = H5Screate_simple(rank, dim_size, NULL);
2217 filespace = H5Dget_space(dset_id);
2218
2219 for(int i = 0; i < blocks; ++ i) {
2220 offset[0] = node_offset[0] + i * dim_size[0];
2221 H5Sselect_hyperslab(filespace, H5S_SELECT_SET, const_cast<const hsize_t*>(offset), NULL,
2222 const_cast<const hsize_t*>(dim_size), NULL);
2223 // write
2224 H5Dwrite(dset_id, datatype, memspace, filespace, plist_id, buf + i * total_size);
2225 } // for
2226
2227 // cleaning up
2228 delete [] dim_size;
2229 delete [] total_count;
2230 delete [] offset;
2231 delete [] node_offset;
2232
2233 H5Dclose(dset_id);
2234 H5Sclose(filespace);
2235 H5Sclose(memspace);
2236 H5Pclose(plist_id);
2237 } // writeLattice()
2238
2239 //float lattice color matrix:
2240 template<>
2241 void HDF5Writer::write< PScalar< PColorMatrix< RComplex<REAL32>, 3> > >(const std::string& name, const LatticeColorMatrixF3& field, const HDF5Base::writemode& mode){
2242 StopWatch swatch_prepare, swatch_reorder, swatch_write, swatch_datatypes;
2243
2244 //before writing is performed, check if dataset exists:
2245 if(profile) swatch_prepare.start();
2246 writePrepare(name,mode);
2247 if(profile) swatch_prepare.stop();
2248
2249 //color matrix datatype:
2250 if(profile) swatch_datatypes.start();
2251 hid_t complex_id, colmat_id;
2252 bool exists=objectExists(file_id,".ComplexFloat");
2253 if(!exists){
2254 complex_id=createComplexType(sizeof(REAL32));
2255 commitType(".ComplexFloat",complex_id);
2256 }
2257 else{
2258 complex_id=H5Topen(file_id,".ComplexFloat",H5P_DEFAULT);
2259 if(complex_id<0){
2260 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2261 }
2262 }
2263 exists=objectExists(file_id,".ColorMatrixFloat3");
2264 if(!exists){
2265 colmat_id=createColorMatrixType(complex_id,Nc);
2266 commitType(".ColorMatrixFloat3",colmat_id);
2267 }
2268 else{
2269 colmat_id=H5Topen(file_id,".ColorMatrixFloat3",H5P_DEFAULT);
2270 if(colmat_id<0){
2271 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2272 }
2273 }
2274 if(profile) swatch_datatypes.stop();
2275
2276 //get node information:
2277 if(profile) swatch_reorder.start();
2278 const int mynode=Layout::nodeNumber();
2279 const int nodeSites = Layout::sitesOnNode();
2280
2281 //copy buffer into data
2282 size_t float_size=sizeof(REAL32);
2283 size_t obj_size=sizeof(ColorMatrixF3)/float_size;
2284 REAL32* buf=new REAL32[nodeSites*obj_size];
2285 /*#pragma omp parallel for firstprivate(nodeSites,obj_size,float_size) shared(buf,field)
2286 for(unsigned int run=0; run<nodeSites; run++){
2287 memcpy(reinterpret_cast<char*>(buf+run*obj_size),&(field.elem(reordermap[run])),float_size*obj_size);
2288 }*/
2289 CvtToHost(reinterpret_cast<void*>(buf),field,nodeSites,float_size*obj_size);
2290 if(profile) swatch_reorder.stop();
2291
2292 //write out the stuff:
2293 if(profile) swatch_write.start();
2294 writeLattice(name,colmat_id,1,reinterpret_cast<char*>(buf));
2295
2296 //clean up
2297 H5Tclose(colmat_id);
2298 H5Tclose(complex_id);
2299 delete [] buf;
2300 if(profile) swatch_write.stop();
2301
2302 if(profile){
2303 QDPIO::cout << "HDF5-I/O statistics. Write:" << std::endl;
2304 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
2305 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
2306 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
2307 QDPIO::cout << "\t write: " << swatch_write.getTimeInSeconds() << " s." << std::endl;
2308 QDPIO::cout << "\t MB written: " << Layout::vol()*sizeof(ColorMatrixF3)/1024/1024 << std::endl;
2309 }
2310 }
2311
2312 //double lattice color matrix:
2313 template<>
2314 void HDF5Writer::write< PScalar< PColorMatrix< RComplex<REAL64>, 3> > >(const std::string& name, const LatticeColorMatrixD3& field, const HDF5Base::writemode& mode)
2315 {
2316 StopWatch swatch_prepare, swatch_reorder, swatch_write, swatch_datatypes;
2317
2318 //before writing is performed, check if dataset exists:
2319 if(profile) swatch_prepare.start();
2320 writePrepare(name,mode);
2321 if(profile) swatch_prepare.stop();
2322
2323 //color matrix datatype:
2324 if(profile) swatch_datatypes.start();
2325 hid_t complex_id, colmat_id;
2326 bool exists=objectExists(file_id,".ComplexDouble");
2327 if(!exists){
2328 complex_id=createComplexType(sizeof(REAL64));
2329 commitType(".ComplexDouble",complex_id);
2330 }
2331 else{
2332 complex_id=H5Topen(file_id,".ComplexDouble",H5P_DEFAULT);
2333 if(complex_id<0){
2334 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2335 }
2336 }
2337 exists=objectExists(file_id,".ColorMatrixDouble3");
2338 if(!exists){
2339 colmat_id=createColorMatrixType(complex_id,Nc);
2340 commitType(".ColorMatrixDouble3",colmat_id);
2341 }
2342 else{
2343 colmat_id=H5Topen(file_id,".ColorMatrixDouble3",H5P_DEFAULT);
2344 if(colmat_id<0){
2345 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2346 }
2347 }
2348 if(profile) swatch_datatypes.stop();
2349
2350 //get node information:
2351 if(profile) swatch_reorder.start();
2352 const int mynode=Layout::nodeNumber();
2353 const int nodeSites = Layout::sitesOnNode();
2354
2355 //copy buffer into data
2356 size_t float_size=sizeof(REAL64);
2357 size_t obj_size=sizeof(ColorMatrixD3)/float_size;
2358 REAL64* buf=new REAL64[nodeSites*obj_size];
2359 /*#pragma omp parallel for firstprivate(nodeSites,obj_size,float_size) shared(buf,field)
2360 for(unsigned int run=0; run<nodeSites; run++){
2361 memcpy(reinterpret_cast<char*>(buf+run*obj_size),&(field.elem(reordermap[run])),float_size*obj_size);
2362 }*/
2363 CvtToHost(reinterpret_cast<void*>(buf),field,nodeSites,float_size*obj_size);
2364 if(profile) swatch_reorder.stop();
2365
2366 //write out the stuff:
2367 if(profile) swatch_write.start();
2368 writeLattice(name,colmat_id,1,reinterpret_cast<char*>(buf));
2369
2370 //clean up
2371 H5Tclose(colmat_id);
2372 H5Tclose(complex_id);
2373 delete [] buf;
2374 if(profile) swatch_write.stop();
2375
2376 if(profile){
2377 QDPIO::cout << "HDF5-I/O statistics. Write:" << std::endl;
2378 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
2379 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
2380 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
2381 QDPIO::cout << "\t write: " << swatch_write.getTimeInSeconds() << " s." << std::endl;
2382 QDPIO::cout << "\t MB written: " << Layout::vol()*sizeof(ColorMatrixD3)/1024/1024 << std::endl;
2383 }
2384 }
2385
2386 //float lattice propagator:
2387 template<>
2388 void HDF5Writer::write< PSpinMatrix< PColorMatrix< RComplex<REAL32>, 3>, 4> >(const std::string& name, const LatticePropagatorF3& field, const HDF5Base::writemode& mode){
2389 StopWatch swatch_prepare, swatch_reorder, swatch_write, swatch_datatypes;
2390 QDPIO::cout << "\t USING LatticePropagatorF3 writer" << std::endl;
2391 //before writing is performed, check if dataset exists:
2392 if(profile) swatch_prepare.start();
2393 writePrepare(name,mode);
2394 if(profile) swatch_prepare.stop();
2395
2396 //color matrix datatype:
2397 if(profile) swatch_datatypes.start();
2398 hid_t complex_id, colmat_id, prop_id;
2399 bool exists=objectExists(file_id,".ComplexFloat");
2400 if(!exists){
2401 complex_id=createComplexType(sizeof(REAL32));
2402 commitType(".ComplexFloat",complex_id);
2403 }
2404 else{
2405 complex_id=H5Topen(file_id,".ComplexFloat",H5P_DEFAULT);
2406 if(complex_id<0){
2407 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2408 }
2409 }
2410 exists=objectExists(file_id,".ColorMatrixFloat3");
2411 if(!exists){
2412 colmat_id=createColorMatrixType(complex_id,Nc);
2413 commitType(".ColorMatrixFloat3",colmat_id);
2414 }
2415 else{
2416 colmat_id=H5Topen(file_id,".ColorMatrixFloat3",H5P_DEFAULT);
2417 if(colmat_id<0){
2418 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2419 }
2420 }
2421 exists=objectExists(file_id,".PropagatorFloat3");
2422 if(!exists){
2423 prop_id=createPropagatorType(colmat_id,Ns);
2424 commitType(".PropagatorFloat3",prop_id);
2425 }
2426 else{
2427 prop_id=H5Topen(file_id,".PropagatorFloat3",H5P_DEFAULT);
2428 if(prop_id<0){
2429 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2430 }
2431 }
2432 if(profile) swatch_datatypes.stop();
2433
2434 //get node information:
2435 if(profile) swatch_reorder.start();
2436 const int mynode=Layout::nodeNumber();
2437 const int nodeSites = Layout::sitesOnNode();
2438
2439 //copy buffer into data
2440 size_t float_size=sizeof(REAL32);
2441 size_t obj_size=sizeof(PropagatorF3)/float_size;
2442 REAL32* buf=new REAL32[nodeSites*obj_size];
2443 /*#pragma omp parallel for firstprivate(nodeSites,obj_size,float_size) shared(buf,field)
2444 for(unsigned int run=0; run<nodeSites; run++){
2445 memcpy(reinterpret_cast<char*>(buf+run*obj_size),&(field.elem(reordermap[run])),float_size*obj_size);
2446 }*/
2447 CvtToHost(reinterpret_cast<void*>(buf),field,nodeSites,float_size*obj_size);
2448 if(profile) swatch_reorder.stop();
2449
2450 //write out the stuff:
2451 if(profile) swatch_write.start();
2452 writeLattice(name,prop_id,1,reinterpret_cast<char*>(buf));
2453
2454 //clean up
2455 H5Tclose(colmat_id);
2456 H5Tclose(complex_id);
2457 delete [] buf;
2458 if(profile) swatch_write.stop();
2459
2460 if(profile){
2461 QDPIO::cout << "HDF5-I/O statistics. Write:" << std::endl;
2462 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
2463 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
2464 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
2465 QDPIO::cout << "\t write: " << swatch_write.getTimeInSeconds() << " s." << std::endl;
2466 QDPIO::cout << "\t MB written: " << Layout::vol()*sizeof(PropagatorF3)/1024/1024 << std::endl;
2467 }
2468 }
2469
2470 //double lattice propagator:
2471 template<>
2472 void HDF5Writer::write< PSpinMatrix< PColorMatrix< RComplex<REAL64>, 3>, 4> >(const std::string& name, const LatticePropagatorD3& field, const HDF5Base::writemode& mode)
2473 {
2474 StopWatch swatch_prepare, swatch_reorder, swatch_write, swatch_datatypes;
2475
2476 //before writing is performed, check if dataset exists:
2477 if(profile) swatch_prepare.start();
2478 writePrepare(name,mode);
2479 if(profile) swatch_prepare.stop();
2480
2481 //color matrix datatype:
2482 if(profile) swatch_datatypes.start();
2483 hid_t complex_id, colmat_id, prop_id;
2484 bool exists=objectExists(file_id,".ComplexDouble");
2485 if(!exists){
2486 complex_id=createComplexType(sizeof(REAL64));
2487 commitType(".ComplexDouble",complex_id);
2488 }
2489 else{
2490 complex_id=H5Topen(file_id,".ComplexDouble",H5P_DEFAULT);
2491 if(complex_id<0){
2492 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2493 }
2494 }
2495 exists=objectExists(file_id,".ColorMatrixDouble3");
2496 if(!exists){
2497 colmat_id=createColorMatrixType(complex_id,Nc);
2498 commitType(".ColorMatrixDouble3",colmat_id);
2499 }
2500 else{
2501 colmat_id=H5Topen(file_id,".ColorMatrixDouble3",H5P_DEFAULT);
2502 if(colmat_id<0){
2503 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2504 }
2505 }
2506 exists=objectExists(file_id,".PropagatorDouble3");
2507 if(!exists){
2508 prop_id=createPropagatorType(colmat_id,Ns);
2509 commitType(".PropagatorDouble3",prop_id);
2510 }
2511 else{
2512 prop_id=H5Topen(file_id,".PropagatorDouble3",H5P_DEFAULT);
2513 if(prop_id<0){
2514 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2515 }
2516 }
2517 if(profile) swatch_datatypes.stop();
2518
2519 //get node information:
2520 if(profile) swatch_reorder.start();
2521 const int mynode=Layout::nodeNumber();
2522 const int nodeSites = Layout::sitesOnNode();
2523
2524 //copy buffer into data
2525 size_t float_size=sizeof(REAL64);
2526 size_t obj_size=sizeof(PropagatorD3)/float_size;
2527 REAL64* buf=new REAL64[nodeSites*obj_size];
2528 /*#pragma omp parallel for firstprivate(nodeSites,obj_size,float_size) shared(buf,field)
2529 for(unsigned int run=0; run<nodeSites; run++){
2530 memcpy(reinterpret_cast<char*>(buf+run*obj_size),&(field.elem(reordermap[run])),float_size*obj_size);
2531 }*/
2532 CvtToHost(reinterpret_cast<void*>(buf),field,nodeSites,float_size*obj_size);
2533 if(profile) swatch_reorder.stop();
2534
2535 //write out the stuff:
2536 if(profile) swatch_write.start();
2537 writeLattice(name,prop_id,1,reinterpret_cast<char*>(buf));
2538
2539 //clean up
2540 H5Tclose(colmat_id);
2541 H5Tclose(complex_id);
2542 delete [] buf;
2543 if(profile) swatch_write.stop();
2544
2545 if(profile){
2546 QDPIO::cout << "HDF5-I/O statistics. Write:" << std::endl;
2547 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
2548 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
2549 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
2550 QDPIO::cout << "\t write: " << swatch_write.getTimeInSeconds() << " s." << std::endl;
2551 QDPIO::cout << "\t MB written: " << Layout::vol()*sizeof(PropagatorD3)/1024/1024 << std::endl;
2552 }
2553 }
2554
2555 //write chroma configuration:
2556 template<>
2557 void HDF5Writer::write< PScalar< PColorMatrix< RComplex<REAL64>, 3> > >(const std::string& name, const multi1d<LatticeColorMatrixD3>& field, const HDF5Base::writemode& mode)
2558 {
2559 StopWatch swatch_prepare, swatch_reorder, swatch_write, swatch_datatypes;
2560
2561 //before writing is performed, check if dataset exists:
2562 if(profile) swatch_prepare.start();
2563 writePrepare(name,mode);
2564 if(profile) swatch_prepare.stop();
2565
2566 //color matrix datatype:
2567 if(profile) swatch_datatypes.start();
2568 hid_t complex_id, colmat_id;
2569 bool exists=objectExists(file_id,".ComplexDouble");
2570 if(!exists){
2571 complex_id=createComplexType(sizeof(REAL64));
2572 commitType(".ComplexDouble",complex_id);
2573 }
2574 else{
2575 complex_id=H5Topen(file_id,".ComplexDouble",H5P_DEFAULT);
2576 if(complex_id<0){
2577 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2578 }
2579 }
2580 exists=objectExists(file_id,".ColorMatrixDouble3");
2581 if(!exists){
2582 colmat_id=createColorMatrixType(complex_id,Nc);
2583 commitType(".ColorMatrixDouble3",colmat_id);
2584 }
2585 else{
2586 colmat_id=H5Topen(file_id,".ColorMatrixDouble3",H5P_DEFAULT);
2587 if(colmat_id<0){
2588 HDF5_error_exit("HDF5Writer::write: error, cannot open committed Datatype!");
2589 }
2590 }
2591 if(profile) swatch_datatypes.stop();
2592
2593 //get node information:
2594 if(profile) swatch_reorder.start();
2595 const int mynode=Layout::nodeNumber();
2596 const int nodeSites = Layout::sitesOnNode();
2597
2598 //copy buffer into data
2599 size_t float_size=sizeof(REAL64);
2600 size_t obj_size=sizeof(ColorMatrixD3)/float_size;
2601 size_t tot_size = nodeSites*field.size()*obj_size;
2602 REAL64* buf=new REAL64[tot_size];
2603 unsigned int fsize=field.size();
2604 /*#pragma omp parallel for firstprivate(nodeSites,fsize,obj_size,float_size) shared(buf,field)
2605 for(unsigned int run=0; run<nodeSites; run++){
2606 for(unsigned int dd=0; dd<fsize; dd++){
2607 memcpy(reinterpret_cast<char*>(buf+(dd+fsize*run)*obj_size),&(field[dd].elem(reordermap[run])),float_size*obj_size);
2608 }
2609 }*/
2610 CvtToHost(reinterpret_cast<void*>(buf),field,nodeSites,fsize,float_size*obj_size);
2611 if(profile) swatch_reorder.stop();
2612
2613 //write out the stuff:
2614 if(profile) swatch_write.start();
2615 writeLattice(name,colmat_id,field.size(),reinterpret_cast<char*>(buf));
2616
2617 //clean up
2618 H5Tclose(colmat_id);
2619 H5Tclose(complex_id);
2620 delete [] buf;
2621 if(profile) swatch_write.stop();
2622
2623 if(profile){
2624 QDPIO::cout << "HDF5-I/O statistics. Write:" << std::endl;
2625 QDPIO::cout << "\t preparing: " << swatch_prepare.getTimeInSeconds() << " s." << std::endl;
2626 QDPIO::cout << "\t datatype-handling: " << swatch_datatypes.getTimeInSeconds() << " s." << std::endl;
2627 QDPIO::cout << "\t reordering: " << swatch_reorder.getTimeInSeconds() << " s." << std::endl;
2628 QDPIO::cout << "\t write: " << swatch_write.getTimeInSeconds() << " s." << std::endl;
2629 QDPIO::cout << "\t MB written: " << Layout::vol()*field.size()*sizeof(ColorMatrixD3)/1024/1024 << std::endl;
2630 }
2631 }
2632
2633 //write FUEL configuration
2634 void HDF5Writer::writeQlua(const std::string& name, const multi1d<LatticeColorMatrixD3>& field, const HDF5Base::writemode& mode){
2635 StopWatch swatch_complete;
2636
2637 QDPIO::cout<< "Writing Qlua config..." << std::flush;
2638 if(profile) swatch_complete.start();
2639 if(field.size()!=Nd) HDF5_error_exit("HDF5Writer::writeQlua: passed vector is not a gauge field!");
2640 if(objectExists(current_group,name)){
2641 if(!(mode&HDF5Base::trunc)){
2642 HDF5_error_exit("HDF5Writer::writeQlua: error, object "+name+" does already exist!");
2643 }
2644 else{
2645 //check if the object is a group and if it contains the datasets 0 to Nd-1:
2646 H5O_info_t objinfo;
2647 herr_t errhandle=H5Oget_info_by_name(current_group,name.c_str(),&objinfo,H5P_DEFAULT);
2648 if(objinfo.type!=H5O_TYPE_GROUP){
2649 HDF5_error_exit("HDF5Writer::writeQlua: error, "+name+" exists but it is not a Qlua config!");
2650 }
2651 //check if datasets from 0 to Nd-1 exist:
2652 for(unsigned int i=0; i<Nd; i++){
2653 std::stringstream stream;
2654 stream << name << "/" << i;
2655 std::string dname=stream.str();
2656 if(!objectExists(current_group,dname)){
2657 HDF5_error_exit("HDF5Writer::writeQlua: error, "+name+" exists but it is not a Qlua config! Dataset "+dname+" was not found!");
2658 }
2659 }
2660 deleteAllAttributes(name);
2661 H5Ldelete(current_group,name.c_str(),H5P_DEFAULT);
2662 }
2663 }
2664
2665 push(name);
2666 for(unsigned int i=0; i<Nd; i++){
2667 std::stringstream stream;
2668 stream << i;
2669 std::string dname(stream.str());
2670 write(dname,field[i],HDF5Base::ate);
2671 writeAttribute(dname,".kind","LatticeColorMatrix");
2672 }
2673 pop();
2674 if(profile) swatch_complete.stop();
2675 QDPIO::cout<< "done!" << std::endl;
2676
2677 if(profile){
2678 QDPIO::cout << "HDF5-I/O statistics for Qlua-write: " << std::endl;
2679 QDPIO::cout << "\t total: " << swatch_complete.getTimeInSeconds() << " s." << std::endl;
2680 QDPIO::cout << "\t MB written: " << Layout::vol()*field.size()*sizeof(ColorMatrixD3)/1024/1024 << std::endl;
2681 }
2682 }
2683}
2684#endif
2685
HDF5(const long int &stripesizee=-1, const long int &maxalign=0)
OScalar< PScalar< PScalar< RComplex< REAL32 > > > > ComplexF
OScalar< PScalar< PColorMatrix< RComplex< REAL32 >, 3 > > > ColorMatrixF3
OLattice< PSpinMatrix< PColorMatrix< RComplex< REAL64 >, 3 >, 4 > > LatticeDiracPropagatorD3
OScalar< PScalar< PColorMatrix< RComplex< REAL64 >, 3 > > > ColorMatrixD3
OScalar< PSpinMatrix< PColorMatrix< RComplex< REAL64 >, 3 >, Ns > > PropagatorD3
OScalar< PScalar< PScalar< RComplex< REAL64 > > > > ComplexD
OScalar< PSpinMatrix< PColorMatrix< RComplex< REAL32 >, 3 >, Ns > > PropagatorF3
OScalar< PSpinMatrix< PColorMatrix< RComplex< REAL64 >, 3 >, 4 > > DiracPropagatorD3
OLattice< PSpinMatrix< PColorMatrix< RComplex< REAL32 >, 3 >, 4 > > LatticeDiracPropagatorF3
OLattice< PScalar< PColorMatrix< RComplex< REAL32 >, 3 > > > LatticeColorMatrixF3
OLattice< PScalar< PColorMatrix< RComplex< REAL64 >, 3 > > > LatticeColorMatrixD3
OLattice< PSpinMatrix< PColorMatrix< RComplex< REAL32 >, 3 >, Ns > > LatticePropagatorF3
OLattice< PSpinMatrix< PColorMatrix< RComplex< REAL64 >, 3 >, Ns > > LatticePropagatorD3
OScalar< PSpinMatrix< PColorMatrix< RComplex< REAL32 >, 3 >, 4 > > DiracPropagatorF3
void write(BinaryWriter &bin, const std::string &output)
Definition qdp_io.cc:1204
void read(BinaryReader &bin, std::string &input, size_t maxBytes)
Definition qdp_io.cc:778
const int Nd
Definition qdp_params.h:24
void close(QDPFileReader &qsw)
Close a QDPFileReader.
Definition qdp_qdpio.cc:207
Yet another random number generator.
bool get_global(ctype &global, const ctype &local)
void push(XMLWriter &xml, const std::string &s)
Push a group name.
Definition qdp_xmlio.cc:749
void pop(XMLWriter &xml)
Pop a group name.
Definition qdp_xmlio.cc:752
multi1d< int > crtesn(int ipos, const multi1d< int > &latt_size)
Decompose a lexicographic site into coordinates.
Primary include file for QDP.
#define local
Definition qdp_crc32.cc:51
unsigned long long ullong
Definition qdp_hdf5.h:24
double REAL64
float REAL32
REAL32 REAL