EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
90 struct FindDataset_t
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, std::string GroupName)
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  // recursively look for groups
405  size_t pos = Name.find("/");
406  if (pos != std::string::npos)
407  {
408  std::string NewGroupName = Name.substr(0, pos);
409  if (GroupName != "")
410  NewGroupName = GroupName + "/" + NewGroupName;
411  std::string NewName = Name.substr(pos + 1);
412  return IsContained(NewName, NewGroupName);
413  }
414 
415  GroupName = "/" + GroupName;
416 
417  //int idx_f =
418  H5Giterate(file_id_, GroupName.c_str(), NULL, FindDataset, (void*)&data);
419 
420  return(data.found);
421 }
422 
423 // ============================ //
424 // Epetra_BlockMap / Epetra_Map //
425 // ============================ //
426 
427 // ==========================================================================
428 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_BlockMap& BlockMap)
429 {
430  int NumMyPoints = BlockMap.NumMyPoints();
431  int NumMyElements = BlockMap.NumMyElements();
432  int NumGlobalElements = BlockMap.NumGlobalElements();
433  int NumGlobalPoints = BlockMap.NumGlobalPoints();
434  int* MyGlobalElements = BlockMap.MyGlobalElements();
435  int* ElementSizeList = BlockMap.ElementSizeList();
436 
437  std::vector<int> NumMyElements_v(Comm_.NumProc());
438  Comm_.GatherAll(&NumMyElements, &NumMyElements_v[0], 1);
439 
440  std::vector<int> NumMyPoints_v(Comm_.NumProc());
441  Comm_.GatherAll(&NumMyPoints, &NumMyPoints_v[0], 1);
442 
443  Write(Name, "MyGlobalElements", NumMyElements, NumGlobalElements,
444  H5T_NATIVE_INT, MyGlobalElements);
445  Write(Name, "ElementSizeList", NumMyElements, NumGlobalElements,
446  H5T_NATIVE_INT, ElementSizeList);
447  Write(Name, "NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
448 
449  // need to know how many processors currently host this map
450  Write(Name, "NumProc", Comm_.NumProc());
451  // write few more data about this map
452  Write(Name, "NumGlobalPoints", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalPoints);
453  Write(Name, "NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalElements);
454  Write(Name, "IndexBase", BlockMap.IndexBase());
455  Write(Name, "__type__", "Epetra_BlockMap");
456 }
457 
458 // ==========================================================================
459 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_BlockMap*& BlockMap)
460 {
461  int NumGlobalElements, NumGlobalPoints, IndexBase, NumProc;
462 
463  ReadBlockMapProperties(GroupName, NumGlobalElements, NumGlobalPoints,
464  IndexBase, NumProc);
465 
466  std::vector<int> NumMyPoints_v(Comm_.NumProc());
467  std::vector<int> NumMyElements_v(Comm_.NumProc());
468 
469  Read(GroupName, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
470  Read(GroupName, "NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
471  int NumMyElements = NumMyElements_v[Comm_.MyPID()];
472 // int NumMyPoints = NumMyPoints_v[Comm_.MyPID()];
473 
474  if (NumProc != Comm_.NumProc())
475  throw(Exception(__FILE__, __LINE__,
476  "requested map not compatible with current number of",
477  "processors, " + toString(Comm_.NumProc()) +
478  " vs. " + toString(NumProc)));
479 
480  std::vector<int> MyGlobalElements(NumMyElements);
481  std::vector<int> ElementSizeList(NumMyElements);
482 
483  Read(GroupName, "MyGlobalElements", NumMyElements, NumGlobalElements,
484  H5T_NATIVE_INT, &MyGlobalElements[0]);
485 
486  Read(GroupName, "ElementSizeList", NumMyElements, NumGlobalElements,
487  H5T_NATIVE_INT, &ElementSizeList[0]);
488 
489  BlockMap = new Epetra_BlockMap(NumGlobalElements, NumMyElements,
490  &MyGlobalElements[0], &ElementSizeList[0],
491  IndexBase, Comm_);
492 }
493 
494 // ==========================================================================
495 void EpetraExt::HDF5::ReadBlockMapProperties(const std::string& GroupName,
496  int& NumGlobalElements,
497  int& NumGlobalPoints,
498  int& IndexBase,
499  int& NumProc)
500 {
501  if (!IsContained(GroupName))
502  throw(Exception(__FILE__, __LINE__,
503  "requested group " + GroupName + " not found"));
504 
505  std::string Label;
506  Read(GroupName, "__type__", Label);
507 
508  if (Label != "Epetra_BlockMap")
509  throw(Exception(__FILE__, __LINE__,
510  "requested group " + GroupName + " is not an Epetra_BlockMap",
511  "__type__ = " + Label));
512 
513  Read(GroupName, "NumGlobalElements", NumGlobalElements);
514  Read(GroupName, "NumGlobalPoints", NumGlobalPoints);
515  Read(GroupName, "IndexBase", IndexBase);
516  Read(GroupName, "NumProc", NumProc);
517 }
518 
519 // ==========================================================================
520 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_Map& Map)
521 {
522  int MySize = Map.NumMyElements();
523  int GlobalSize = Map.NumGlobalElements();
524  int* MyGlobalElements = Map.MyGlobalElements();
525 
526  std::vector<int> NumMyElements(Comm_.NumProc());
527  Comm_.GatherAll(&MySize, &NumMyElements[0], 1);
528 
529  Write(Name, "MyGlobalElements", MySize, GlobalSize,
530  H5T_NATIVE_INT, MyGlobalElements);
531  Write(Name, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements[0]);
532  Write(Name, "NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &GlobalSize);
533  Write(Name, "NumProc", Comm_.NumProc());
534  Write(Name, "IndexBase", Map.IndexBase());
535  Write(Name, "__type__", "Epetra_Map");
536 }
537 
538 // ==========================================================================
539 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_Map*& Map)
540 {
541  int NumGlobalElements, IndexBase, NumProc;
542 
543  ReadMapProperties(GroupName, NumGlobalElements, IndexBase, NumProc);
544 
545  std::vector<int> NumMyElements_v(Comm_.NumProc());
546 
547  Read(GroupName, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
548  int NumMyElements = NumMyElements_v[Comm_.MyPID()];
549 
550  if (NumProc != Comm_.NumProc())
551  throw(Exception(__FILE__, __LINE__,
552  "requested map not compatible with current number of",
553  "processors, " + toString(Comm_.NumProc()) +
554  " vs. " + toString(NumProc)));
555 
556  std::vector<int> MyGlobalElements(NumMyElements);
557 
558  Read(GroupName, "MyGlobalElements", NumMyElements, NumGlobalElements,
559  H5T_NATIVE_INT, &MyGlobalElements[0]);
560 
561  Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], IndexBase, Comm_);
562 }
563 
564 // ==========================================================================
565 void EpetraExt::HDF5::ReadMapProperties(const std::string& GroupName,
566  int& NumGlobalElements,
567  int& IndexBase,
568  int& NumProc)
569 {
570  if (!IsContained(GroupName))
571  throw(Exception(__FILE__, __LINE__,
572  "requested group " + GroupName + " not found"));
573 
574  std::string Label;
575  Read(GroupName, "__type__", Label);
576 
577  if (Label != "Epetra_Map")
578  throw(Exception(__FILE__, __LINE__,
579  "requested group " + GroupName + " is not an Epetra_Map",
580  "__type__ = " + Label));
581 
582  Read(GroupName, "NumGlobalElements", NumGlobalElements);
583  Read(GroupName, "IndexBase", IndexBase);
584  Read(GroupName, "NumProc", NumProc);
585 }
586 
587 // ================ //
588 // Epetra_IntVector //
589 // ================ //
590 
591 // ==========================================================================
592 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_IntVector& x)
593 {
594  if (x.Map().LinearMap())
595  {
596  Write(Name, "GlobalLength", x.GlobalLength());
597  Write(Name, "Values", x.Map().NumMyElements(), x.Map().NumGlobalElements(),
598  H5T_NATIVE_INT, x.Values());
599  }
600  else
601  {
602  // need to build a linear map first, the import data, then
603  // finally write them
604  const Epetra_BlockMap& OriginalMap = x.Map();
605  Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
606  Epetra_Import Importer(LinearMap, OriginalMap);
607 
608  Epetra_IntVector LinearX(LinearMap);
609  LinearX.Import(x, Importer, Insert);
610 
611  Write(Name, "GlobalLength", x.GlobalLength());
612  Write(Name, "Values", LinearMap.NumMyElements(), LinearMap.NumGlobalElements(),
613  H5T_NATIVE_INT, LinearX.Values());
614  }
615  Write(Name, "__type__", "Epetra_IntVector");
616 }
617 
618 // ==========================================================================
619 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_IntVector*& X)
620 {
621  // gets the length of the std::vector
622  int GlobalLength;
623  ReadIntVectorProperties(GroupName, GlobalLength);
624 
625  // creates a new linear map and use it to read data
626  Epetra_Map LinearMap(GlobalLength, 0, Comm_);
627  X = new Epetra_IntVector(LinearMap);
628 
629  Read(GroupName, "Values", LinearMap.NumMyElements(),
630  LinearMap.NumGlobalElements(), H5T_NATIVE_INT, X->Values());
631 }
632 
633 // ==========================================================================
634 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
635  Epetra_IntVector*& X)
636 {
637  // gets the length of the std::vector
638  int GlobalLength;
639  ReadIntVectorProperties(GroupName, GlobalLength);
640 
641  if (Map.LinearMap())
642  {
643  X = new Epetra_IntVector(Map);
644  // simply read stuff and go home
645  Read(GroupName, "Values", Map.NumMyElements(), Map.NumGlobalElements(),
646  H5T_NATIVE_INT, X->Values());
647 
648  }
649  else
650  {
651  // we need to first create a linear map, read the std::vector,
652  // then import it to the actual nonlinear map
653  Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
654  Epetra_IntVector LinearX(LinearMap);
655 
656  Read(GroupName, "Values", LinearMap.NumMyElements(),
657  LinearMap.NumGlobalElements(),
658  H5T_NATIVE_INT, LinearX.Values());
659 
660  Epetra_Import Importer(Map, LinearMap);
661  X = new Epetra_IntVector(Map);
662  X->Import(LinearX, Importer, Insert);
663  }
664 }
665 
666 // ==========================================================================
667 void EpetraExt::HDF5::ReadIntVectorProperties(const std::string& GroupName,
668  int& GlobalLength)
669 {
670  if (!IsContained(GroupName))
671  throw(Exception(__FILE__, __LINE__,
672  "requested group " + GroupName + " not found"));
673 
674  std::string Label;
675  Read(GroupName, "__type__", Label);
676 
677  if (Label != "Epetra_IntVector")
678  throw(Exception(__FILE__, __LINE__,
679  "requested group " + GroupName + " is not an Epetra_IntVector",
680  "__type__ = " + Label));
681 
682  Read(GroupName, "GlobalLength", GlobalLength);
683 }
684 
685 // =============== //
686 // Epetra_CrsGraph //
687 // =============== //
688 
689 // ==========================================================================
690 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_CrsGraph& Graph)
691 {
692  if (!Graph.Filled())
693  throw(Exception(__FILE__, __LINE__,
694  "input Epetra_CrsGraph is not FillComplete()'d"));
695 
696  // like RowMatrix, only without values
697  int MySize = Graph.NumMyNonzeros();
698  int GlobalSize = Graph.NumGlobalNonzeros();
699  std::vector<int> ROW(MySize);
700  std::vector<int> COL(MySize);
701 
702  int count = 0;
703  int* RowIndices;
704  int NumEntries;
705 
706  for (int i = 0; i < Graph.NumMyRows(); ++i)
707  {
708  Graph.ExtractMyRowView(i, NumEntries, RowIndices);
709  for (int j = 0; j < NumEntries; ++j)
710  {
711  ROW[count] = Graph.GRID(i);
712  COL[count] = Graph.GCID(RowIndices[j]);
713  ++count;
714  }
715  }
716 
717  Write(Name, "ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
718  Write(Name, "COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
719  Write(Name, "MaxNumIndices", Graph.MaxNumIndices());
720  Write(Name, "NumGlobalRows", Graph.NumGlobalRows());
721  Write(Name, "NumGlobalCols", Graph.NumGlobalCols());
722  Write(Name, "NumGlobalNonzeros", Graph.NumGlobalNonzeros());
723  Write(Name, "NumGlobalDiagonals", Graph.NumGlobalDiagonals());
724  Write(Name, "__type__", "Epetra_CrsGraph");
725 }
726 
727 // ==========================================================================
728 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_CrsGraph*& Graph)
729 {
730  int NumGlobalRows, NumGlobalCols;
731  int NumGlobalNonzeros, NumGlobalDiagonals, MaxNumIndices;
732 
733  ReadCrsGraphProperties(GroupName, NumGlobalRows,
734  NumGlobalCols, NumGlobalNonzeros,
735  NumGlobalDiagonals, MaxNumIndices);
736 
737  Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
738  Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
739 
740  Read(GroupName, DomainMap, RangeMap, Graph);
741 }
742 
743 // ==========================================================================
744 void EpetraExt::HDF5::ReadCrsGraphProperties(const std::string& GroupName,
745  int& NumGlobalRows,
746  int& NumGlobalCols,
747  int& NumGlobalNonzeros,
748  int& NumGlobalDiagonals,
749  int& MaxNumIndices)
750 {
751  if (!IsContained(GroupName))
752  throw(Exception(__FILE__, __LINE__,
753  "requested group " + GroupName + " not found"));
754 
755  std::string Label;
756  Read(GroupName, "__type__", Label);
757 
758  if (Label != "Epetra_CrsGraph")
759  throw(Exception(__FILE__, __LINE__,
760  "requested group " + GroupName + " is not an Epetra_CrsGraph",
761  "__type__ = " + Label));
762 
763  Read(GroupName, "NumGlobalRows", NumGlobalRows);
764  Read(GroupName, "NumGlobalCols", NumGlobalCols);
765  Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
766  Read(GroupName, "NumGlobalDiagonals", NumGlobalDiagonals);
767  Read(GroupName, "MaxNumIndices", MaxNumIndices);
768 }
769 
770 // ==========================================================================
771 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& DomainMap,
772  const Epetra_Map& RangeMap, Epetra_CrsGraph*& Graph)
773 {
774  if (!IsContained(GroupName))
775  throw(Exception(__FILE__, __LINE__,
776  "requested group " + GroupName + " not found in database"));
777 
778  std::string Label;
779  Read(GroupName, "__type__", Label);
780 
781  if (Label != "Epetra_CrsGraph")
782  throw(Exception(__FILE__, __LINE__,
783  "requested group " + GroupName + " is not an Epetra_CrsGraph",
784  "__type__ = " + Label));
785 
786  int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
787  Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
788  Read(GroupName, "NumGlobalRows", NumGlobalRows);
789  Read(GroupName, "NumGlobalCols", NumGlobalCols);
790 
791  // linear distribution for nonzeros
792  int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
793  if (Comm_.MyPID() == 0)
794  NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
795 
796  std::vector<int> ROW(NumMyNonzeros);
797  Read(GroupName, "ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
798 
799  std::vector<int> COL(NumMyNonzeros);
800  Read(GroupName, "COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
801 
802  Epetra_FECrsGraph* Graph2 = new Epetra_FECrsGraph(Copy, DomainMap, 0);
803 
804  for (int i = 0; i < NumMyNonzeros; )
805  {
806  int count = 1;
807  while (ROW[i + count] == ROW[i])
808  ++count;
809 
810  Graph2->InsertGlobalIndices(1, &ROW[i], count, &COL[i]);
811  i += count;
812  }
813 
814  Graph2->FillComplete(DomainMap, RangeMap);
815 
816  Graph = Graph2;
817 }
818 
819 // ================ //
820 // Epetra_RowMatrix //
821 // ================ //
822 
823 // ==========================================================================
824 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_RowMatrix& Matrix)
825 {
826  if (!Matrix.Filled())
827  throw(Exception(__FILE__, __LINE__,
828  "input Epetra_RowMatrix is not FillComplete()'d"));
829 
830  int MySize = Matrix.NumMyNonzeros();
831  int GlobalSize = Matrix.NumGlobalNonzeros();
832  std::vector<int> ROW(MySize);
833  std::vector<int> COL(MySize);
834  std::vector<double> VAL(MySize);
835 
836  int count = 0;
837  int Length = Matrix.MaxNumEntries();
838  std::vector<int> Indices(Length);
839  std::vector<double> Values(Length);
840  int NumEntries;
841 
842  for (int i = 0; i < Matrix.NumMyRows(); ++i)
843  {
844  Matrix.ExtractMyRowCopy(i, Length, NumEntries, &Values[0], &Indices[0]);
845  for (int j = 0; j < NumEntries; ++j)
846  {
847  ROW[count] = Matrix.RowMatrixRowMap().GID(i);
848  COL[count] = Matrix.RowMatrixColMap().GID(Indices[j]);
849  VAL[count] = Values[j];
850  ++count;
851  }
852  }
853 
854  Write(Name, "ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
855  Write(Name, "COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
856  Write(Name, "VAL", MySize, GlobalSize, H5T_NATIVE_DOUBLE, &VAL[0]);
857  Write(Name, "NumGlobalRows", Matrix.NumGlobalRows());
858  Write(Name, "NumGlobalCols", Matrix.NumGlobalCols());
859  Write(Name, "NumGlobalNonzeros", Matrix.NumGlobalNonzeros());
860  Write(Name, "NumGlobalDiagonals", Matrix.NumGlobalDiagonals());
861  Write(Name, "MaxNumEntries", Matrix.MaxNumEntries());
862  Write(Name, "NormOne", Matrix.NormOne());
863  Write(Name, "NormInf", Matrix.NormInf());
864  Write(Name, "__type__", "Epetra_RowMatrix");
865 }
866 
867 // ==========================================================================
868 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_CrsMatrix*& A)
869 {
870  int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
871  int NumGlobalDiagonals, MaxNumEntries;
872  double NormOne, NormInf;
873 
874  ReadCrsMatrixProperties(GroupName, NumGlobalRows, NumGlobalCols,
875  NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
876  NormOne, NormInf);
877 
878  // build simple linear maps for domain and range space
879  Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
880  Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
881 
882  Read(GroupName, DomainMap, RangeMap, A);
883 }
884 
885 // ==========================================================================
886 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& DomainMap,
887  const Epetra_Map& RangeMap, Epetra_CrsMatrix*& A)
888 {
889  int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
890  int NumGlobalDiagonals, MaxNumEntries;
891  double NormOne, NormInf;
892 
893  ReadCrsMatrixProperties(GroupName, NumGlobalRows, NumGlobalCols,
894  NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
895  NormOne, NormInf);
896 
897  int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
898  if (Comm_.MyPID() == 0)
899  NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
900 
901  std::vector<int> ROW(NumMyNonzeros);
902  Read(GroupName, "ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
903 
904  std::vector<int> COL(NumMyNonzeros);
905  Read(GroupName, "COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
906 
907  std::vector<double> VAL(NumMyNonzeros);
908  Read(GroupName, "VAL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_DOUBLE, &VAL[0]);
909 
910  Epetra_FECrsMatrix* A2 = new Epetra_FECrsMatrix(Copy, DomainMap, 0);
911 
912  for (int i = 0; i < NumMyNonzeros; )
913  {
914  int count = 1;
915  while (ROW[i + count] == ROW[i])
916  ++count;
917 
918  A2->InsertGlobalValues(1, &ROW[i], count, &COL[i], &VAL[i]);
919  i += count;
920  }
921 
922  A2->GlobalAssemble(DomainMap, RangeMap);
923 
924  A = A2;
925 }
926 
927 // ==========================================================================
928 void EpetraExt::HDF5::ReadCrsMatrixProperties(const std::string& GroupName,
929  int& NumGlobalRows,
930  int& NumGlobalCols,
931  int& NumGlobalNonzeros,
932  int& NumGlobalDiagonals,
933  int& MaxNumEntries,
934  double& NormOne,
935  double& NormInf)
936 {
937  if (!IsContained(GroupName))
938  throw(Exception(__FILE__, __LINE__,
939  "requested group " + GroupName + " not found"));
940 
941  std::string Label;
942  Read(GroupName, "__type__", Label);
943 
944  if (Label != "Epetra_RowMatrix")
945  throw(Exception(__FILE__, __LINE__,
946  "requested group " + GroupName + " is not an Epetra_RowMatrix",
947  "__type__ = " + Label));
948 
949  Read(GroupName, "NumGlobalRows", NumGlobalRows);
950  Read(GroupName, "NumGlobalCols", NumGlobalCols);
951  Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
952  Read(GroupName, "NumGlobalDiagonals", NumGlobalDiagonals);
953  Read(GroupName, "MaxNumEntries", MaxNumEntries);
954  Read(GroupName, "NormOne", NormOne);
955  Read(GroupName, "NormInf", NormInf);
956 }
957 
958 // ============= //
959 // ParameterList //
960 // ============= //
961 
962 // ==========================================================================
963 void EpetraExt::HDF5::Write(const std::string& GroupName, const Teuchos::ParameterList& params)
964 {
965  if (!IsOpen())
966  throw(Exception(__FILE__, __LINE__, "no file open yet"));
967 
968  hid_t group_id = H5Gcreate(file_id_, GroupName.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
969 
970  // Iterate through parameter list
971  WriteParameterListRecursive(params, group_id);
972 
973  // Finalize hdf5 file
974  status = H5Gclose(group_id);
975  CHECK_STATUS(status);
976 
977  Write(GroupName, "__type__", "Teuchos::ParameterList");
978 }
979 
980 // ==========================================================================
981 void EpetraExt::HDF5::Read(const std::string& GroupName, Teuchos::ParameterList& params)
982 {
983  if (!IsOpen())
984  throw(Exception(__FILE__, __LINE__, "no file open yet"));
985 
986  std::string Label;
987  Read(GroupName, "__type__", Label);
988 
989  if (Label != "Teuchos::ParameterList")
990  throw(Exception(__FILE__, __LINE__,
991  "requested group " + GroupName + " is not a Teuchos::ParameterList",
992  "__type__ = " + Label));
993 
994  // Open hdf5 file
995  hid_t group_id; /* identifiers */
996  herr_t status;
997 
998  // open group in the root group using absolute name.
999  group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1000  CHECK_HID(group_id);
1001 
1002  // Iterate through parameter list
1003  std::string xxx = "/" + GroupName;
1004  status = H5Giterate(group_id, xxx.c_str() , NULL, f_operator, &params);
1005  CHECK_STATUS(status);
1006 
1007  // Finalize hdf5 file
1008  status = H5Gclose(group_id);
1009  CHECK_STATUS(status);
1010 }
1011 
1012 // ================== //
1013 // Epetra_MultiVector //
1014 // ================== //
1015 
1016 // ==========================================================================
1017 void EpetraExt::HDF5::Write(const std::string& GroupName, const Epetra_MultiVector& X, bool writeTranspose)
1018 {
1019 
1020  if (!IsOpen())
1021  throw(Exception(__FILE__, __LINE__, "no file open yet"));
1022 
1023  hid_t group_id, dset_id;
1024  hid_t filespace_id, memspace_id;
1025  herr_t status;
1026 
1027  // need a linear distribution to use hyperslabs
1028  Teuchos::RCP<Epetra_MultiVector> LinearX;
1029 
1030  if (X.Map().LinearMap())
1031  LinearX = Teuchos::rcp(const_cast<Epetra_MultiVector*>(&X), false);
1032  else
1033  {
1034  Epetra_Map LinearMap(X.GlobalLength(), X.Map().IndexBase(), Comm_);
1035  LinearX = Teuchos::rcp(new Epetra_MultiVector(LinearMap, X.NumVectors()));
1036  Epetra_Import Importer(LinearMap, X.Map());
1037  LinearX->Import(X, Importer, Insert);
1038  }
1039 
1040  int NumVectors = X.NumVectors();
1041  int GlobalLength = X.GlobalLength();
1042 
1043  // Whether or not we do writeTranspose or not is
1044  // handled by one of the components of q_dimsf, offset and count.
1045  // They are determined by indexT
1046  int indexT(0);
1047  if (writeTranspose) indexT = 1;
1048 
1049  hsize_t q_dimsf[] = {static_cast<hsize_t>(GlobalLength), static_cast<hsize_t>(GlobalLength)};
1050  q_dimsf[indexT] = NumVectors;
1051 
1052  filespace_id = H5Screate_simple(2, q_dimsf, NULL);
1053 
1054  if (!IsContained(GroupName))
1055  CreateGroup(GroupName);
1056 
1057  group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1058 
1059  // Create the dataset with default properties and close filespace_id.
1060  dset_id = H5Dcreate(group_id, "Values", H5T_NATIVE_DOUBLE, filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1061 
1062  // Create property list for collective dataset write.
1063  plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1064 #ifdef HAVE_MPI
1065  H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1066 #endif
1067 
1068 
1069  // Select hyperslab in the file.
1070  hsize_t offset[] = {static_cast<hsize_t>(LinearX->Map().GID(0)-X.Map().IndexBase()),
1071  static_cast<hsize_t>(LinearX->Map().GID(0)-X.Map().IndexBase())};
1072  hsize_t stride[] = {1, 1};
1073  hsize_t count[] = {static_cast<hsize_t>(LinearX->MyLength()),
1074  static_cast<hsize_t>(LinearX->MyLength())};
1075  hsize_t block[] = {1, 1};
1076 
1077  // write vectors one by one
1078  for (int n(0); n < NumVectors; ++n)
1079  {
1080  // Select hyperslab in the file.
1081  offset[indexT] = n;
1082  count [indexT] = 1;
1083 
1084  filespace_id = H5Dget_space(dset_id);
1085  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1086  count, block);
1087 
1088  // Each process defines dataset in memory and writes it to the hyperslab in the file.
1089  hsize_t dimsm[] = {static_cast<hsize_t>(LinearX->MyLength())};
1090  memspace_id = H5Screate_simple(1, dimsm, NULL);
1091 
1092  // Write hyperslab
1093  status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1094  plist_id_, LinearX->operator[](n));
1095  CHECK_STATUS(status);
1096  H5Sclose(memspace_id);
1097  }
1098  H5Gclose(group_id);
1099  H5Sclose(filespace_id);
1100  H5Dclose(dset_id);
1101  H5Pclose(plist_id_);
1102 
1103  Write(GroupName, "GlobalLength", GlobalLength);
1104  Write(GroupName, "NumVectors", NumVectors);
1105  Write(GroupName, "__type__", "Epetra_MultiVector");
1106 }
1107 
1108 // ==========================================================================
1109 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
1110  Epetra_MultiVector*& X, bool writeTranspose)
1111 {
1112  // first read it with linear distribution
1113  Epetra_MultiVector* LinearX;
1114  Read(GroupName, LinearX, writeTranspose, Map.IndexBase());
1115 
1116  // now build the importer to the actual one
1117  Epetra_Import Importer(Map, LinearX->Map());
1118  X = new Epetra_MultiVector(Map, LinearX->NumVectors());
1119  X->Import(*LinearX, Importer, Insert);
1120 
1121  // delete memory
1122  delete LinearX;
1123 }
1124 
1125 // ==========================================================================
1126 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_MultiVector*& LinearX,
1127  bool readTranspose, const int& indexBase)
1128 {
1129  int GlobalLength, NumVectors;
1130 
1131  ReadMultiVectorProperties(GroupName, GlobalLength, NumVectors);
1132 
1133  hid_t group_id;
1134  hid_t memspace_id;
1135 
1136  // Whether or not we do readTranspose or not is
1137  // handled by one of the components of q_dimsf, offset and count.
1138  // They are determined by indexT
1139  int indexT(0);
1140  if (readTranspose) indexT = 1;
1141 
1142  hsize_t q_dimsf[] = {static_cast<hsize_t>(GlobalLength), static_cast<hsize_t>(GlobalLength)};
1143  q_dimsf[indexT] = NumVectors;
1144 
1145  hid_t filespace_id = H5Screate_simple(2, q_dimsf, NULL);
1146 
1147  if (!IsContained(GroupName))
1148  throw(Exception(__FILE__, __LINE__,
1149  "requested group " + GroupName + " not found"));
1150 
1151  group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1152 
1153  // Create the dataset with default properties and close filespace_id.
1154  hid_t dset_id = H5Dopen(group_id, "Values", H5P_DEFAULT);
1155 
1156  // Create property list for collective dataset write.
1157  plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1158 #ifdef HAVE_MPI
1159  H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1160 #endif
1161  H5Pclose(plist_id_);
1162 
1163  Epetra_Map LinearMap(GlobalLength, indexBase, Comm());
1164  LinearX = new Epetra_MultiVector(LinearMap, NumVectors);
1165 
1166  // Select hyperslab in the file.
1167  hsize_t offset[] = {static_cast<hsize_t>(LinearMap.GID(0) - indexBase), static_cast<hsize_t>(LinearMap.GID(0) - indexBase)};
1168  hsize_t stride[] = {1, 1};
1169 
1170  // If readTranspose is false, we can read the data in one shot.
1171  // It would actually be possible to skip this first part and
1172  if (!readTranspose)
1173  {
1174  // Select hyperslab in the file.
1175  hsize_t count[] = {1, 1};
1176  hsize_t block[] = {static_cast<hsize_t>(LinearX->MyLength()), static_cast<hsize_t>(LinearX->MyLength())};
1177 
1178  offset[indexT] = 0;
1179  count [indexT] = NumVectors;
1180  block [indexT] = 1;
1181 
1182  filespace_id = H5Dget_space(dset_id);
1183  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1184  count, block);
1185 
1186  // Each process defines dataset in memory and writes it to the hyperslab in the file.
1187  hsize_t dimsm[] = {NumVectors * static_cast<hsize_t>(LinearX->MyLength())};
1188  memspace_id = H5Screate_simple(1, dimsm, NULL);
1189 
1190  // Write hyperslab
1191  CHECK_STATUS(H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1192  H5P_DEFAULT, LinearX->Values()));
1193 
1194  } else {
1195  // doing exactly the same as in write
1196 
1197  // Select hyperslab in the file.
1198  hsize_t count[] = {static_cast<hsize_t>(LinearX->MyLength()),
1199  static_cast<hsize_t>(LinearX->MyLength())};
1200  hsize_t block[] = {1, 1};
1201 
1202  // write vectors one by one
1203  for (int n(0); n < NumVectors; ++n)
1204  {
1205  // Select hyperslab in the file.
1206  offset[indexT] = n;
1207  count [indexT] = 1;
1208 
1209  filespace_id = H5Dget_space(dset_id);
1210  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1211  count, block);
1212 
1213  // Each process defines dataset in memory and writes it to the hyperslab in the file.
1214  hsize_t dimsm[] = {static_cast<hsize_t>(LinearX->MyLength())};
1215  memspace_id = H5Screate_simple(1, dimsm, NULL);
1216 
1217  // Read hyperslab
1218  CHECK_STATUS(H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1219  H5P_DEFAULT, LinearX->operator[](n)));
1220 
1221  }
1222  }
1223 
1224  CHECK_STATUS(H5Gclose(group_id));
1225  CHECK_STATUS(H5Sclose(memspace_id));
1226  CHECK_STATUS(H5Sclose(filespace_id));
1227  CHECK_STATUS(H5Dclose(dset_id));
1228 }
1229 
1230 // ==========================================================================
1231 void EpetraExt::HDF5::ReadMultiVectorProperties(const std::string& GroupName,
1232  int& GlobalLength,
1233  int& NumVectors)
1234 {
1235  if (!IsContained(GroupName))
1236  throw(Exception(__FILE__, __LINE__,
1237  "requested group " + GroupName + " not found"));
1238 
1239  std::string Label;
1240  Read(GroupName, "__type__", Label);
1241 
1242  if (Label != "Epetra_MultiVector")
1243  throw(Exception(__FILE__, __LINE__,
1244  "requested group " + GroupName + " is not an Epetra_MultiVector",
1245  "__type__ = " + Label));
1246 
1247  Read(GroupName, "GlobalLength", GlobalLength);
1248  Read(GroupName, "NumVectors", NumVectors);
1249 }
1250 
1251 // ========================= //
1252 // EpetraExt::DistArray<int> //
1253 // ========================= //
1254 
1255 // ==========================================================================
1256 void EpetraExt::HDF5::Write(const std::string& GroupName, const DistArray<int>& x)
1257 {
1258  if (x.Map().LinearMap())
1259  {
1260  Write(GroupName, "Values", x.Map().NumMyElements() * x.RowSize(),
1261  x.Map().NumGlobalElements() * x.RowSize(),
1262  H5T_NATIVE_INT, x.Values());
1263  }
1264  else
1265  {
1266  // need to build a linear map first, the import data, then
1267  // finally write them
1268  const Epetra_BlockMap& OriginalMap = x.Map();
1269  Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
1270  Epetra_Import Importer(LinearMap, OriginalMap);
1271 
1272  EpetraExt::DistArray<int> LinearX(LinearMap, x.RowSize());
1273  LinearX.Import(x, Importer, Insert);
1274 
1275  Write(GroupName, "Values", LinearMap.NumMyElements() * x.RowSize(),
1276  LinearMap.NumGlobalElements() * x.RowSize(),
1277  H5T_NATIVE_INT, LinearX.Values());
1278  }
1279 
1280  Write(GroupName, "__type__", "EpetraExt::DistArray<int>");
1281  Write(GroupName, "GlobalLength", x.GlobalLength());
1282  Write(GroupName, "RowSize", x.RowSize());
1283 }
1284 
1285 // ==========================================================================
1286 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
1287  DistArray<int>*& X)
1288 {
1289  // gets the length of the std::vector
1290  int GlobalLength, RowSize;
1291  ReadIntDistArrayProperties(GroupName, GlobalLength, RowSize);
1292 
1293  if (Map.LinearMap())
1294  {
1295  X = new EpetraExt::DistArray<int>(Map, RowSize);
1296  // simply read stuff and go home
1297  Read(GroupName, "Values", Map.NumMyElements() * RowSize,
1298  Map.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->Values());
1299  }
1300  else
1301  {
1302  // we need to first create a linear map, read the std::vector,
1303  // then import it to the actual nonlinear map
1304  Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
1305  EpetraExt::DistArray<int> LinearX(LinearMap, RowSize);
1306 
1307  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1308  LinearMap.NumGlobalElements() * RowSize,
1309  H5T_NATIVE_INT, LinearX.Values());
1310 
1311  Epetra_Import Importer(Map, LinearMap);
1312  X = new EpetraExt::DistArray<int>(Map, RowSize);
1313  X->Import(LinearX, Importer, Insert);
1314  }
1315 }
1316 
1317 // ==========================================================================
1318 void EpetraExt::HDF5::Read(const std::string& GroupName, DistArray<int>*& X)
1319 {
1320  // gets the length of the std::vector
1321  int GlobalLength, RowSize;
1322  ReadIntDistArrayProperties(GroupName, GlobalLength, RowSize);
1323 
1324  // creates a new linear map and use it to read data
1325  Epetra_Map LinearMap(GlobalLength, 0, Comm_);
1326  X = new EpetraExt::DistArray<int>(LinearMap, RowSize);
1327 
1328  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1329  LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->Values());
1330 }
1331 
1332 // ==========================================================================
1333 void EpetraExt::HDF5::ReadIntDistArrayProperties(const std::string& GroupName,
1334  int& GlobalLength,
1335  int& RowSize)
1336 {
1337  if (!IsContained(GroupName))
1338  throw(Exception(__FILE__, __LINE__,
1339  "requested group " + GroupName + " not found"));
1340 
1341  std::string Label;
1342  Read(GroupName, "__type__", Label);
1343 
1344  if (Label != "EpetraExt::DistArray<int>")
1345  throw(Exception(__FILE__, __LINE__,
1346  "requested group " + GroupName + " is not an EpetraExt::DistArray<int>",
1347  "__type__ = " + Label));
1348 
1349  Read(GroupName, "GlobalLength", GlobalLength);
1350  Read(GroupName, "RowSize", RowSize);
1351 }
1352 
1353 // ============================ //
1354 // EpetraExt::DistArray<double> //
1355 // ============================ //
1356 
1357 // ==========================================================================
1358 void EpetraExt::HDF5::Write(const std::string& GroupName, const DistArray<double>& x)
1359 {
1360  if (x.Map().LinearMap())
1361  {
1362  Write(GroupName, "Values", x.Map().NumMyElements() * x.RowSize(),
1363  x.Map().NumGlobalElements() * x.RowSize(),
1364  H5T_NATIVE_DOUBLE, x.Values());
1365  }
1366  else
1367  {
1368  // need to build a linear map first, the import data, then
1369  // finally write them
1370  const Epetra_BlockMap& OriginalMap = x.Map();
1371  Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
1372  Epetra_Import Importer(LinearMap, OriginalMap);
1373 
1374  EpetraExt::DistArray<double> LinearX(LinearMap, x.RowSize());
1375  LinearX.Import(x, Importer, Insert);
1376 
1377  Write(GroupName, "Values", LinearMap.NumMyElements() * x.RowSize(),
1378  LinearMap.NumGlobalElements() * x.RowSize(),
1379  H5T_NATIVE_DOUBLE, LinearX.Values());
1380  }
1381 
1382  Write(GroupName, "__type__", "EpetraExt::DistArray<double>");
1383  Write(GroupName, "GlobalLength", x.GlobalLength());
1384  Write(GroupName, "RowSize", x.RowSize());
1385 }
1386 
1387 // ==========================================================================
1388 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
1389  DistArray<double>*& X)
1390 {
1391  // gets the length of the std::vector
1392  int GlobalLength, RowSize;
1393  ReadDoubleDistArrayProperties(GroupName, GlobalLength, RowSize);
1394 
1395  if (Map.LinearMap())
1396  {
1397  X = new EpetraExt::DistArray<double>(Map, RowSize);
1398  // simply read stuff and go home
1399  Read(GroupName, "Values", Map.NumMyElements() * RowSize,
1400  Map.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->Values());
1401  }
1402  else
1403  {
1404  // we need to first create a linear map, read the std::vector,
1405  // then import it to the actual nonlinear map
1406  Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
1407  EpetraExt::DistArray<double> LinearX(LinearMap, RowSize);
1408 
1409  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1410  LinearMap.NumGlobalElements() * RowSize,
1411  H5T_NATIVE_DOUBLE, LinearX.Values());
1412 
1413  Epetra_Import Importer(Map, LinearMap);
1414  X = new EpetraExt::DistArray<double>(Map, RowSize);
1415  X->Import(LinearX, Importer, Insert);
1416  }
1417 }
1418 
1419 // ==========================================================================
1420 void EpetraExt::HDF5::Read(const std::string& GroupName, DistArray<double>*& X)
1421 {
1422  // gets the length of the std::vector
1423  int GlobalLength, RowSize;
1424  ReadDoubleDistArrayProperties(GroupName, GlobalLength, RowSize);
1425 
1426  // creates a new linear map and use it to read data
1427  Epetra_Map LinearMap(GlobalLength, 0, Comm_);
1428  X = new EpetraExt::DistArray<double>(LinearMap, RowSize);
1429 
1430  Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1431  LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->Values());
1432 }
1433 //
1434 // ==========================================================================
1435 void EpetraExt::HDF5::ReadDoubleDistArrayProperties(const std::string& GroupName,
1436  int& GlobalLength,
1437  int& RowSize)
1438 {
1439  if (!IsContained(GroupName))
1440  throw(Exception(__FILE__, __LINE__,
1441  "requested group " + GroupName + " not found"));
1442 
1443  std::string Label;
1444  Read(GroupName, "__type__", Label);
1445 
1446  if (Label != "EpetraExt::DistArray<double>")
1447  throw(Exception(__FILE__, __LINE__,
1448  "requested group " + GroupName + " is not an EpetraExt::DistArray<double>",
1449  "__type__ = " + Label));
1450 
1451  Read(GroupName, "GlobalLength", GlobalLength);
1452  Read(GroupName, "RowSize", RowSize);
1453 }
1454 
1455 // ======================= //
1456 // basic serial data types //
1457 // ======================= //
1458 
1459 // ==========================================================================
1460 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1461  int what)
1462 {
1463  if (!IsContained(GroupName))
1464  CreateGroup(GroupName);
1465 
1466  hid_t filespace_id = H5Screate(H5S_SCALAR);
1467  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1468  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_INT,
1469  filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1470 
1471  herr_t status = H5Dwrite(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
1472  H5P_DEFAULT, &what);
1473  CHECK_STATUS(status);
1474 
1475  // Close/release resources.
1476  H5Dclose(dset_id);
1477  H5Gclose(group_id);
1478  H5Sclose(filespace_id);
1479 }
1480 
1481 // ==========================================================================
1482 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1483  double what)
1484 {
1485  if (!IsContained(GroupName))
1486  CreateGroup(GroupName);
1487 
1488  hid_t filespace_id = H5Screate(H5S_SCALAR);
1489  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1490  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_DOUBLE,
1491  filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1492 
1493  herr_t status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL,
1494  filespace_id, H5P_DEFAULT, &what);
1495  CHECK_STATUS(status);
1496 
1497  // Close/release resources.
1498  H5Dclose(dset_id);
1499  H5Gclose(group_id);
1500  H5Sclose(filespace_id);
1501 }
1502 
1503 // ==========================================================================
1504 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName, int& data)
1505 {
1506  if (!IsContained(GroupName))
1507  throw(Exception(__FILE__, __LINE__,
1508  "requested group " + GroupName + " not found"));
1509 
1510  // Create group in the root group using absolute name.
1511  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1512 
1513  hid_t filespace_id = H5Screate(H5S_SCALAR);
1514  hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1515 
1516  herr_t status = H5Dread(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
1517  H5P_DEFAULT, &data);
1518  CHECK_STATUS(status);
1519 
1520  H5Sclose(filespace_id);
1521  H5Dclose(dset_id);
1522  H5Gclose(group_id);
1523 }
1524 
1525 // ==========================================================================
1526 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName, double& data)
1527 {
1528  if (!IsContained(GroupName))
1529  throw(Exception(__FILE__, __LINE__,
1530  "requested group " + GroupName + " not found"));
1531 
1532  // Create group in the root group using absolute name.
1533  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1534 
1535  hid_t filespace_id = H5Screate(H5S_SCALAR);
1536  hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1537 
1538  herr_t status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, filespace_id,
1539  H5P_DEFAULT, &data);
1540  CHECK_STATUS(status);
1541 
1542  H5Sclose(filespace_id);
1543  H5Dclose(dset_id);
1544  H5Gclose(group_id);
1545 }
1546 
1547 // ==========================================================================
1548 void EpetraExt::HDF5::Write(const std::string& GroupName,
1549  const std::string& DataSetName,
1550  const std::string& data)
1551 {
1552  if (!IsContained(GroupName))
1553  CreateGroup(GroupName);
1554 
1555  hsize_t len = 1;
1556 
1557  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1558 
1559  hid_t dataspace_id = H5Screate_simple(1, &len, NULL);
1560 
1561  hid_t atype = H5Tcopy(H5T_C_S1);
1562  H5Tset_size(atype, data.size() + 1);
1563 
1564  hid_t dataset_id = H5Dcreate(group_id, DataSetName.c_str(), atype,
1565  dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1566  CHECK_STATUS(H5Dwrite(dataset_id, atype, H5S_ALL, H5S_ALL,
1567  H5P_DEFAULT, data.c_str()));
1568 
1569  CHECK_STATUS(H5Tclose(atype));
1570 
1571  CHECK_STATUS(H5Dclose(dataset_id));
1572 
1573  CHECK_STATUS(H5Sclose(dataspace_id));
1574 
1575  CHECK_STATUS(H5Gclose(group_id));
1576 }
1577 
1578 // ==========================================================================
1579 void EpetraExt::HDF5::Read(const std::string& GroupName,
1580  const std::string& DataSetName,
1581  std::string& data)
1582 {
1583  if (!IsContained(GroupName))
1584  throw(Exception(__FILE__, __LINE__,
1585  "requested group " + GroupName + " not found"));
1586 
1587  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1588 
1589  hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1590 
1591  hid_t datatype_id = H5Dget_type(dataset_id);
1592 // size_t typesize_id = H5Tget_size(datatype_id);
1593  H5T_class_t typeclass_id = H5Tget_class(datatype_id);
1594 
1595  if(typeclass_id != H5T_STRING)
1596  throw(Exception(__FILE__, __LINE__,
1597  "requested group " + GroupName + " is not a std::string"));
1598  char data2[160];
1599  CHECK_STATUS(H5Dread(dataset_id, datatype_id, H5S_ALL, H5S_ALL,
1600  H5P_DEFAULT, data2));
1601  data = data2;
1602 
1603  H5Dclose(dataset_id);
1604  H5Gclose(group_id);
1605 }
1606 
1607 // ============= //
1608 // serial arrays //
1609 // ============= //
1610 
1611 // ==========================================================================
1612 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1613  const hid_t type, const int Length,
1614  const void* data)
1615 {
1616  if (!IsContained(GroupName))
1617  CreateGroup(GroupName);
1618 
1619  hsize_t dimsf = Length;
1620 
1621  hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
1622 
1623  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1624 
1625  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type,
1626  filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1627 
1628  herr_t status = H5Dwrite(dset_id, type, H5S_ALL, H5S_ALL,
1629  H5P_DEFAULT, data);
1630  CHECK_STATUS(status);
1631 
1632  // Close/release resources.
1633  H5Dclose(dset_id);
1634  H5Gclose(group_id);
1635  H5Sclose(filespace_id);
1636 }
1637 
1638 // ==========================================================================
1639 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName,
1640  const hid_t type, const int Length,
1641  void* data)
1642 {
1643  if (!IsContained(GroupName))
1644  throw(Exception(__FILE__, __LINE__,
1645  "requested group " + GroupName + " not found"));
1646 
1647  hsize_t dimsf = Length;
1648 
1649  hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
1650 
1651  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1652 
1653  hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1654 
1655  herr_t status = H5Dread(dset_id, type, H5S_ALL, filespace_id,
1656  H5P_DEFAULT, data);
1657  CHECK_STATUS(status);
1658 
1659  // Close/release resources.
1660  H5Dclose(dset_id);
1661  H5Gclose(group_id);
1662  H5Sclose(filespace_id);
1663 }
1664 
1665 // ================== //
1666 // distributed arrays //
1667 // ================== //
1668 
1669 // ==========================================================================
1670 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
1671  int MySize, int GlobalSize, hid_t type, const void* data)
1672 {
1673  int Offset;
1674  Comm_.ScanSum(&MySize, &Offset, 1);
1675  Offset -= MySize;
1676 
1677  hsize_t MySize_t = MySize;
1678  hsize_t GlobalSize_t = GlobalSize;
1679  hsize_t Offset_t = Offset;
1680 
1681  hid_t filespace_id = H5Screate_simple(1, &GlobalSize_t, NULL);
1682 
1683  // Create the dataset with default properties and close filespace.
1684  if (!IsContained(GroupName))
1685  CreateGroup(GroupName);
1686 
1687  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1688  hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type, filespace_id,
1689  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1690 
1691  H5Sclose(filespace_id);
1692 
1693  // Each process defines dataset in memory and writes it to the hyperslab
1694  // in the file.
1695 
1696  hid_t memspace_id = H5Screate_simple(1, &MySize_t, NULL);
1697 
1698  // Select hyperslab in the file.
1699  filespace_id = H5Dget_space(dset_id);
1700  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL, &MySize_t, NULL);
1701 
1702  plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1703 #ifdef HAVE_MPI
1704  H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1705 #endif
1706 
1707  status = H5Dwrite(dset_id, type, memspace_id, filespace_id,
1708  plist_id_, data);
1709  CHECK_STATUS(status);
1710 
1711  // Close/release resources.
1712  H5Dclose(dset_id);
1713  H5Gclose(group_id);
1714  H5Sclose(filespace_id);
1715  H5Sclose(memspace_id);
1716  H5Pclose(plist_id_);
1717 }
1718 
1719 // ==========================================================================
1720 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName,
1721  int MySize, int GlobalSize,
1722  const hid_t type, void* data)
1723 {
1724  if (!IsOpen())
1725  throw(Exception(__FILE__, __LINE__, "no file open yet"));
1726 
1727  hsize_t MySize_t = MySize;
1728 
1729  // offset
1730  int itmp;
1731  Comm_.ScanSum(&MySize, &itmp, 1);
1732  hsize_t Offset_t = itmp - MySize;
1733 
1734  hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1735  hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1736  //hid_t space_id = H5Screate_simple(1, &Offset_t, 0);
1737 
1738  // Select hyperslab in the file.
1739  hid_t filespace_id = H5Dget_space(dataset_id);
1740  H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL,
1741  &MySize_t, NULL);
1742 
1743  hid_t mem_dataspace = H5Screate_simple (1, &MySize_t, NULL);
1744 
1745  herr_t status = H5Dread(dataset_id, type, mem_dataspace, filespace_id,
1746  H5P_DEFAULT, data);
1747  CHECK_STATUS(status);
1748 
1749  H5Sclose(mem_dataspace);
1750  H5Gclose(group_id);
1751  //H5Sclose(space_id);
1752  H5Dclose(dataset_id);
1753 // H5Dclose(filespace_id);
1754 }
1755 
1756 
1757 #endif
1758 
Epetra_MultiVector
EpetraExt::HDF5::ReadDoubleDistArrayProperties
void ReadDoubleDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
Definition: EpetraExt_HDF5.cpp:1435
Epetra_BlockMap::NumGlobalElements
int NumGlobalElements() const
Epetra_RowMatrix.h
Teuchos::ParameterList::name
const std::string & name() const
Epetra_CrsGraph::NumGlobalRows
int NumGlobalRows() const
EpetraExt::HDF5::Write
void Write(const std::string &GroupName, const std::string &DataSetName, int data)
Write an integer in group GroupName using the given DataSetName.
Definition: EpetraExt_HDF5.cpp:1460
Epetra_Map
EpetraExt_Exception.h
Epetra_FECrsGraph::InsertGlobalIndices
int InsertGlobalIndices(int numRows, const int *rows, int numCols, const int *cols)
Teuchos::ParameterList::end
ConstIterator end() const
Epetra_CrsGraph::FillComplete
int FillComplete()
EpetraExt::HDF5::Write
void Write(const std::string &GroupName, const Epetra_BlockMap &Map)
Write a block map to group GroupName.
Definition: EpetraExt_HDF5.cpp:428
Epetra_RowMatrix::RowMatrixRowMap
virtual const Epetra_Map & RowMatrixRowMap() const =0
Epetra_RowMatrix::NormOne
virtual double NormOne() const =0
Epetra_IntVector
Epetra_BlockMap::MyGlobalElements
int MyGlobalElements(int *MyGlobalElementList) const
Teuchos::ParameterList::get
T & get(ParameterList &l, const std::string &name)
Epetra_BlockMap.h
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
CHECK_HID
#define CHECK_HID(hid_t)
Definition: EpetraExt_HDF5.cpp:78
EpetraExt::DistArray::GlobalLength
int GlobalLength() const
Returns the global length of the array.
Definition: EpetraExt_DistArray.h:119
Epetra_FECrsGraph
EpetraExt::HDF5::ReadMapProperties
void ReadMapProperties(const std::string &GroupName, int &NumGlobalElements, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_Map.
Definition: EpetraExt_HDF5.cpp:565
Epetra_CrsGraph::Filled
bool Filled() const
WriteParameterListRecursive
static void WriteParameterListRecursive(const Teuchos::ParameterList &params, hid_t group_id)
Definition: EpetraExt_HDF5.cpp:107
EpetraExt::HDF5::IsContained
bool IsContained(std::string Name, std::string GroupName="")
Return true if Name is contained in the database.
Definition: EpetraExt_HDF5.cpp:395
Epetra_CrsGraph::NumMyRows
int NumMyRows() const
Epetra_BlockMap::ElementSizeList
int * ElementSizeList() const
Epetra_Import.h
EpetraExt::HDF5::Read
void Read(const std::string &GroupName, const std::string &DataSetName, int &data)
Read an integer from group /GroupName/DataSetName.
Definition: EpetraExt_HDF5.cpp:1504
Epetra_CrsGraph.h
Insert
Teuchos::ParameterList::sublist
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Epetra_MpiComm
Epetra_Import
Epetra_IntVector::Values
int * Values() const
Epetra_IntVector::GlobalLength
int GlobalLength() const
Epetra_SerialComm.h
Epetra_RowMatrix::NumGlobalNonzeros
virtual int NumGlobalNonzeros() const =0
Epetra_FECrsMatrix::InsertGlobalValues
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
EpetraExt_ConfigDefs.h
Epetra_BlockMap::IndexBase
int IndexBase() const
EpetraExt_DistArray.h
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Epetra_RowMatrix::MaxNumEntries
virtual int MaxNumEntries() const =0
Epetra_FECrsMatrix
Epetra_CrsGraph::NumGlobalCols
int NumGlobalCols() const
EpetraExt::DistArray::Values
T * Values()
Returns a pointer to the internally stored data (non-const version).
Definition: EpetraExt_DistArray.h:198
Epetra_Comm
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos_ParameterList.hpp
EpetraExt::HDF5::ReadBlockMapProperties
void ReadBlockMapProperties(const std::string &GroupName, int &NumGlobalElements, int &NumGlobalPoints, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_BlockMap.
Definition: EpetraExt_HDF5.cpp:495
EpetraExt::HDF5::ReadIntVectorProperties
void ReadIntVectorProperties(const std::string &GroupName, int &GlobalLength)
Read basic properties of specified Epetra_IntVector.
Definition: EpetraExt_HDF5.cpp:667
EpetraExt::HDF5::ReadCrsGraphProperties
void ReadCrsGraphProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumGlobalNonzeros, int &NumGlobalDiagonals, int &MaxNumIndices)
Read basic properties of specified Epetra_CrsGraph.
Definition: EpetraExt_HDF5.cpp:744
Teuchos::ParameterList::isSublist
bool isSublist(const std::string &name) const
Epetra_BlockMap::GID
int GID(int LID) const
Epetra_SrcDistObject::Map
virtual const Epetra_BlockMap & Map() const =0
Epetra_RowMatrix::NormInf
virtual double NormInf() const =0
Epetra_RowMatrix::NumMyRows
virtual int NumMyRows() const =0
Epetra_BlockMap
Epetra_FECrsMatrix::GlobalAssemble
int GlobalAssemble(bool callFillComplete=true, Epetra_CombineMode combineMode=Add, bool save_off_and_reuse_map_exporter=false)
Epetra_CrsGraph::GCID
int GCID(int LCID_in) const
Epetra_RowMatrix::Filled
virtual bool Filled() const =0
Epetra_CrsGraph::NumGlobalNonzeros
int NumGlobalNonzeros() const
Epetra_MpiComm::Comm
MPI_Comm Comm() const
Epetra_RowMatrix::NumGlobalDiagonals
virtual int NumGlobalDiagonals() const =0
Epetra_RowMatrix::NumGlobalCols
virtual int NumGlobalCols() const =0
Teuchos::ParameterList::begin
ConstIterator begin() const
EpetraExt::DistArray::RowSize
int RowSize() const
Returns the row size, that is, the amount of data associated with each element.
Definition: EpetraExt_DistArray.h:125
EpetraExt_Utils.h
Epetra_CrsGraph::ExtractMyRowView
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Epetra_SerialComm
Epetra_BlockMap::LinearMap
bool LinearMap() const
EpetraExt::HDF5::Open
void Open(const std::string FileName, int AccessType=H5F_ACC_RDWR)
Open specified file with given access type.
Definition: EpetraExt_HDF5.cpp:365
Epetra_MpiComm.h
Epetra_IntVector.h
CHECK_STATUS
#define CHECK_STATUS(status)
Definition: EpetraExt_HDF5.cpp:83
Epetra_CrsGraph::GRID
int GRID(int LRID_in) const
f_operator
static herr_t f_operator(hid_t loc_id, const char *name, void *opdata)
Definition: EpetraExt_HDF5.cpp:210
Epetra_FECrsGraph.h
Epetra_CrsMatrix
Epetra_CrsGraph::MaxNumIndices
int MaxNumIndices() const
Epetra_Map.h
Epetra_MultiVector.h
Teuchos::ParameterList::isType
bool isType(const std::string &name) const
EpetraExt::HDF5::HDF5
HDF5(const Epetra_Comm &Comm)
Constructor.
Definition: EpetraExt_HDF5.cpp:272
Epetra_RowMatrix::ExtractMyRowCopy
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
EpetraExt::HDF5::Create
void Create(const std::string FileName)
Create a new file.
Definition: EpetraExt_HDF5.cpp:278
Epetra_RowMatrix::RowMatrixColMap
virtual const Epetra_Map & RowMatrixColMap() const =0
Epetra_BlockMap::NumGlobalPoints
int NumGlobalPoints() const
EpetraExt::HDF5::ReadIntDistArrayProperties
void ReadIntDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
Definition: EpetraExt_HDF5.cpp:1333
Epetra_BlockMap::NumMyPoints
int NumMyPoints() const
Epetra_RowMatrix
EpetraExt::HDF5::ReadCrsMatrixProperties
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.
Definition: EpetraExt_HDF5.cpp:928
Epetra_CrsMatrix.h
FindDataset
static herr_t FindDataset(hid_t loc_id, const char *name, void *opdata)
Definition: EpetraExt_HDF5.cpp:96
Epetra_CrsGraph::NumMyNonzeros
int NumMyNonzeros() const
FindDataset_t::found
bool found
Definition: EpetraExt_HDF5.cpp:93
Epetra_FECrsMatrix.h
n
int n
Teuchos_RCP.hpp
EpetraExt_HDF5.h
Epetra_RowMatrix::NumGlobalRows
virtual int NumGlobalRows() const =0
FindDataset_t::name
std::string name
Definition: EpetraExt_HDF5.cpp:92
EpetraExt::HDF5::ReadMultiVectorProperties
void ReadMultiVectorProperties(const std::string &GroupName, int &GlobalLength, int &NumVectors)
Read basic properties of specified Epetra_MultiVector.
Definition: EpetraExt_HDF5.cpp:1231
Epetra_RowMatrix::NumMyNonzeros
virtual int NumMyNonzeros() const =0
EpetraExt::toString
std::string toString(const int &x)
Definition: EpetraExt_Utils.cpp:48
Epetra_CrsGraph
Epetra_CrsGraph::NumGlobalDiagonals
int NumGlobalDiagonals() const
FindDataset_t
Definition: EpetraExt_HDF5.cpp:90
EpetraExt::DistArray
DistArray&lt;T&gt;: A class to store row-oriented multivectors of type T.
Definition: EpetraExt_DistArray.h:85
Generated by   doxygen 1.8.5