EpetraExt Package Browser (Single Doxygen Collection)  Development
EpetraExt_HDF5.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // EpetraExt: Epetra Extended - Linear Algebra Services Package
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #include "EpetraExt_ConfigDefs.h"
45 
46 
47 #ifdef HAVE_EPETRAEXT_HDF5
48 
49 #include "EpetraExt_HDF5.h"
50 #ifdef HAVE_MPI
51 # include "mpi.h"
52 # include <H5FDmpio.h>
53 # include "Epetra_MpiComm.h"
54 #endif // HAVE_MPI
55 
56 // The user could have passed in an Epetra_Comm that is either an
57 // Epetra_MpiComm or an Epetra_SerialComm. The latter could happen
58 // whether or not we built Trilinos with MPI. Thus, we need to
59 // include this header regardless.
60 #include "Epetra_SerialComm.h"
61 
62 #include "Teuchos_ParameterList.hpp"
63 #include "Teuchos_RCP.hpp"
64 #include "Epetra_Map.h"
65 #include "Epetra_BlockMap.h"
66 #include "Epetra_CrsGraph.h"
67 #include "Epetra_FECrsGraph.h"
68 #include "Epetra_RowMatrix.h"
69 #include "Epetra_CrsMatrix.h"
70 #include "Epetra_FECrsMatrix.h"
71 #include "Epetra_IntVector.h"
72 #include "Epetra_MultiVector.h"
73 #include "Epetra_Import.h"
74 #include "EpetraExt_Exception.h"
75 #include "EpetraExt_Utils.h"
76 #include "EpetraExt_DistArray.h"
77 
78 #define CHECK_HID(hid_t) \
79  { if (hid_t < 0) \
80  throw(EpetraExt::Exception(__FILE__, __LINE__, \
81  "hid_t is negative")); }
82 
83 #define CHECK_STATUS(status) \
84  { if (status < 0) \
85  throw(EpetraExt::Exception(__FILE__, __LINE__, \
86  "function H5Giterater returned a negative value")); }
87 
88 // ==========================================================================
89 // data container and iterators to find a dataset with a given name
91 {
92  std::string name;
93  bool found;
94 };
95 
96 static herr_t FindDataset(hid_t loc_id, const char *name, void *opdata)
97 {
98  std::string& token = ((FindDataset_t*)opdata)->name;
99  if (token == name)
100  ((FindDataset_t*)opdata)->found = true;
101 
102  return(0);
103 }
104 
105 // ==========================================================================
106 // This function copied from Roman Geus' FEMAXX code
107 static void WriteParameterListRecursive(const Teuchos::ParameterList& params,
108  hid_t group_id)
109 {
110  Teuchos::ParameterList::ConstIterator it = params.begin();
111  for (; it != params.end(); ++ it)
112  {
113  std::string key(params.name(it));
114  if (params.isSublist(key))
115  {
116  // Sublist
117 
118  // Create subgroup for sublist.
119  hid_t child_group_id = H5Gcreate(group_id, key.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
120  WriteParameterListRecursive(params.sublist(key), child_group_id);
121  H5Gclose(child_group_id);
122  }
123  else
124  {
125  //
126  // Regular parameter
127  //
128 
129  // Create dataspace/dataset.
130  herr_t status;
131  hsize_t one = 1;
132  hid_t dataspace_id, dataset_id;
133  bool found = false; // to avoid a compiler error on MAC OS X GCC 4.0.0
134 
135  // Write the dataset.
136  if (params.isType<std::string>(key))
137  {
138  std::string value = params.get<std::string>(key);
139  hsize_t len = value.size() + 1;
140  dataspace_id = H5Screate_simple(1, &len, NULL);
141  dataset_id = H5Dcreate(group_id, key.c_str(), H5T_C_S1, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
142  status = H5Dwrite(dataset_id, H5T_C_S1, H5S_ALL, H5S_ALL,
143  H5P_DEFAULT, value.c_str());
144  CHECK_STATUS(status);
145  status = H5Dclose(dataset_id);
146  CHECK_STATUS(status);
147  status = H5Sclose(dataspace_id);
148  CHECK_STATUS(status);
149  found = true;
150  }
151 
152  if (params.isType<bool>(key))
153  {
154  // Use H5T_NATIVE_USHORT to store a bool value
155  unsigned short value = params.get<bool>(key) ? 1 : 0;
156  dataspace_id = H5Screate_simple(1, &one, NULL);
157  dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_USHORT, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
158  status = H5Dwrite(dataset_id, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL,
159  H5P_DEFAULT, &value);
160  CHECK_STATUS(status);
161  status = H5Dclose(dataset_id);
162  CHECK_STATUS(status);
163  status = H5Sclose(dataspace_id);
164  CHECK_STATUS(status);
165  found = true;
166  }
167 
168  if (params.isType<int>(key))
169  {
170  int value = params.get<int>(key);
171  dataspace_id = H5Screate_simple(1, &one, NULL);
172  dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_INT, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
173  status = H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
174  H5P_DEFAULT, &value);
175  CHECK_STATUS(status);
176  status = H5Dclose(dataset_id);
177  CHECK_STATUS(status);
178  status = H5Sclose(dataspace_id);
179  CHECK_STATUS(status);
180  found = true;
181  }
182 
183  if (params.isType<double>(key))
184  {
185  double value = params.get<double>(key);
186  dataspace_id = H5Screate_simple(1, &one, NULL);
187  dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
188  status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
189  H5P_DEFAULT, &value);
190  CHECK_STATUS(status);
191  status = H5Dclose(dataset_id);
192  CHECK_STATUS(status);
193  status = H5Sclose(dataspace_id);
194  CHECK_STATUS(status);
195  found = true;
196  }
197 
198  if (!found)
199  {
200  throw(EpetraExt::Exception(__FILE__, __LINE__,
201  "type for parameter " + key + " not supported"));
202  }
203  }
204  }
205 }
206 
207 // ==========================================================================
208 // Recursive Operator function called by H5Giterate for each entity in group.
209 // This function copied from Roman Geus' FEMAXX code
210 static herr_t f_operator(hid_t loc_id, const char *name, void *opdata)
211 {
212  H5G_stat_t statbuf;
213  hid_t dataset_id, space_id, type_id;
214  Teuchos::ParameterList* sublist;
215  Teuchos::ParameterList* params = (Teuchos::ParameterList*)opdata;
216  /*
217  * Get type of the object and display its name and type.
218  * The name of the object is passed to this function by
219  * the Library. Some magic :-)
220  */
221  H5Gget_objinfo(loc_id, name, 0, &statbuf);
222  if (strcmp(name, "__type__") == 0)
223  return(0); // skip list identifier
224 
225  switch (statbuf.type) {
226  case H5G_GROUP:
227  sublist = &params->sublist(name);
228  H5Giterate(loc_id, name , NULL, f_operator, sublist);
229  break;
230  case H5G_DATASET:
231  hsize_t len;
232  dataset_id = H5Dopen(loc_id, name, H5P_DEFAULT);
233  space_id = H5Dget_space(dataset_id);
234  if (H5Sget_simple_extent_ndims(space_id) != 1)
235  throw(EpetraExt::Exception(__FILE__, __LINE__,
236  "dimensionality of parameters must be 1."));
237  H5Sget_simple_extent_dims(space_id, &len, NULL);
238  type_id = H5Dget_type(dataset_id);
239  if (H5Tequal(type_id, H5T_NATIVE_DOUBLE) > 0) {
240  double value;
241  H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
242  params->set(name, value);
243  } else if (H5Tequal(type_id, H5T_NATIVE_INT) > 0) {
244  int value;
245  H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
246  params->set(name, value);
247  } else if (H5Tequal(type_id, H5T_C_S1) > 0) {
248  char* buf = new char[len];
249  H5Dread(dataset_id, H5T_C_S1, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
250  params->set(name, std::string(buf));
251  delete[] buf;
252  } else if (H5Tequal(type_id, H5T_NATIVE_USHORT) > 0) {
253  unsigned short value;
254  H5Dread(dataset_id, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
255  params->set(name, value != 0 ? true : false);
256  } else {
257  throw(EpetraExt::Exception(__FILE__, __LINE__,
258  "unsupported datatype")); // FIXME
259  }
260  H5Tclose(type_id);
261  H5Sclose(space_id);
262  H5Dclose(dataset_id);
263  break;
264  default:
265  throw(EpetraExt::Exception(__FILE__, __LINE__,
266  "unsupported datatype")); // FIXME
267  }
268  return 0;
269 }
270 
271 // ==========================================================================
272 EpetraExt::HDF5::HDF5(const Epetra_Comm& Comm) :
273  Comm_(Comm),
274  IsOpen_(false)
275 {}
276 
277 // ==========================================================================
278 void EpetraExt::HDF5::Create(const std::string FileName)
279 {
280  if (IsOpen())
281  throw(Exception(__FILE__, __LINE__,
282  "an HDF5 is already open, first close the current one",
283  "using method Close(), then open/create a new one"));
284 
285  FileName_ = FileName;
286 
287  // Set up file access property list with parallel I/O access
288  plist_id_ = H5Pcreate(H5P_FILE_ACCESS);
289 #ifdef HAVE_MPI
290  {
291  // Tell HDF5 what MPI communicator to use for parallel file access
292  // for the above property list.
293  //
294  // HAVE_MPI is defined, so we know that Trilinos was built with
295  // MPI. However, we don't know whether Comm_ wraps an MPI
296  // communicator. Comm_ could very well be a serial communicator.
297  // Try a dynamic cast to Epetra_MpiComm to find out.
298  MPI_Comm mpiComm = MPI_COMM_NULL; // Hopefully not for long
299 
300  // Is Comm_ an Epetra_MpiComm?
301  const Epetra_MpiComm* mpiWrapper =
302  dynamic_cast<const Epetra_MpiComm*> (&Comm_);
303  if (mpiWrapper != NULL) {
304  mpiComm = mpiWrapper->Comm();
305  }
306  else {
307  // Is Comm_ an Epetra_SerialComm?
308  const Epetra_SerialComm* serialWrapper =
309  dynamic_cast<const Epetra_SerialComm*> (&Comm_);
310 
311  if (serialWrapper != NULL) {
312  // Comm_ is an Epetra_SerialComm. This means that even though
313  // Trilinos was built with MPI, the user who instantiated the
314  // HDF5 class wants only the calling process to access HDF5.
315  // The right communicator to use in that case is
316  // MPI_COMM_SELF.
317  mpiComm = MPI_COMM_SELF;
318  } else {
319  // Comm_ must be some other subclass of Epetra_Comm.
320  // We don't know how to get an MPI communicator out of it.
321  const char* const errMsg = "EpetraExt::HDF5::Create: This HDF5 object "
322  "was created with an Epetra_Comm instance which is neither an "
323  "Epetra_MpiComm nor a Epetra_SerialComm. As a result, we do not "
324  "know how to get an MPI communicator from it. Our HDF5 class only "
325  "understands Epetra_Comm objects which are instances of one of these "
326  "two subclasses.";
327  throw EpetraExt::Exception (__FILE__, __LINE__, errMsg);
328  }
329  }
330 
331  // By this point, mpiComm should be something other than
332  // MPI_COMM_NULL. Otherwise, Comm_ wraps MPI_COMM_NULL.
333  if (mpiComm == MPI_COMM_NULL) {
334  const char* const errMsg = "EpetraExt::HDF5::Create: The Epetra_Comm "
335  "object with which this HDF5 instance was created wraps MPI_COMM_NULL, "
336  "which is an invalid MPI communicator. HDF5 requires a valid MPI "
337  "communicator.";
338  throw EpetraExt::Exception (__FILE__, __LINE__, errMsg);
339  }
340 
341  // Tell HDF5 what MPI communicator to use for parallel file access
342  // for the above property list. For details, see e.g.,
343  //
344  // http://www.hdfgroup.org/HDF5/doc/UG/08_TheFile.html
345  //
346  // [last accessed 06 Oct 2011]
347  H5Pset_fapl_mpio(plist_id_, mpiComm, MPI_INFO_NULL);
348  }
349 #endif
350 
351 #if 0
352  unsigned int boh = H5Z_FILTER_MAX;
353  H5Pset_filter(plist_id_, H5Z_FILTER_DEFLATE, H5Z_FILTER_MAX, 0, &boh);
354 #endif
355 
356  // create the file collectively and release property list identifier.
357  file_id_ = H5Fcreate(FileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT,
358  plist_id_);
359  H5Pclose(plist_id_);
360 
361  IsOpen_ = true;
362 }
363 
364 // ==========================================================================
365 void EpetraExt::HDF5::Open(const std::string FileName, int AccessType)
366 {
367  if (IsOpen())
368  throw(Exception(__FILE__, __LINE__,
369  "an HDF5 is already open, first close the current one",
370  "using method Close(), then open/create a new one"));
371 
372  FileName_ = FileName;
373 
374  // Set up file access property list with parallel I/O access
375  plist_id_ = H5Pcreate(H5P_FILE_ACCESS);
376 
377 #ifdef HAVE_MPI
378  // Create property list for collective dataset write.
379  const Epetra_MpiComm* MpiComm ( dynamic_cast<const Epetra_MpiComm*> (&Comm_) );
380 
381  if (MpiComm == 0)
382  H5Pset_fapl_mpio(plist_id_, MPI_COMM_WORLD, MPI_INFO_NULL);
383  else
384  H5Pset_fapl_mpio(plist_id_, MpiComm->Comm(), MPI_INFO_NULL);
385 #endif
386 
387  // create the file collectively and release property list identifier.
388  file_id_ = H5Fopen(FileName.c_str(), AccessType, plist_id_);
389  H5Pclose(plist_id_);
390 
391  IsOpen_ = true;
392 }
393 
394 // ==========================================================================
395 bool EpetraExt::HDF5::IsContained(std::string Name)
396 {
397  if (!IsOpen())
398  throw(Exception(__FILE__, __LINE__, "no file open yet"));
399 
400  FindDataset_t data;
401  data.name = Name;
402  data.found = false;
403 
404  //int idx_f =
405  H5Giterate(file_id_, "/", NULL, FindDataset, (void*)&data);
406 
407  return(data.found);
408 }
409 
410 // ============================ //
411 // Epetra_BlockMap / Epetra_Map //
412 // ============================ //
413 
414 // ==========================================================================
415 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_BlockMap& BlockMap)
416 {
417  int NumMyPoints = BlockMap.NumMyPoints();
418  int NumMyElements = BlockMap.NumMyElements();
419  int NumGlobalElements = BlockMap.NumGlobalElements();
420  int NumGlobalPoints = BlockMap.NumGlobalPoints();
421  int* MyGlobalElements = BlockMap.MyGlobalElements();
422  int* ElementSizeList = BlockMap.ElementSizeList();
423 
424  std::vector<int> NumMyElements_v(Comm_.NumProc());
425  Comm_.GatherAll(&NumMyElements, &NumMyElements_v[0], 1);
426 
427  std::vector<int> NumMyPoints_v(Comm_.NumProc());
428  Comm_.GatherAll(&NumMyPoints, &NumMyPoints_v[0], 1);
429 
430  Write(Name, "MyGlobalElements", NumMyElements, NumGlobalElements,
431  H5T_NATIVE_INT, MyGlobalElements);
432  Write(Name, "ElementSizeList", NumMyElements, NumGlobalElements,
433  H5T_NATIVE_INT, ElementSizeList);
434  Write(Name, "NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
435 
436  // need to know how many processors currently host this map
437  Write(Name, "NumProc", Comm_.NumProc());
438  // write few more data about this map
439  Write(Name, "NumGlobalPoints", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalPoints);
440  Write(Name, "NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalElements);
441  Write(Name, "IndexBase", BlockMap.IndexBase());
442  Write(Name, "__type__", "Epetra_BlockMap");
443 }
444 
445 // ==========================================================================
446 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_BlockMap*& BlockMap)
447 {
448  int NumGlobalElements, NumGlobalPoints, IndexBase, NumProc;
449 
450  ReadBlockMapProperties(GroupName, NumGlobalElements, NumGlobalPoints,
451  IndexBase, NumProc);
452 
453  std::vector<int> NumMyPoints_v(Comm_.NumProc());
454  std::vector<int> NumMyElements_v(Comm_.NumProc());
455 
456  Read(GroupName, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
457  Read(GroupName, "NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
458  int NumMyElements = NumMyElements_v[Comm_.MyPID()];
459 // int NumMyPoints = NumMyPoints_v[Comm_.MyPID()];
460 
461  if (NumProc != Comm_.NumProc())
462  throw(Exception(__FILE__, __LINE__,
463  "requested map not compatible with current number of",
464  "processors, " + toString(Comm_.NumProc()) +
465  " vs. " + toString(NumProc)));
466 
467  std::vector<int> MyGlobalElements(NumMyElements);
468  std::vector<int> ElementSizeList(NumMyElements);
469 
470  Read(GroupName, "MyGlobalElements", NumMyElements, NumGlobalElements,
471  H5T_NATIVE_INT, &MyGlobalElements[0]);
472 
473  Read(GroupName, "ElementSizeList", NumMyElements, NumGlobalElements,
474  H5T_NATIVE_INT, &ElementSizeList[0]);
475 
476  BlockMap = new Epetra_BlockMap(NumGlobalElements, NumMyElements,
477  &MyGlobalElements[0], &ElementSizeList[0],
478  IndexBase, Comm_);
479 }
480 
481 // ==========================================================================
482 void EpetraExt::HDF5::ReadBlockMapProperties(const std::string& GroupName,
483  int& NumGlobalElements,
484  int& NumGlobalPoints,
485  int& IndexBase,
486  int& NumProc)
487 {
488  if (!IsContained(GroupName))
489  throw(Exception(__FILE__, __LINE__,
490  "requested group " + GroupName + " not found"));
491 
492  std::string Label;
493  Read(GroupName, "__type__", Label);
494 
495  if (Label != "Epetra_BlockMap")
496  throw(Exception(__FILE__, __LINE__,
497  "requested group " + GroupName + " is not an Epetra_BlockMap",
498  "__type__ = " + Label));
499 
500  Read(GroupName, "NumGlobalElements", NumGlobalElements);
501  Read(GroupName, "NumGlobalPoints", NumGlobalPoints);
502  Read(GroupName, "IndexBase", IndexBase);
503  Read(GroupName, "NumProc", NumProc);
504 }
505 
506 // ==========================================================================
507 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_Map& Map)
508 {
509  int MySize = Map.NumMyElements();
510  int GlobalSize = Map.NumGlobalElements();
511  int* MyGlobalElements = Map.MyGlobalElements();
512 
513  std::vector<int> NumMyElements(Comm_.NumProc());
514  Comm_.GatherAll(&MySize, &NumMyElements[0], 1);
515 
516  Write(Name, "MyGlobalElements", MySize, GlobalSize,
517  H5T_NATIVE_INT, MyGlobalElements);
518  Write(Name, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements[0]);
519  Write(Name, "NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &GlobalSize);
520  Write(Name, "NumProc", Comm_.NumProc());
521  Write(Name, "IndexBase", Map.IndexBase());
522  Write(Name, "__type__", "Epetra_Map");
523 }
524 
525 // ==========================================================================
526 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_Map*& Map)
527 {
528  int NumGlobalElements, IndexBase, NumProc;
529 
530  ReadMapProperties(GroupName, NumGlobalElements, IndexBase, NumProc);
531 
532  std::vector<int> NumMyElements_v(Comm_.NumProc());
533 
534  Read(GroupName, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
535  int NumMyElements = NumMyElements_v[Comm_.MyPID()];
536 
537  if (NumProc != Comm_.NumProc())
538  throw(Exception(__FILE__, __LINE__,
539  "requested map not compatible with current number of",
540  "processors, " + toString(Comm_.NumProc()) +
541  " vs. " + toString(NumProc)));
542 
543  std::vector<int> MyGlobalElements(NumMyElements);
544 
545  Read(GroupName, "MyGlobalElements", NumMyElements, NumGlobalElements,
546  H5T_NATIVE_INT, &MyGlobalElements[0]);
547 
548  Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], IndexBase, Comm_);
549 }
550 
551 // ==========================================================================
552 void EpetraExt::HDF5::ReadMapProperties(const std::string& GroupName,
553  int& NumGlobalElements,
554  int& IndexBase,
555  int& NumProc)
556 {
557  if (!IsContained(GroupName))
558  throw(Exception(__FILE__, __LINE__,
559  "requested group " + GroupName + " not found"));
560 
561  std::string Label;
562  Read(GroupName, "__type__", Label);
563 
564  if (Label != "Epetra_Map")
565  throw(Exception(__FILE__, __LINE__,
566  "requested group " + GroupName + " is not an Epetra_Map",
567  "__type__ = " + Label));
568 
569  Read(GroupName, "NumGlobalElements", NumGlobalElements);
570  Read(GroupName, "IndexBase", IndexBase);
571  Read(GroupName, "NumProc", NumProc);
572 }
573 
574 // ================ //
575 // Epetra_IntVector //
576 // ================ //
577 
578 // ==========================================================================
579 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_IntVector& x)
580 {
581  if (x.Map().LinearMap())
582  {
583  Write(Name, "GlobalLength", x.GlobalLength());
584  Write(Name, "Values", x.Map().NumMyElements(), x.Map().NumGlobalElements(),
585  H5T_NATIVE_INT, x.Values());
586  }
587  else
588  {
589  // need to build a linear map first, the import data, then
590  // finally write them
591  const Epetra_BlockMap& OriginalMap = x.Map();
592  Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
593  Epetra_Import Importer(LinearMap, OriginalMap);
594 
595  Epetra_IntVector LinearX(LinearMap);
596  LinearX.Import(x, Importer, Insert);
597 
598  Write(Name, "GlobalLength", x.GlobalLength());
599  Write(Name, "Values", LinearMap.NumMyElements(), LinearMap.NumGlobalElements(),
600  H5T_NATIVE_INT, LinearX.Values());
601  }
602  Write(Name, "__type__", "Epetra_IntVector");
603 }
604 
605 // ==========================================================================
606 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_IntVector*& X)
607 {
608  // gets the length of the std::vector
609  int GlobalLength;
610  ReadIntVectorProperties(GroupName, GlobalLength);
611 
612  // creates a new linear map and use it to read data
613  Epetra_Map LinearMap(GlobalLength, 0, Comm_);
614  X = new Epetra_IntVector(LinearMap);
615 
616  Read(GroupName, "Values", LinearMap.NumMyElements(),
617  LinearMap.NumGlobalElements(), H5T_NATIVE_INT, X->Values());
618 }
619 
620 // ==========================================================================
621 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
622  Epetra_IntVector*& X)
623 {
624  // gets the length of the std::vector
625  int GlobalLength;
626  ReadIntVectorProperties(GroupName, GlobalLength);
627 
628  if (Map.LinearMap())
629  {
630  X = new Epetra_IntVector(Map);
631  // simply read stuff and go home
632  Read(GroupName, "Values", Map.NumMyElements(), Map.NumGlobalElements(),
633  H5T_NATIVE_INT, X->Values());
634 
635  }
636  else
637  {
638  // we need to first create a linear map, read the std::vector,
639  // then import it to the actual nonlinear map
640  Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
641  Epetra_IntVector LinearX(LinearMap);
642 
643  Read(GroupName, "Values", LinearMap.NumMyElements(),
644  LinearMap.NumGlobalElements(),
645  H5T_NATIVE_INT, LinearX.Values());
646 
647  Epetra_Import Importer(Map, LinearMap);
648  X = new Epetra_IntVector(Map);
649  X->Import(LinearX, Importer, Insert);
650  }
651 }
652 
653 // ==========================================================================
654 void EpetraExt::HDF5::ReadIntVectorProperties(const std::string& GroupName,
655  int& GlobalLength)
656 {
657  if (!IsContained(GroupName))
658  throw(Exception(__FILE__, __LINE__,
659  "requested group " + GroupName + " not found"));
660 
661  std::string Label;
662  Read(GroupName, "__type__", Label);
663 
664  if (Label != "Epetra_IntVector")
665  throw(Exception(__FILE__, __LINE__,
666  "requested group " + GroupName + " is not an Epetra_IntVector",
667  "__type__ = " + Label));
668 
669  Read(GroupName, "GlobalLength", GlobalLength);
670 }
671 
672 // =============== //
673 // Epetra_CrsGraph //
674 // =============== //
675 
676 // ==========================================================================
677 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_CrsGraph& Graph)
678 {
679  if (!Graph.Filled())
680  throw(Exception(__FILE__, __LINE__,
681  "input Epetra_CrsGraph is not FillComplete()'d"));
682 
683  // like RowMatrix, only without values
684  int MySize = Graph.NumMyNonzeros();
685  int GlobalSize = Graph.NumGlobalNonzeros();
686  std::vector<int> ROW(MySize);
687  std::vector<int> COL(MySize);
688 
689  int count = 0;
690  int* RowIndices;
691  int NumEntries;
692 
693  for (int i = 0; i < Graph.NumMyRows(); ++i)
694  {
695  Graph.ExtractMyRowView(i, NumEntries, RowIndices);
696  for (int j = 0; j < NumEntries; ++j)
697  {
698  ROW[count] = Graph.GRID(i);
699  COL[count] = Graph.GCID(RowIndices[j]);
700  ++count;
701  }
702  }
703 
704  Write(Name, "ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
705  Write(Name, "COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
706  Write(Name, "MaxNumIndices", Graph.MaxNumIndices());
707  Write(Name, "NumGlobalRows", Graph.NumGlobalRows());
708  Write(Name, "NumGlobalCols", Graph.NumGlobalCols());
709  Write(Name, "NumGlobalNonzeros", Graph.NumGlobalNonzeros());
710  Write(Name, "NumGlobalDiagonals", Graph.NumGlobalDiagonals());
711  Write(Name, "__type__", "Epetra_CrsGraph");
712 }
713 
714 // ==========================================================================
715 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_CrsGraph*& Graph)
716 {
717  int NumGlobalRows, NumGlobalCols;
718  int NumGlobalNonzeros, NumGlobalDiagonals, MaxNumIndices;
719 
720  ReadCrsGraphProperties(GroupName, NumGlobalRows,
721  NumGlobalCols, NumGlobalNonzeros,
722  NumGlobalDiagonals, MaxNumIndices);
723 
724  Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
725  Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
726 
727  Read(GroupName, DomainMap, RangeMap, Graph);
728 }
729 
730 // ==========================================================================
731 void EpetraExt::HDF5::ReadCrsGraphProperties(const std::string& GroupName,
732  int& NumGlobalRows,
733  int& NumGlobalCols,
734  int& NumGlobalNonzeros,
735  int& NumGlobalDiagonals,
736  int& MaxNumIndices)
737 {
738  if (!IsContained(GroupName))
739  throw(Exception(__FILE__, __LINE__,
740  "requested group " + GroupName + " not found"));
741 
742  std::string Label;
743  Read(GroupName, "__type__", Label);
744 
745  if (Label != "Epetra_CrsGraph")
746  throw(Exception(__FILE__, __LINE__,
747  "requested group " + GroupName + " is not an Epetra_CrsGraph",
748  "__type__ = " + Label));
749 
750  Read(GroupName, "NumGlobalRows", NumGlobalRows);
751  Read(GroupName, "NumGlobalCols", NumGlobalCols);
752  Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
753  Read(GroupName, "NumGlobalDiagonals", NumGlobalDiagonals);
754  Read(GroupName, "MaxNumIndices", MaxNumIndices);
755 }
756 
757 // ==========================================================================
758 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& DomainMap,
759  const Epetra_Map& RangeMap, Epetra_CrsGraph*& Graph)
760 {
761  if (!IsContained(GroupName))
762  throw(Exception(__FILE__, __LINE__,
763  "requested group " + GroupName + " not found in database"));
764 
765  std::string Label;
766  Read(GroupName, "__type__", Label);
767 
768  if (Label != "Epetra_CrsGraph")
769  throw(Exception(__FILE__, __LINE__,
770  "requested group " + GroupName + " is not an Epetra_CrsGraph",
771  "__type__ = " + Label));
772 
773  int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
774  Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
775  Read(GroupName, "NumGlobalRows", NumGlobalRows);
776  Read(GroupName, "NumGlobalCols", NumGlobalCols);
777 
778  // linear distribution for nonzeros
779  int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
780  if (Comm_.MyPID() == 0)
781  NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
782 
783  std::vector<int> ROW(NumMyNonzeros);
784  Read(GroupName, "ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
785 
786  std::vector<int> COL(NumMyNonzeros);
787  Read(GroupName, "COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
788 
789  Epetra_FECrsGraph* Graph2 = new Epetra_FECrsGraph(Copy, DomainMap, 0);
790 
791  for (int i = 0; i < NumMyNonzeros; )
792  {
793  int count = 1;
794  while (ROW[i + count] == ROW[i])
795  ++count;
796 
797  Graph2->InsertGlobalIndices(1, &ROW[i], count, &COL[i]);
798  i += count;
799  }
800 
801  Graph2->FillComplete(DomainMap, RangeMap);
802 
803  Graph = Graph2;
804 }
805 
806 // ================ //
807 // Epetra_RowMatrix //
808 // ================ //
809 
810 // ==========================================================================
811 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_RowMatrix& Matrix)
812 {
813  if (!Matrix.Filled())
814  throw(Exception(__FILE__, __LINE__,
815  "input Epetra_RowMatrix is not FillComplete()'d"));
816 
817  int MySize = Matrix.NumMyNonzeros();
818  int GlobalSize = Matrix.NumGlobalNonzeros();
819  std::vector<int> ROW(MySize);
820  std::vector<int> COL(MySize);
821  std::vector<double> VAL(MySize);
822 
823  int count = 0;
824  int Length = Matrix.MaxNumEntries();
825  std::vector<int> Indices(Length);
826  std::vector<double> Values(Length);
827  int NumEntries;
828 
829  for (int i = 0; i < Matrix.NumMyRows(); ++i)
830  {
831  Matrix.ExtractMyRowCopy(i, Length, NumEntries, &Values[0], &Indices[0]);
832  for (int j = 0; j < NumEntries; ++j)
833  {
834  ROW[count] = Matrix.RowMatrixRowMap().GID(i);
835  COL[count] = Matrix.RowMatrixColMap().GID(Indices[j]);
836  VAL[count] = Values[j];
837  ++count;
838  }
839  }
840 
841  Write(Name, "ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
842  Write(Name, "COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
843  Write(Name, "VAL", MySize, GlobalSize, H5T_NATIVE_DOUBLE, &VAL[0]);
844  Write(Name, "NumGlobalRows", Matrix.NumGlobalRows());
845  Write(Name, "NumGlobalCols", Matrix.NumGlobalCols());
846  Write(Name, "NumGlobalNonzeros", Matrix.NumGlobalNonzeros());
847  Write(Name, "NumGlobalDiagonals", Matrix.NumGlobalDiagonals());
848  Write(Name, "MaxNumEntries", Matrix.MaxNumEntries());
849  Write(Name, "NormOne", Matrix.NormOne());
850  Write(Name, "NormInf", Matrix.NormInf());
851  Write(Name, "__type__", "Epetra_RowMatrix");
852 }
853 
854 // ==========================================================================
855 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_CrsMatrix*& A)
856 {
857  int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
858  int NumGlobalDiagonals, MaxNumEntries;
859  double NormOne, NormInf;
860 
861  ReadCrsMatrixProperties(GroupName, NumGlobalRows, NumGlobalCols,
862  NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
863  NormOne, NormInf);
864 
865  // build simple linear maps for domain and range space
866  Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
867  Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
868 
869  Read(GroupName, DomainMap, RangeMap, A);
870 }
871 
872 // ==========================================================================
873 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& DomainMap,
874  const Epetra_Map& RangeMap, Epetra_CrsMatrix*& A)
875 {
876  int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
877  int NumGlobalDiagonals, MaxNumEntries;
878  double NormOne, NormInf;
879 
880  ReadCrsMatrixProperties(GroupName, NumGlobalRows, NumGlobalCols,
881  NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
882  NormOne, NormInf);
883 
884  int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
885  if (Comm_.MyPID() == 0)
886  NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
887 
888  std::vector<int> ROW(NumMyNonzeros);
889  Read(GroupName, "ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
890 
891  std::vector<int> COL(NumMyNonzeros);
892  Read(GroupName, "COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
893 
894  std::vector<double> VAL(NumMyNonzeros);
895  Read(GroupName, "VAL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_DOUBLE, &VAL[0]);
896 
897  Epetra_FECrsMatrix* A2 = new Epetra_FECrsMatrix(Copy, DomainMap, 0);
898 
899  for (int i = 0; i < NumMyNonzeros; )
900  {
901  int count = 1;
902  while (ROW[i + count] == ROW[i])
903  ++count;
904 
905  A2->InsertGlobalValues(1, &ROW[i], count, &COL[i], &VAL[i]);
906  i += count;
907  }
908 
909  A2->GlobalAssemble(DomainMap, RangeMap);
910 
911  A = A2;
912 }
913 
914 // ==========================================================================
915 void EpetraExt::HDF5::ReadCrsMatrixProperties(const std::string& GroupName,
916  int& NumGlobalRows,
917  int& NumGlobalCols,
918  int& NumGlobalNonzeros,
919  int& NumGlobalDiagonals,
920  int& MaxNumEntries,
921  double& NormOne,
922  double& NormInf)
923 {
924  if (!IsContained(GroupName))
925  throw(Exception(__FILE__, __LINE__,
926  "requested group " + GroupName + " not found"));
927 
928  std::string Label;
929  Read(GroupName, "__type__", Label);
930 
931  if (Label != "Epetra_RowMatrix")
932  throw(Exception(__FILE__, __LINE__,
933  "requested group " + GroupName + " is not an Epetra_RowMatrix",
934  "__type__ = " + Label));
935 
936  Read(GroupName, "NumGlobalRows", NumGlobalRows);
937  Read(GroupName, "NumGlobalCols", NumGlobalCols);
938  Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
939  Read(GroupName, "NumGlobalDiagonals", NumGlobalDiagonals);
940  Read(GroupName, "MaxNumEntries", MaxNumEntries);
941  Read(GroupName, "NormOne", NormOne);
942  Read(GroupName, "NormInf", NormInf);
943 }
944 
945 // ============= //
946 // ParameterList //
947 // ============= //
948 
949 // ==========================================================================
950 void EpetraExt::HDF5::Write(const std::string& GroupName, const Teuchos::ParameterList& params)
951 {
952  if (!IsOpen())
953  throw(Exception(__FILE__, __LINE__, "no file open yet"));
954 
955  hid_t group_id = H5Gcreate(file_id_, GroupName.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
956 
957  // Iterate through parameter list
958  WriteParameterListRecursive(params, group_id);
959 
960  // Finalize hdf5 file
961  status = H5Gclose(group_id);
962  CHECK_STATUS(status);
963 
964  Write(GroupName, "__type__", "Teuchos::ParameterList");
965 }
966 
967 // ==========================================================================
968 void EpetraExt::HDF5::Read(const std::string& GroupName, Teuchos::ParameterList& params)
969 {
970  if (!IsOpen())
971  throw(Exception(__FILE__, __LINE__, "no file open yet"));
972 
973  std::string Label;
974  Read(GroupName, "__type__", Label);
975 
976  if (Label != "Teuchos::ParameterList")
977  throw(Exception(__FILE__, __LINE__,
978  "requested group " + GroupName + " is not a Teuchos::ParameterList",
979  "__type__ = " + Label));
980 
981  // Open hdf5 file
982  hid_t group_id; /* identifiers */
983  herr_t status;
984 
985  // open group in the root group using absolute name.
986  group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
987  CHECK_HID(group_id);
988 
989  // Iterate through parameter list
990  std::string xxx = "/" + GroupName;
991  status = H5Giterate(group_id, xxx.c_str() , NULL, f_operator, &params);
992  CHECK_STATUS(status);
993 
994  // Finalize hdf5 file
995  status = H5Gclose(group_id);
996  CHECK_STATUS(status);
997 }
998 
999 // ================== //
1000 // Epetra_MultiVector //
1001 // ================== //
1002 
1003 // ==========================================================================
1004 void EpetraExt::HDF5::Write(const std::string& GroupName, const Epetra_MultiVector& X, bool writeTranspose)
1005 {
1006 
1007  if (!IsOpen())
1008  throw(Exception(__FILE__, __LINE__, "no file open yet"));
1009 
1010  hid_t group_id, dset_id;
1011  hid_t filespace_id, memspace_id;
1012  herr_t status;
1013 
1014  // need a linear distribution to use hyperslabs
1015  Teuchos::RCP<Epetra_MultiVector> LinearX;
1016 
1017  if (X.Map().LinearMap())
1018  LinearX = Teuchos::rcp(const_cast<Epetra_MultiVector*>(&X), false);
1019  else
1020  {
1021  Epetra_Map LinearMap(X.GlobalLength(), X.Map().IndexBase(), Comm_);
1022  LinearX = Teuchos::rcp(new Epetra_MultiVector(LinearMap, X.NumVectors()));
1023  Epetra_Import Importer(LinearMap, X.Map());
1024  LinearX->Import(X, Importer, Insert);
1025  }
1026 
1027  int NumVectors = X.NumVectors();
1028  int GlobalLength = X.GlobalLength();
1029 
1030  // Whether or not we do writeTranspose or not is
1031  // handled by one of the components of q_dimsf, offset and count.
1032  // They are determined by indexT
1033  int indexT(0);
1034  if (writeTranspose) indexT = 1;
1035 
1036  hsize_t q_dimsf[] = {static_cast<hsize_t>(GlobalLength), static_cast<hsize_t>(GlobalLength)};
1037  q_dimsf[indexT] = NumVectors;
1038 
1039  filespace_id = H5Screate_simple(2, q_dimsf, NULL);
1040 
1041  if (!IsContained(GroupName))
1042  CreateGroup(GroupName);
1043 
1044  group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1045 
1046  // Create the dataset with default properties and close filespace_id.
1047  dset_id = H5Dcreate(group_id, "Values", H5T_NATIVE_DOUBLE, filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1048 
1049  // Create property list for collective dataset write.
1050  plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1051 #ifdef HAVE_MPI
1052  H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1053 #endif
1054 
1055 
1056  // Select hyperslab in the file.
1057  hsize_t offset[] = {static_cast<hsize_t>(LinearX->Map().GID(0)-X.Map().IndexBase()),
1058  static_cast<hsize_t>(LinearX->Map().GID(0)-X.Map().IndexBase())};
1059  hsize_t stride[] = {1, 1};
1060  hsize_t count[] = {static_cast<hsize_t>(LinearX->MyLength()),
1061  static_cast<hsize_t>(LinearX->MyLength())};
1062  hsize_t block[] = {1, 1};
1063 
1064  // write vectors one by one
1065  for (int n(0); n < NumVectors; ++n)
1066  {
1067  // Select hyperslab in the file.
1068  offset[indexT] = n;
1069  count [indexT] = 1;
1070 
1071  filespace_id = H5Dget_space(dset_id);
1072  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1073  count, block);
1074 
1075  // Each process defines dataset in memory and writes it to the hyperslab in the file.
1076  hsize_t dimsm[] = {static_cast<hsize_t>(LinearX->MyLength())};
1077  memspace_id = H5Screate_simple(1, dimsm, NULL);
1078 
1079  // Write hyperslab
1080  status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1081  plist_id_, LinearX->operator[](n));
1082  CHECK_STATUS(status);
1083  }
1084  H5Gclose(group_id);
1085  H5Sclose(memspace_id);
1086  H5Sclose(filespace_id);
1087  H5Dclose(dset_id);
1088  H5Pclose(plist_id_);
1089 
1090  Write(GroupName, "GlobalLength", GlobalLength);
1091  Write(GroupName, "NumVectors", NumVectors);
1092  Write(GroupName, "__type__", "Epetra_MultiVector");
1093 }
1094 
1095 // ==========================================================================
1096 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
1097  Epetra_MultiVector*& X, bool writeTranspose)
1098 {
1099  // first read it with linear distribution
1100  Epetra_MultiVector* LinearX;
1101  Read(GroupName, LinearX, writeTranspose, Map.IndexBase());
1102 
1103  // now build the importer to the actual one
1104  Epetra_Import Importer(Map, LinearX->Map());
1105  X = new Epetra_MultiVector(Map, LinearX->NumVectors());
1106  X->Import(*LinearX, Importer, Insert);
1107 
1108  // delete memory
1109  delete LinearX;
1110 }
1111 
1112 // ==========================================================================
1113 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_MultiVector*& LinearX,
1114  bool readTranspose, const int& indexBase)
1115 {
1116  int GlobalLength, NumVectors;
1117 
1118  ReadMultiVectorProperties(GroupName, GlobalLength, NumVectors);
1119 
1120  hid_t group_id;
1121  hid_t memspace_id;
1122 
1123  // Whether or not we do readTranspose or not is
1124  // handled by one of the components of q_dimsf, offset and count.
1125  // They are determined by indexT
1126  int indexT(0);
1127  if (readTranspose) indexT = 1;
1128 
1129  hsize_t q_dimsf[] = {static_cast<hsize_t>(GlobalLength), static_cast<hsize_t>(GlobalLength)};
1130  q_dimsf[indexT] = NumVectors;
1131 
1132  hid_t filespace_id = H5Screate_simple(2, q_dimsf, NULL);
1133 
1134  if (!IsContained(GroupName))
1135  throw(Exception(__FILE__, __LINE__,
1136  "requested group " + GroupName + " not found"));
1137 
1138  group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1139 
1140  // Create the dataset with default properties and close filespace_id.
1141  hid_t dset_id = H5Dopen(group_id, "Values", H5P_DEFAULT);
1142 
1143  // Create property list for collective dataset write.
1144  plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1145 #ifdef HAVE_MPI
1146  H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1147 #endif
1148  H5Pclose(plist_id_);
1149 
1150  Epetra_Map LinearMap(GlobalLength, indexBase, Comm());
1151  LinearX = new Epetra_MultiVector(LinearMap, NumVectors);
1152 
1153  // Select hyperslab in the file.
1154  hsize_t offset[] = {static_cast<hsize_t>(LinearMap.GID(0) - indexBase), static_cast<hsize_t>(LinearMap.GID(0) - indexBase)};
1155  hsize_t stride[] = {1, 1};
1156 
1157  // If readTranspose is false, we can read the data in one shot.
1158  // It would actually be possible to skip this first part and
1159  if (!readTranspose)
1160  {
1161  // Select hyperslab in the file.
1162  hsize_t count[] = {1, 1};
1163  hsize_t block[] = {static_cast<hsize_t>(LinearX->MyLength()), static_cast<hsize_t>(LinearX->MyLength())};
1164 
1165  offset[indexT] = 0;
1166  count [indexT] = NumVectors;
1167  block [indexT] = 1;
1168 
1169  filespace_id = H5Dget_space(dset_id);
1170  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1171  count, block);
1172 
1173  // Each process defines dataset in memory and writes it to the hyperslab in the file.
1174  hsize_t dimsm[] = {NumVectors * static_cast<hsize_t>(LinearX->MyLength())};
1175  memspace_id = H5Screate_simple(1, dimsm, NULL);
1176 
1177  // Write hyperslab
1178  CHECK_STATUS(H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1179  H5P_DEFAULT, LinearX->Values()));
1180 
1181  } else {
1182  // doing exactly the same as in write
1183 
1184  // Select hyperslab in the file.
1185  hsize_t count[] = {static_cast<hsize_t>(LinearX->MyLength()),
1186  static_cast<hsize_t>(LinearX->MyLength())};
1187  hsize_t block[] = {1, 1};
1188 
1189  // write vectors one by one
1190  for (int n(0); n < NumVectors; ++n)
1191  {
1192  // Select hyperslab in the file.
1193  offset[indexT] = n;
1194  count [indexT] = 1;
1195 
1196  filespace_id = H5Dget_space(dset_id);
1197  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1198  count, block);
1199 
1200  // Each process defines dataset in memory and writes it to the hyperslab in the file.
1201  hsize_t dimsm[] = {static_cast<hsize_t>(LinearX->MyLength())};
1202  memspace_id = H5Screate_simple(1, dimsm, NULL);
1203 
1204  // Read hyperslab
1205  CHECK_STATUS(H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1206  H5P_DEFAULT, LinearX->operator[](n)));
1207 
1208  }
1209  }
1210 
1211  CHECK_STATUS(H5Gclose(group_id));
1212  CHECK_STATUS(H5Sclose(memspace_id));
1213  CHECK_STATUS(H5Sclose(filespace_id));
1214  CHECK_STATUS(H5Dclose(dset_id));
1215 }
1216 
1217 // ==========================================================================
1218 void EpetraExt::HDF5::ReadMultiVectorProperties(const std::string& GroupName,
1219  int& GlobalLength,
1220  int& NumVectors)
1221 {
1222  if (!IsContained(GroupName))
1223  throw(Exception(__FILE__, __LINE__,
1224  "requested group " + GroupName + " not found"));
1225 
1226  std::string Label;
1227  Read(GroupName, "__type__", Label);
1228 
1229  if (Label != "Epetra_MultiVector")
1230  throw(Exception(__FILE__, __LINE__,
1231  "requested group " + GroupName + " is not an Epetra_MultiVector",
1232  "__type__ = " + Label));
1233 
1234  Read(GroupName, "GlobalLength", GlobalLength);
1235  Read(GroupName, "NumVectors", NumVectors);
1236 }
1237 
1238 // ========================= //
1239 // EpetraExt::DistArray<int> //
1240 // ========================= //
1241 
1242 // ==========================================================================
1243 void EpetraExt::HDF5::Write(const std::string& GroupName, const DistArray<int>& x)
1244 {
1245  if (x.Map().LinearMap())
1246  {
1247  Write(GroupName, "Values", x.Map().NumMyElements() * x.RowSize(),
1248  x.Map().NumGlobalElements() * x.RowSize(),
1249  H5T_NATIVE_INT, x.Values());
1250  }
1251  else
1252  {
1253  // need to build a linear map first, the import data, then
1254  // finally write them
1255  const Epetra_BlockMap& OriginalMap = x.Map();
1256  Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
1257  Epetra_Import Importer(LinearMap, OriginalMap);
1258 
1259  EpetraExt::DistArray<int> LinearX(LinearMap, x.RowSize());
1260  LinearX.Import(x, Importer, Insert);
1261 
1262  Write(GroupName, "Values", LinearMap.NumMyElements() * x.RowSize(),
1263  LinearMap.NumGlobalElements() * x.RowSize(),
1264  H5T_NATIVE_INT, LinearX.Values());
1265  }
1266 
1267  Write(GroupName, "__type__", "EpetraExt::DistArray<int>");
1268  Write(GroupName, "GlobalLength", x.GlobalLength());
1269  Write(GroupName, "RowSize", x.RowSize());
1270 }
1271 
1272 // ==========================================================================
1273 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
1274  DistArray<int>*& X)
1275 {
1276  // gets the length of the std::vector
1277  int GlobalLength, RowSize;
1278  ReadIntDistArrayProperties(GroupName, GlobalLength, RowSize);
1279 
1280  if (Map.LinearMap())
1281  {
1282  X = new EpetraExt::DistArray<int>(Map, RowSize);
1283  // simply read stuff and go home
1284  Read(GroupName, "Values", Map.NumMyElements() * RowSize,
1285  Map.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->Values());
1286  }
1287  else
1288  {
1289  // we need to first create a linear map, read the std::vector,
1290  // then import it to the actual nonlinear map
1291  Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
1292  EpetraExt::DistArray<int> LinearX(LinearMap, RowSize);
1293 
1294  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1295  LinearMap.NumGlobalElements() * RowSize,
1296  H5T_NATIVE_INT, LinearX.Values());
1297 
1298  Epetra_Import Importer(Map, LinearMap);
1299  X = new EpetraExt::DistArray<int>(Map, RowSize);
1300  X->Import(LinearX, Importer, Insert);
1301  }
1302 }
1303 
1304 // ==========================================================================
1305 void EpetraExt::HDF5::Read(const std::string& GroupName, DistArray<int>*& X)
1306 {
1307  // gets the length of the std::vector
1308  int GlobalLength, RowSize;
1309  ReadIntDistArrayProperties(GroupName, GlobalLength, RowSize);
1310 
1311  // creates a new linear map and use it to read data
1312  Epetra_Map LinearMap(GlobalLength, 0, Comm_);
1313  X = new EpetraExt::DistArray<int>(LinearMap, RowSize);
1314 
1315  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1316  LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->Values());
1317 }
1318 
1319 // ==========================================================================
1320 void EpetraExt::HDF5::ReadIntDistArrayProperties(const std::string& GroupName,
1321  int& GlobalLength,
1322  int& RowSize)
1323 {
1324  if (!IsContained(GroupName))
1325  throw(Exception(__FILE__, __LINE__,
1326  "requested group " + GroupName + " not found"));
1327 
1328  std::string Label;
1329  Read(GroupName, "__type__", Label);
1330 
1331  if (Label != "EpetraExt::DistArray<int>")
1332  throw(Exception(__FILE__, __LINE__,
1333  "requested group " + GroupName + " is not an EpetraExt::DistArray<int>",
1334  "__type__ = " + Label));
1335 
1336  Read(GroupName, "GlobalLength", GlobalLength);
1337  Read(GroupName, "RowSize", RowSize);
1338 }
1339 
1340 // ============================ //
1341 // EpetraExt::DistArray<double> //
1342 // ============================ //
1343 
1344 // ==========================================================================
1345 void EpetraExt::HDF5::Write(const std::string& GroupName, const DistArray<double>& x)
1346 {
1347  if (x.Map().LinearMap())
1348  {
1349  Write(GroupName, "Values", x.Map().NumMyElements() * x.RowSize(),
1350  x.Map().NumGlobalElements() * x.RowSize(),
1351  H5T_NATIVE_DOUBLE, x.Values());
1352  }
1353  else
1354  {
1355  // need to build a linear map first, the import data, then
1356  // finally write them
1357  const Epetra_BlockMap& OriginalMap = x.Map();
1358  Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
1359  Epetra_Import Importer(LinearMap, OriginalMap);
1360 
1361  EpetraExt::DistArray<double> LinearX(LinearMap, x.RowSize());
1362  LinearX.Import(x, Importer, Insert);
1363 
1364  Write(GroupName, "Values", LinearMap.NumMyElements() * x.RowSize(),
1365  LinearMap.NumGlobalElements() * x.RowSize(),
1366  H5T_NATIVE_DOUBLE, LinearX.Values());
1367  }
1368 
1369  Write(GroupName, "__type__", "EpetraExt::DistArray<double>");
1370  Write(GroupName, "GlobalLength", x.GlobalLength());
1371  Write(GroupName, "RowSize", x.RowSize());
1372 }
1373 
1374 // ==========================================================================
1375 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
1376  DistArray<double>*& X)
1377 {
1378  // gets the length of the std::vector
1379  int GlobalLength, RowSize;
1380  ReadDoubleDistArrayProperties(GroupName, GlobalLength, RowSize);
1381 
1382  if (Map.LinearMap())
1383  {
1384  X = new EpetraExt::DistArray<double>(Map, RowSize);
1385  // simply read stuff and go home
1386  Read(GroupName, "Values", Map.NumMyElements() * RowSize,
1387  Map.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->Values());
1388  }
1389  else
1390  {
1391  // we need to first create a linear map, read the std::vector,
1392  // then import it to the actual nonlinear map
1393  Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
1394  EpetraExt::DistArray<double> LinearX(LinearMap, RowSize);
1395 
1396  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1397  LinearMap.NumGlobalElements() * RowSize,
1398  H5T_NATIVE_DOUBLE, LinearX.Values());
1399 
1400  Epetra_Import Importer(Map, LinearMap);
1401  X = new EpetraExt::DistArray<double>(Map, RowSize);
1402  X->Import(LinearX, Importer, Insert);
1403  }
1404 }
1405 
1406 // ==========================================================================
1407 void EpetraExt::HDF5::Read(const std::string& GroupName, DistArray<double>*& X)
1408 {
1409  // gets the length of the std::vector
1410  int GlobalLength, RowSize;
1411  ReadDoubleDistArrayProperties(GroupName, GlobalLength, RowSize);
1412 
1413  // creates a new linear map and use it to read data
1414  Epetra_Map LinearMap(GlobalLength, 0, Comm_);
1415  X = new EpetraExt::DistArray<double>(LinearMap, RowSize);
1416 
1417  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1418  LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->Values());
1419 }
1420 //
1421 // ==========================================================================
1422 void EpetraExt::HDF5::ReadDoubleDistArrayProperties(const std::string& GroupName,
1423  int& GlobalLength,
1424  int& RowSize)
1425 {
1426  if (!IsContained(GroupName))
1427  throw(Exception(__FILE__, __LINE__,
1428  "requested group " + GroupName + " not found"));
1429 
1430  std::string Label;
1431  Read(GroupName, "__type__", Label);
1432 
1433  if (Label != "EpetraExt::DistArray<double>")
1434  throw(Exception(__FILE__, __LINE__,
1435  "requested group " + GroupName + " is not an EpetraExt::DistArray<double>",
1436  "__type__ = " + Label));
1437 
1438  Read(GroupName, "GlobalLength", GlobalLength);
1439  Read(GroupName, "RowSize", RowSize);
1440 }
1441 
1442 // ======================= //
1443 // basic serial data types //
1444 // ======================= //
1445 
1446 // ==========================================================================
1447 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1448  int what)
1449 {
1450  if (!IsContained(GroupName))
1451  CreateGroup(GroupName);
1452 
1453  hid_t filespace_id = H5Screate(H5S_SCALAR);
1454  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1455  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_INT,
1456  filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1457 
1458  herr_t status = H5Dwrite(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
1459  H5P_DEFAULT, &what);
1460  CHECK_STATUS(status);
1461 
1462  // Close/release resources.
1463  H5Dclose(dset_id);
1464  H5Gclose(group_id);
1465  H5Sclose(filespace_id);
1466 }
1467 
1468 // ==========================================================================
1469 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1470  double what)
1471 {
1472  if (!IsContained(GroupName))
1473  CreateGroup(GroupName);
1474 
1475  hid_t filespace_id = H5Screate(H5S_SCALAR);
1476  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1477  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_DOUBLE,
1478  filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1479 
1480  herr_t status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL,
1481  filespace_id, H5P_DEFAULT, &what);
1482  CHECK_STATUS(status);
1483 
1484  // Close/release resources.
1485  H5Dclose(dset_id);
1486  H5Gclose(group_id);
1487  H5Sclose(filespace_id);
1488 }
1489 
1490 // ==========================================================================
1491 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName, int& data)
1492 {
1493  if (!IsContained(GroupName))
1494  throw(Exception(__FILE__, __LINE__,
1495  "requested group " + GroupName + " not found"));
1496 
1497  // Create group in the root group using absolute name.
1498  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1499 
1500  hid_t filespace_id = H5Screate(H5S_SCALAR);
1501  hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1502 
1503  herr_t status = H5Dread(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
1504  H5P_DEFAULT, &data);
1505  CHECK_STATUS(status);
1506 
1507  H5Sclose(filespace_id);
1508  H5Dclose(dset_id);
1509  H5Gclose(group_id);
1510 }
1511 
1512 // ==========================================================================
1513 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName, double& data)
1514 {
1515  if (!IsContained(GroupName))
1516  throw(Exception(__FILE__, __LINE__,
1517  "requested group " + GroupName + " not found"));
1518 
1519  // Create group in the root group using absolute name.
1520  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1521 
1522  hid_t filespace_id = H5Screate(H5S_SCALAR);
1523  hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1524 
1525  herr_t status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, filespace_id,
1526  H5P_DEFAULT, &data);
1527  CHECK_STATUS(status);
1528 
1529  H5Sclose(filespace_id);
1530  H5Dclose(dset_id);
1531  H5Gclose(group_id);
1532 }
1533 
1534 // ==========================================================================
1535 void EpetraExt::HDF5::Write(const std::string& GroupName,
1536  const std::string& DataSetName,
1537  const std::string& data)
1538 {
1539  if (!IsContained(GroupName))
1540  CreateGroup(GroupName);
1541 
1542  hsize_t len = 1;
1543 
1544  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1545 
1546  hid_t dataspace_id = H5Screate_simple(1, &len, NULL);
1547 
1548  hid_t atype = H5Tcopy(H5T_C_S1);
1549  H5Tset_size(atype, data.size() + 1);
1550 
1551  hid_t dataset_id = H5Dcreate(group_id, DataSetName.c_str(), atype,
1552  dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1553  CHECK_STATUS(H5Dwrite(dataset_id, atype, H5S_ALL, H5S_ALL,
1554  H5P_DEFAULT, data.c_str()));
1555 
1556  CHECK_STATUS(H5Tclose(atype));
1557 
1558  CHECK_STATUS(H5Dclose(dataset_id));
1559 
1560  CHECK_STATUS(H5Sclose(dataspace_id));
1561 
1562  CHECK_STATUS(H5Gclose(group_id));
1563 }
1564 
1565 // ==========================================================================
1566 void EpetraExt::HDF5::Read(const std::string& GroupName,
1567  const std::string& DataSetName,
1568  std::string& data)
1569 {
1570  if (!IsContained(GroupName))
1571  throw(Exception(__FILE__, __LINE__,
1572  "requested group " + GroupName + " not found"));
1573 
1574  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1575 
1576  hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1577 
1578  hid_t datatype_id = H5Dget_type(dataset_id);
1579 // size_t typesize_id = H5Tget_size(datatype_id);
1580  H5T_class_t typeclass_id = H5Tget_class(datatype_id);
1581 
1582  if(typeclass_id != H5T_STRING)
1583  throw(Exception(__FILE__, __LINE__,
1584  "requested group " + GroupName + " is not a std::string"));
1585  char data2[160];
1586  CHECK_STATUS(H5Dread(dataset_id, datatype_id, H5S_ALL, H5S_ALL,
1587  H5P_DEFAULT, data2));
1588  data = data2;
1589 
1590  H5Dclose(dataset_id);
1591  H5Gclose(group_id);
1592 }
1593 
1594 // ============= //
1595 // serial arrays //
1596 // ============= //
1597 
1598 // ==========================================================================
1599 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1600  const int type, const int Length,
1601  void* data)
1602 {
1603  if (!IsContained(GroupName))
1604  CreateGroup(GroupName);
1605 
1606  hsize_t dimsf = Length;
1607 
1608  hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
1609 
1610  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1611 
1612  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type,
1613  filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1614 
1615  herr_t status = H5Dwrite(dset_id, type, H5S_ALL, H5S_ALL,
1616  H5P_DEFAULT, data);
1617  CHECK_STATUS(status);
1618 
1619  // Close/release resources.
1620  H5Dclose(dset_id);
1621  H5Gclose(group_id);
1622  H5Sclose(filespace_id);
1623 }
1624 
1625 // ==========================================================================
1626 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName,
1627  const int type, const int Length,
1628  void* data)
1629 {
1630  if (!IsContained(GroupName))
1631  throw(Exception(__FILE__, __LINE__,
1632  "requested group " + GroupName + " not found"));
1633 
1634  hsize_t dimsf = Length;
1635 
1636  hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
1637 
1638  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1639 
1640  hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1641 
1642  herr_t status = H5Dread(dset_id, type, H5S_ALL, filespace_id,
1643  H5P_DEFAULT, data);
1644  CHECK_STATUS(status);
1645 
1646  // Close/release resources.
1647  H5Dclose(dset_id);
1648  H5Gclose(group_id);
1649  H5Sclose(filespace_id);
1650 }
1651 
1652 // ================== //
1653 // distributed arrays //
1654 // ================== //
1655 
1656 // ==========================================================================
1657 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1658  int MySize, int GlobalSize, int type, const void* data)
1659 {
1660  int Offset;
1661  Comm_.ScanSum(&MySize, &Offset, 1);
1662  Offset -= MySize;
1663 
1664  hsize_t MySize_t = MySize;
1665  hsize_t GlobalSize_t = GlobalSize;
1666  hsize_t Offset_t = Offset;
1667 
1668  hid_t filespace_id = H5Screate_simple(1, &GlobalSize_t, NULL);
1669 
1670  // Create the dataset with default properties and close filespace.
1671  if (!IsContained(GroupName))
1672  CreateGroup(GroupName);
1673 
1674  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1675  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type, filespace_id,
1676  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1677 
1678  H5Sclose(filespace_id);
1679 
1680  // Each process defines dataset in memory and writes it to the hyperslab
1681  // in the file.
1682 
1683  hid_t memspace_id = H5Screate_simple(1, &MySize_t, NULL);
1684 
1685  // Select hyperslab in the file.
1686  filespace_id = H5Dget_space(dset_id);
1687  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL, &MySize_t, NULL);
1688 
1689  plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1690 #ifdef HAVE_MPI
1691  H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1692 #endif
1693 
1694  status = H5Dwrite(dset_id, type, memspace_id, filespace_id,
1695  plist_id_, data);
1696  CHECK_STATUS(status);
1697 
1698  // Close/release resources.
1699  H5Dclose(dset_id);
1700  H5Gclose(group_id);
1701  H5Sclose(filespace_id);
1702  H5Sclose(memspace_id);
1703  H5Pclose(plist_id_);
1704 }
1705 
1706 // ==========================================================================
1707 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName,
1708  int MySize, int GlobalSize,
1709  const int type, void* data)
1710 {
1711  if (!IsOpen())
1712  throw(Exception(__FILE__, __LINE__, "no file open yet"));
1713 
1714  hsize_t MySize_t = MySize;
1715 
1716  // offset
1717  int itmp;
1718  Comm_.ScanSum(&MySize, &itmp, 1);
1719  hsize_t Offset_t = itmp - MySize;
1720 
1721  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1722  hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1723  //hid_t space_id = H5Screate_simple(1, &Offset_t, 0);
1724 
1725  // Select hyperslab in the file.
1726  hid_t filespace_id = H5Dget_space(dataset_id);
1727  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL,
1728  &MySize_t, NULL);
1729 
1730  hid_t mem_dataspace = H5Screate_simple (1, &MySize_t, NULL);
1731 
1732  herr_t status = H5Dread(dataset_id, type, mem_dataspace, filespace_id,
1733  H5P_DEFAULT, data);
1734  CHECK_STATUS(status);
1735 
1736  H5Sclose(mem_dataspace);
1737  H5Gclose(group_id);
1738  //H5Sclose(space_id);
1739  H5Dclose(dataset_id);
1740 // H5Dclose(filespace_id);
1741 }
1742 
1743 
1744 #endif
1745 
void ReadDoubleDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
void Write(const std::string &GroupName, const std::string &DataSetName, int data)
Write an integer in group GroupName using the given DataSetName.
#define CHECK_HID(hid_t)
void ReadMapProperties(const std::string &GroupName, int &NumGlobalElements, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_Map.
static void WriteParameterListRecursive(const Teuchos::ParameterList &params, hid_t group_id)
void Read(const std::string &GroupName, const std::string &DataSetName, int &data)
Read an integer from group /GroupName/DataSetName.
T * Values()
Returns a pointer to the internally stored data (non-const version).
int RowSize() const
Returns the row size, that is, the amount of data associated with each element.
void ReadBlockMapProperties(const std::string &GroupName, int &NumGlobalElements, int &NumGlobalPoints, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_BlockMap.
void ReadIntVectorProperties(const std::string &GroupName, int &GlobalLength)
Read basic properties of specified Epetra_IntVector.
void ReadCrsGraphProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumGlobalNonzeros, int &NumGlobalDiagonals, int &MaxNumIndices)
Read basic properties of specified Epetra_CrsGraph.
bool IsContained(const std::string Name)
Return true if Name is contained in the database.
int GlobalLength() const
Returns the global length of the array.
void Open(const std::string FileName, int AccessType=H5F_ACC_RDWR)
Open specified file with given access type.
#define CHECK_STATUS(status)
static herr_t f_operator(hid_t loc_id, const char *name, void *opdata)
HDF5(const Epetra_Comm &Comm)
Constructor.
void Create(const std::string FileName)
Create a new file.
void ReadIntDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
void ReadCrsMatrixProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumNonzeros, int &NumGlobalDiagonals, int &MaxNumEntries, double &NormOne, double &NormInf)
Read basic properties of specified Epetra_CrsMatrix.
static herr_t FindDataset(hid_t loc_id, const char *name, void *opdata)
std::string name
void ReadMultiVectorProperties(const std::string &GroupName, int &GlobalLength, int &NumVectors)
Read basic properties of specified Epetra_MultiVector.
std::string toString(const int &x)
DistArray<T>: A class to store row-oriented multivectors of type T.