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
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
90struct FindDataset_t
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// ==========================================================================
272EpetraExt::HDF5::HDF5(const Epetra_Comm& Comm) :
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); // CHECK: ALLOW MPI_COMM_WORLD
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,
635 Epetra_IntVector*& X)
636{
637 // gets the length of the std::vector
638 int GlobalLength;
639 ReadIntVectorProperties(GroupName, GlobalLength);
640
641 if (Map.LinearMap())
642 {
643 X = new Epetra_IntVector(Map);
644 // simply read stuff and go home
645 Read(GroupName, "Values", Map.NumMyElements(), Map.NumGlobalElements(),
646 H5T_NATIVE_INT, X->Values());
647
648 }
649 else
650 {
651 // we need to first create a linear map, read the std::vector,
652 // then import it to the actual nonlinear map
653 Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
654 Epetra_IntVector LinearX(LinearMap);
655
656 Read(GroupName, "Values", LinearMap.NumMyElements(),
657 LinearMap.NumGlobalElements(),
658 H5T_NATIVE_INT, LinearX.Values());
659
660 Epetra_Import Importer(Map, LinearMap);
661 X = new Epetra_IntVector(Map);
662 X->Import(LinearX, Importer, Insert);
663 }
664}
665
666// ==========================================================================
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,
1389 DistArray<double>*& X)
1390{
1391 // gets the length of the std::vector
1392 int GlobalLength, RowSize;
1393 ReadDoubleDistArrayProperties(GroupName, GlobalLength, RowSize);
1394
1395 if (Map.LinearMap())
1396 {
1397 X = new EpetraExt::DistArray<double>(Map, RowSize);
1398 // simply read stuff and go home
1399 Read(GroupName, "Values", Map.NumMyElements() * RowSize,
1400 Map.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->Values());
1401 }
1402 else
1403 {
1404 // we need to first create a linear map, read the std::vector,
1405 // then import it to the actual nonlinear map
1406 Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
1407 EpetraExt::DistArray<double> LinearX(LinearMap, RowSize);
1408
1409 Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize,
1410 LinearMap.NumGlobalElements() * RowSize,
1411 H5T_NATIVE_DOUBLE, LinearX.Values());
1412
1413 Epetra_Import Importer(Map, LinearMap);
1414 X = new EpetraExt::DistArray<double>(Map, RowSize);
1415 X->Import(LinearX, Importer, Insert);
1416 }
1417}
1418
1419// ==========================================================================
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
EpetraExt_ConfigDefs.h
EpetraExt_DistArray.h
EpetraExt_Exception.h
CHECK_STATUS
#define CHECK_STATUS(status)
Definition EpetraExt_HDF5.cpp:83
CHECK_HID
#define CHECK_HID(hid_t)
Definition EpetraExt_HDF5.cpp:78
WriteParameterListRecursive
static void WriteParameterListRecursive(const Teuchos::ParameterList &params, hid_t group_id)
Definition EpetraExt_HDF5.cpp:107
f_operator
static herr_t f_operator(hid_t loc_id, const char *name, void *opdata)
Definition EpetraExt_HDF5.cpp:210
FindDataset
static herr_t FindDataset(hid_t loc_id, const char *name, void *opdata)
Definition EpetraExt_HDF5.cpp:96
EpetraExt_HDF5.h
EpetraExt_Utils.h
Insert
Insert
Copy
Copy
EpetraExt::DistArray
DistArray<T>: A class to store row-oriented multivectors of type T.
Definition EpetraExt_DistArray.h:92
EpetraExt::DistArray::Values
T * Values()
Returns a pointer to the internally stored data (non-const version).
Definition EpetraExt_DistArray.h:204
EpetraExt::DistArray::RowSize
int RowSize() const
Returns the row size, that is, the amount of data associated with each element.
Definition EpetraExt_DistArray.h:131
EpetraExt::DistArray::GlobalLength
int GlobalLength() const
Returns the global length of the array.
Definition EpetraExt_DistArray.h:125
EpetraExt::HDF5::HDF5
HDF5(const Epetra_Comm &Comm)
Constructor.
Definition EpetraExt_HDF5.cpp:272
EpetraExt::HDF5::ReadMultiVectorProperties
void ReadMultiVectorProperties(const std::string &GroupName, int &GlobalLength, int &NumVectors)
Read basic properties of specified Epetra_MultiVector.
Definition EpetraExt_HDF5.cpp:1231
EpetraExt::HDF5::IsContained
bool IsContained(std::string Name, std::string GroupName="")
Return true if Name is contained in the database.
Definition EpetraExt_HDF5.cpp:395
EpetraExt::HDF5::Write
void Write(const std::string &GroupName, const std::string &DataSetName, int data)
Write an integer in group GroupName using the given DataSetName.
Definition EpetraExt_HDF5.cpp:1460
EpetraExt::HDF5::Open
void Open(const std::string FileName, int AccessType=H5F_ACC_RDWR)
Open specified file with given access type.
Definition EpetraExt_HDF5.cpp:365
EpetraExt::HDF5::ReadCrsGraphProperties
void ReadCrsGraphProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumGlobalNonzeros, int &NumGlobalDiagonals, int &MaxNumIndices)
Read basic properties of specified Epetra_CrsGraph.
Definition EpetraExt_HDF5.cpp:744
EpetraExt::HDF5::ReadCrsMatrixProperties
void ReadCrsMatrixProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumNonzeros, int &NumGlobalDiagonals, int &MaxNumEntries, double &NormOne, double &NormInf)
Read basic properties of specified Epetra_CrsMatrix.
Definition EpetraExt_HDF5.cpp:928
EpetraExt::HDF5::ReadIntDistArrayProperties
void ReadIntDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
Definition EpetraExt_HDF5.cpp:1333
EpetraExt::HDF5::ReadDoubleDistArrayProperties
void ReadDoubleDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
Definition EpetraExt_HDF5.cpp:1435
EpetraExt::HDF5::ReadIntVectorProperties
void ReadIntVectorProperties(const std::string &GroupName, int &GlobalLength)
Read basic properties of specified Epetra_IntVector.
Definition EpetraExt_HDF5.cpp:667
EpetraExt::HDF5::Read
void Read(const std::string &GroupName, const std::string &DataSetName, int &data)
Read an integer from group /GroupName/DataSetName.
Definition EpetraExt_HDF5.cpp:1504
EpetraExt::HDF5::Create
void Create(const std::string FileName)
Create a new file.
Definition EpetraExt_HDF5.cpp:278
EpetraExt::HDF5::ReadMapProperties
void ReadMapProperties(const std::string &GroupName, int &NumGlobalElements, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_Map.
Definition EpetraExt_HDF5.cpp:565
EpetraExt::HDF5::ReadBlockMapProperties
void ReadBlockMapProperties(const std::string &GroupName, int &NumGlobalElements, int &NumGlobalPoints, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_BlockMap.
Definition EpetraExt_HDF5.cpp:495
Epetra_BlockMap
Epetra_BlockMap::MyGlobalElements
int MyGlobalElements(int *MyGlobalElementList) const
Epetra_BlockMap::ElementSizeList
int * ElementSizeList() const
Epetra_BlockMap::NumMyPoints
int NumMyPoints() const
Epetra_BlockMap::GID
int GID(int LID) const
Epetra_BlockMap::LinearMap
bool LinearMap() const
Epetra_BlockMap::IndexBase
int IndexBase() const
Epetra_BlockMap::NumGlobalElements
int NumGlobalElements() const
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Epetra_BlockMap::NumGlobalPoints
int NumGlobalPoints() const
Epetra_Comm
Epetra_CrsGraph
Epetra_CrsGraph::Filled
bool Filled() const
Epetra_CrsGraph::ExtractMyRowView
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Epetra_CrsGraph::MaxNumIndices
int MaxNumIndices() const
Epetra_CrsGraph::GCID
int GCID(int LCID_in) const
Epetra_CrsGraph::FillComplete
int FillComplete()
Epetra_CrsGraph::GRID
int GRID(int LRID_in) const
Epetra_CrsGraph::NumGlobalDiagonals
int NumGlobalDiagonals() const
Epetra_CrsGraph::NumMyRows
int NumMyRows() const
Epetra_CrsGraph::NumGlobalNonzeros
int NumGlobalNonzeros() const
Epetra_CrsGraph::NumGlobalRows
int NumGlobalRows() const
Epetra_CrsGraph::NumMyNonzeros
int NumMyNonzeros() const
Epetra_CrsGraph::NumGlobalCols
int NumGlobalCols() const
Epetra_CrsMatrix
Epetra_DistObject::Map
const Epetra_BlockMap & Map() const
Epetra_DistObject::Import
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Epetra_FECrsGraph
Epetra_FECrsGraph::InsertGlobalIndices
int InsertGlobalIndices(int numRows, const int *rows, int numCols, const int *cols)
Epetra_FECrsMatrix
Epetra_FECrsMatrix::InsertGlobalValues
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_FECrsMatrix::GlobalAssemble
int GlobalAssemble(bool callFillComplete=true, Epetra_CombineMode combineMode=Add, bool save_off_and_reuse_map_exporter=false)
Epetra_Import
Epetra_IntVector
Epetra_IntVector::Values
int * Values() const
Epetra_IntVector::GlobalLength
int GlobalLength() const
Epetra_Map
Epetra_MpiComm
Epetra_MpiComm::Comm
MPI_Comm Comm() const
Epetra_MultiVector
Epetra_MultiVector::GlobalLength
int GlobalLength() const
Epetra_MultiVector::NumVectors
int NumVectors() const
Epetra_MultiVector::MyLength
int MyLength() const
Epetra_MultiVector::Values
double * Values() const
Epetra_RowMatrix
Epetra_RowMatrix::NumMyRows
virtual int NumMyRows() const=0
Epetra_RowMatrix::NumGlobalCols
virtual int NumGlobalCols() const=0
Epetra_RowMatrix::NumGlobalNonzeros
virtual int NumGlobalNonzeros() const=0
Epetra_RowMatrix::RowMatrixColMap
virtual const Epetra_Map & RowMatrixColMap() const=0
Epetra_RowMatrix::RowMatrixRowMap
virtual const Epetra_Map & RowMatrixRowMap() const=0
Epetra_RowMatrix::NumMyNonzeros
virtual int NumMyNonzeros() const=0
Epetra_RowMatrix::NumGlobalRows
virtual int NumGlobalRows() const=0
Epetra_RowMatrix::NumGlobalDiagonals
virtual int NumGlobalDiagonals() const=0
Epetra_RowMatrix::MaxNumEntries
virtual int MaxNumEntries() const=0
Epetra_RowMatrix::ExtractMyRowCopy
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Epetra_RowMatrix::NormOne
virtual double NormOne() const=0
Epetra_RowMatrix::NormInf
virtual double NormInf() const=0
Epetra_RowMatrix::Filled
virtual bool Filled() const=0
Epetra_SerialComm
EpetraExt::toString
std::string toString(const int &x)
Definition EpetraExt_Utils.cpp:48
FindDataset_t
Definition EpetraExt_HDF5.cpp:91
FindDataset_t::found
bool found
Definition EpetraExt_HDF5.cpp:93
FindDataset_t::name
std::string name
Definition EpetraExt_HDF5.cpp:92
Generated by doxygen 1.9.8