EpetraExt Development
Loading...
Searching...
No Matches
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
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
96static 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
107static 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
210static 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// ==========================================================================
273 Comm_(Comm),
274 IsOpen_(false)
275{}
276
277// ==========================================================================
278void 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// ==========================================================================
365void 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// ==========================================================================
395bool 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// ==========================================================================
428void 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// ==========================================================================
459void 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// ==========================================================================
495void 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// ==========================================================================
520void 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// ==========================================================================
539void 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// ==========================================================================
565void 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// ==========================================================================
592void 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// ==========================================================================
619void 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// ==========================================================================
634void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
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// ==========================================================================
667void 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// ==========================================================================
690void 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// ==========================================================================
728void 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// ==========================================================================
744void 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// ==========================================================================
771void 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// ==========================================================================
824void 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// ==========================================================================
868void 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// ==========================================================================
886void 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// ==========================================================================
928void 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// ==========================================================================
963void 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// ==========================================================================
981void 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// ==========================================================================
1017void 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// ==========================================================================
1109void 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// ==========================================================================
1126void 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// ==========================================================================
1231void 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// ==========================================================================
1256void 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// ==========================================================================
1286void 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// ==========================================================================
1318void 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// ==========================================================================
1333void 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// ==========================================================================
1358void 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// ==========================================================================
1388void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
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// ==========================================================================
1420void 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// ==========================================================================
1435void 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// ==========================================================================
1460void 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// ==========================================================================
1482void 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// ==========================================================================
1504void 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// ==========================================================================
1526void 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// ==========================================================================
1548void 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// ==========================================================================
1579void 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// ==========================================================================
1612void 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// ==========================================================================
1639void 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// ==========================================================================
1670void 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// ==========================================================================
1720void 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
#define CHECK_STATUS(status)
#define CHECK_HID(hid_t)
static void WriteParameterListRecursive(const Teuchos::ParameterList &params, hid_t group_id)
static herr_t f_operator(hid_t loc_id, const char *name, void *opdata)
static herr_t FindDataset(hid_t loc_id, const char *name, void *opdata)
Insert
Copy
DistArray<T>: A class to store row-oriented multivectors of type T.
T * Values()
Returns a pointer to the internally stored data (non-const version).
int RowSize() const
Returns the row size, that is, the amount of data associated with each element.
int GlobalLength() const
Returns the global length of the array.
void Write(const std::string &GroupName, const Epetra_BlockMap &Map)
Write a block map to group GroupName.
HDF5(const Epetra_Comm &Comm)
Constructor.
void ReadMultiVectorProperties(const std::string &GroupName, int &GlobalLength, int &NumVectors)
Read basic properties of specified Epetra_MultiVector.
bool IsContained(std::string Name, std::string GroupName="")
Return true if Name is contained in the database.
void Write(const std::string &GroupName, const std::string &DataSetName, int data)
Write an integer in group GroupName using the given DataSetName.
void Open(const std::string FileName, int AccessType=H5F_ACC_RDWR)
Open specified file with given access type.
void ReadCrsGraphProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumGlobalNonzeros, int &NumGlobalDiagonals, int &MaxNumIndices)
Read basic properties of specified Epetra_CrsGraph.
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.
void ReadIntDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
void ReadDoubleDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
void ReadIntVectorProperties(const std::string &GroupName, int &GlobalLength)
Read basic properties of specified Epetra_IntVector.
void Read(const std::string &GroupName, const std::string &DataSetName, int &data)
Read an integer from group /GroupName/DataSetName.
void Create(const std::string FileName)
Create a new file.
void ReadMapProperties(const std::string &GroupName, int &NumGlobalElements, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_Map.
void ReadBlockMapProperties(const std::string &GroupName, int &NumGlobalElements, int &NumGlobalPoints, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_BlockMap.
int MyGlobalElements(int *MyGlobalElementList) const
int * ElementSizeList() const
int NumMyPoints() const
int GID(int LID) const
bool LinearMap() const
int IndexBase() const
int NumGlobalElements() const
int NumMyElements() const
int NumGlobalPoints() const
bool Filled() const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
int MaxNumIndices() const
int GCID(int LCID_in) const
int GRID(int LRID_in) const
int NumGlobalDiagonals() const
int NumMyRows() const
int NumGlobalNonzeros() const
int NumGlobalRows() const
int NumMyNonzeros() const
int NumGlobalCols() const
const Epetra_BlockMap & Map() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int InsertGlobalIndices(int numRows, const int *rows, int numCols, const int *cols)
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int GlobalAssemble(bool callFillComplete=true, Epetra_CombineMode combineMode=Add, bool save_off_and_reuse_map_exporter=false)
int * Values() const
int GlobalLength() const
MPI_Comm Comm() const
int GlobalLength() const
int NumVectors() const
int MyLength() const
double * Values() const
virtual int NumMyRows() const=0
virtual int NumGlobalCols() const=0
virtual int NumGlobalNonzeros() const=0
virtual const Epetra_Map & RowMatrixColMap() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int NumMyNonzeros() const=0
virtual int NumGlobalRows() const=0
virtual int NumGlobalDiagonals() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
virtual double NormOne() const=0
virtual double NormInf() const=0
virtual bool Filled() const=0
std::string toString(const int &x)
std::string name