EpetraExt  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)
int GlobalLength() const
Returns the global length of the array.
void ReadMapProperties(const std::string &GroupName, int &NumGlobalElements, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_Map.
void CreateGroup(const std::string &GroupName)
Create group GroupName.
static void WriteParameterListRecursive(const Teuchos::ParameterList &params, hid_t group_id)
bool IsOpen() const
Return true if a file has already been opened using Open()/Create()
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).
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 RowSize() const
Returns the row size, that is, the amount of data associated with each element.
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.