47 #ifdef HAVE_EPETRAEXT_HDF5 52 # include <H5FDmpio.h> 53 # include "Epetra_MpiComm.h" 60 #include "Epetra_SerialComm.h" 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" 78 #define CHECK_HID(hid_t) \ 80 throw(EpetraExt::Exception(__FILE__, __LINE__, \ 81 "hid_t is negative")); } 83 #define CHECK_STATUS(status) \ 85 throw(EpetraExt::Exception(__FILE__, __LINE__, \ 86 "function H5Giterater returned a negative value")); } 110 Teuchos::ParameterList::ConstIterator it = params.begin();
111 for (; it != params.end(); ++ it)
113 std::string key(params.name(it));
114 if (params.isSublist(key))
119 hid_t child_group_id = H5Gcreate(group_id, key.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
121 H5Gclose(child_group_id);
132 hid_t dataspace_id, dataset_id;
136 if (params.isType<std::string>(key))
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());
145 status = H5Dclose(dataset_id);
147 status = H5Sclose(dataspace_id);
152 if (params.isType<
bool>(key))
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);
161 status = H5Dclose(dataset_id);
163 status = H5Sclose(dataspace_id);
168 if (params.isType<
int>(key))
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);
176 status = H5Dclose(dataset_id);
178 status = H5Sclose(dataspace_id);
183 if (params.isType<
double>(key))
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);
191 status = H5Dclose(dataset_id);
193 status = H5Sclose(dataspace_id);
201 "type for parameter " + key +
" not supported"));
213 hid_t dataset_id, space_id, type_id;
214 Teuchos::ParameterList* sublist;
215 Teuchos::ParameterList* params = (Teuchos::ParameterList*)opdata;
221 H5Gget_objinfo(loc_id, name, 0, &statbuf);
222 if (strcmp(name,
"__type__") == 0)
225 switch (statbuf.type) {
227 sublist = ¶ms->sublist(name);
228 H5Giterate(loc_id, name , NULL,
f_operator, sublist);
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)
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) {
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) {
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));
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);
258 "unsupported datatype"));
262 H5Dclose(dataset_id);
266 "unsupported datatype"));
282 "an HDF5 is already open, first close the current one",
283 "using method Close(), then open/create a new one"));
285 FileName_ = FileName;
288 plist_id_ = H5Pcreate(H5P_FILE_ACCESS);
298 MPI_Comm mpiComm = MPI_COMM_NULL;
301 const Epetra_MpiComm* mpiWrapper =
302 dynamic_cast<const Epetra_MpiComm*
> (&Comm_);
303 if (mpiWrapper != NULL) {
304 mpiComm = mpiWrapper->Comm();
308 const Epetra_SerialComm* serialWrapper =
309 dynamic_cast<const Epetra_SerialComm*
> (&Comm_);
311 if (serialWrapper != NULL) {
317 mpiComm = MPI_COMM_SELF;
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 " 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 " 347 H5Pset_fapl_mpio(plist_id_, mpiComm, MPI_INFO_NULL);
352 unsigned int boh = H5Z_FILTER_MAX;
353 H5Pset_filter(plist_id_, H5Z_FILTER_DEFLATE, H5Z_FILTER_MAX, 0, &boh);
357 file_id_ = H5Fcreate(FileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT,
369 "an HDF5 is already open, first close the current one",
370 "using method Close(), then open/create a new one"));
372 FileName_ = FileName;
375 plist_id_ = H5Pcreate(H5P_FILE_ACCESS);
379 const Epetra_MpiComm* MpiComm ( dynamic_cast<const Epetra_MpiComm*> (&Comm_) );
382 H5Pset_fapl_mpio(plist_id_, MPI_COMM_WORLD, MPI_INFO_NULL);
384 H5Pset_fapl_mpio(plist_id_, MpiComm->Comm(), MPI_INFO_NULL);
388 file_id_ = H5Fopen(FileName.c_str(), AccessType, plist_id_);
398 throw(
Exception(__FILE__, __LINE__,
"no file open yet"));
405 H5Giterate(file_id_,
"/", NULL,
FindDataset, (
void*)&data);
417 int NumMyPoints = BlockMap.NumMyPoints();
418 int NumMyElements = BlockMap.NumMyElements();
419 int NumGlobalElements = BlockMap.NumGlobalElements();
420 int NumGlobalPoints = BlockMap.NumGlobalPoints();
421 int* MyGlobalElements = BlockMap.MyGlobalElements();
422 int* ElementSizeList = BlockMap.ElementSizeList();
424 std::vector<int> NumMyElements_v(Comm_.NumProc());
425 Comm_.GatherAll(&NumMyElements, &NumMyElements_v[0], 1);
427 std::vector<int> NumMyPoints_v(Comm_.NumProc());
428 Comm_.GatherAll(&NumMyPoints, &NumMyPoints_v[0], 1);
430 Write(Name,
"MyGlobalElements", NumMyElements, NumGlobalElements,
431 H5T_NATIVE_INT, MyGlobalElements);
432 Write(Name,
"ElementSizeList", NumMyElements, NumGlobalElements,
433 H5T_NATIVE_INT, ElementSizeList);
434 Write(Name,
"NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
437 Write(Name,
"NumProc", Comm_.NumProc());
439 Write(Name,
"NumGlobalPoints", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalPoints);
440 Write(Name,
"NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalElements);
441 Write(Name,
"IndexBase", BlockMap.IndexBase());
442 Write(Name,
"__type__",
"Epetra_BlockMap");
448 int NumGlobalElements, NumGlobalPoints, IndexBase, NumProc;
453 std::vector<int> NumMyPoints_v(Comm_.NumProc());
454 std::vector<int> NumMyElements_v(Comm_.NumProc());
456 Read(GroupName,
"NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
457 Read(GroupName,
"NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
458 int NumMyElements = NumMyElements_v[Comm_.MyPID()];
461 if (NumProc != Comm_.NumProc())
463 "requested map not compatible with current number of",
464 "processors, " +
toString(Comm_.NumProc()) +
467 std::vector<int> MyGlobalElements(NumMyElements);
468 std::vector<int> ElementSizeList(NumMyElements);
470 Read(GroupName,
"MyGlobalElements", NumMyElements, NumGlobalElements,
471 H5T_NATIVE_INT, &MyGlobalElements[0]);
473 Read(GroupName,
"ElementSizeList", NumMyElements, NumGlobalElements,
474 H5T_NATIVE_INT, &ElementSizeList[0]);
476 BlockMap =
new Epetra_BlockMap(NumGlobalElements, NumMyElements,
477 &MyGlobalElements[0], &ElementSizeList[0],
483 int& NumGlobalElements,
484 int& NumGlobalPoints,
490 "requested group " + GroupName +
" not found"));
493 Read(GroupName,
"__type__", Label);
495 if (Label !=
"Epetra_BlockMap")
497 "requested group " + GroupName +
" is not an Epetra_BlockMap",
498 "__type__ = " + Label));
500 Read(GroupName,
"NumGlobalElements", NumGlobalElements);
501 Read(GroupName,
"NumGlobalPoints", NumGlobalPoints);
502 Read(GroupName,
"IndexBase", IndexBase);
503 Read(GroupName,
"NumProc", NumProc);
509 int MySize = Map.NumMyElements();
510 int GlobalSize = Map.NumGlobalElements();
511 int* MyGlobalElements = Map.MyGlobalElements();
513 std::vector<int> NumMyElements(Comm_.NumProc());
514 Comm_.GatherAll(&MySize, &NumMyElements[0], 1);
516 Write(Name,
"MyGlobalElements", MySize, GlobalSize,
517 H5T_NATIVE_INT, MyGlobalElements);
518 Write(Name,
"NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements[0]);
519 Write(Name,
"NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &GlobalSize);
520 Write(Name,
"NumProc", Comm_.NumProc());
521 Write(Name,
"IndexBase", Map.IndexBase());
522 Write(Name,
"__type__",
"Epetra_Map");
528 int NumGlobalElements, IndexBase, NumProc;
532 std::vector<int> NumMyElements_v(Comm_.NumProc());
534 Read(GroupName,
"NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
535 int NumMyElements = NumMyElements_v[Comm_.MyPID()];
537 if (NumProc != Comm_.NumProc())
539 "requested map not compatible with current number of",
540 "processors, " +
toString(Comm_.NumProc()) +
543 std::vector<int> MyGlobalElements(NumMyElements);
545 Read(GroupName,
"MyGlobalElements", NumMyElements, NumGlobalElements,
546 H5T_NATIVE_INT, &MyGlobalElements[0]);
548 Map =
new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], IndexBase, Comm_);
553 int& NumGlobalElements,
559 "requested group " + GroupName +
" not found"));
562 Read(GroupName,
"__type__", Label);
564 if (Label !=
"Epetra_Map")
566 "requested group " + GroupName +
" is not an Epetra_Map",
567 "__type__ = " + Label));
569 Read(GroupName,
"NumGlobalElements", NumGlobalElements);
570 Read(GroupName,
"IndexBase", IndexBase);
571 Read(GroupName,
"NumProc", NumProc);
581 if (x.Map().LinearMap())
583 Write(Name,
"GlobalLength", x.GlobalLength());
584 Write(Name,
"Values", x.Map().NumMyElements(), x.Map().NumGlobalElements(),
585 H5T_NATIVE_INT, x.Values());
591 const Epetra_BlockMap& OriginalMap = x.Map();
592 Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
593 Epetra_Import Importer(LinearMap, OriginalMap);
595 Epetra_IntVector LinearX(LinearMap);
596 LinearX.Import(x, Importer, Insert);
598 Write(Name,
"GlobalLength", x.GlobalLength());
599 Write(Name,
"Values", LinearMap.NumMyElements(), LinearMap.NumGlobalElements(),
600 H5T_NATIVE_INT, LinearX.Values());
602 Write(Name,
"__type__",
"Epetra_IntVector");
613 Epetra_Map LinearMap(GlobalLength, 0, Comm_);
614 X =
new Epetra_IntVector(LinearMap);
616 Read(GroupName,
"Values", LinearMap.NumMyElements(),
617 LinearMap.NumGlobalElements(), H5T_NATIVE_INT, X->Values());
622 Epetra_IntVector*& X)
630 X =
new Epetra_IntVector(Map);
632 Read(GroupName,
"Values", Map.NumMyElements(), Map.NumGlobalElements(),
633 H5T_NATIVE_INT, X->Values());
640 Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
641 Epetra_IntVector LinearX(LinearMap);
643 Read(GroupName,
"Values", LinearMap.NumMyElements(),
644 LinearMap.NumGlobalElements(),
645 H5T_NATIVE_INT, LinearX.Values());
647 Epetra_Import Importer(Map, LinearMap);
648 X =
new Epetra_IntVector(Map);
649 X->Import(LinearX, Importer, Insert);
659 "requested group " + GroupName +
" not found"));
662 Read(GroupName,
"__type__", Label);
664 if (Label !=
"Epetra_IntVector")
666 "requested group " + GroupName +
" is not an Epetra_IntVector",
667 "__type__ = " + Label));
669 Read(GroupName,
"GlobalLength", GlobalLength);
681 "input Epetra_CrsGraph is not FillComplete()'d"));
684 int MySize = Graph.NumMyNonzeros();
685 int GlobalSize = Graph.NumGlobalNonzeros();
686 std::vector<int> ROW(MySize);
687 std::vector<int> COL(MySize);
693 for (
int i = 0; i < Graph.NumMyRows(); ++i)
695 Graph.ExtractMyRowView(i, NumEntries, RowIndices);
696 for (
int j = 0; j < NumEntries; ++j)
698 ROW[count] = Graph.GRID(i);
699 COL[count] = Graph.GCID(RowIndices[j]);
704 Write(Name,
"ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
705 Write(Name,
"COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
706 Write(Name,
"MaxNumIndices", Graph.MaxNumIndices());
707 Write(Name,
"NumGlobalRows", Graph.NumGlobalRows());
708 Write(Name,
"NumGlobalCols", Graph.NumGlobalCols());
709 Write(Name,
"NumGlobalNonzeros", Graph.NumGlobalNonzeros());
710 Write(Name,
"NumGlobalDiagonals", Graph.NumGlobalDiagonals());
711 Write(Name,
"__type__",
"Epetra_CrsGraph");
717 int NumGlobalRows, NumGlobalCols;
718 int NumGlobalNonzeros, NumGlobalDiagonals, MaxNumIndices;
721 NumGlobalCols, NumGlobalNonzeros,
722 NumGlobalDiagonals, MaxNumIndices);
724 Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
725 Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
727 Read(GroupName, DomainMap, RangeMap, Graph);
734 int& NumGlobalNonzeros,
735 int& NumGlobalDiagonals,
740 "requested group " + GroupName +
" not found"));
743 Read(GroupName,
"__type__", Label);
745 if (Label !=
"Epetra_CrsGraph")
747 "requested group " + GroupName +
" is not an Epetra_CrsGraph",
748 "__type__ = " + Label));
750 Read(GroupName,
"NumGlobalRows", NumGlobalRows);
751 Read(GroupName,
"NumGlobalCols", NumGlobalCols);
752 Read(GroupName,
"NumGlobalNonzeros", NumGlobalNonzeros);
753 Read(GroupName,
"NumGlobalDiagonals", NumGlobalDiagonals);
754 Read(GroupName,
"MaxNumIndices", MaxNumIndices);
759 const Epetra_Map& RangeMap, Epetra_CrsGraph*& Graph)
763 "requested group " + GroupName +
" not found in database"));
766 Read(GroupName,
"__type__", Label);
768 if (Label !=
"Epetra_CrsGraph")
770 "requested group " + GroupName +
" is not an Epetra_CrsGraph",
771 "__type__ = " + Label));
773 int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
774 Read(GroupName,
"NumGlobalNonzeros", NumGlobalNonzeros);
775 Read(GroupName,
"NumGlobalRows", NumGlobalRows);
776 Read(GroupName,
"NumGlobalCols", NumGlobalCols);
779 int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
780 if (Comm_.MyPID() == 0)
781 NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
783 std::vector<int> ROW(NumMyNonzeros);
784 Read(GroupName,
"ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
786 std::vector<int> COL(NumMyNonzeros);
787 Read(GroupName,
"COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
789 Epetra_FECrsGraph* Graph2 =
new Epetra_FECrsGraph(Copy, DomainMap, 0);
791 for (
int i = 0; i < NumMyNonzeros; )
794 while (ROW[i + count] == ROW[i])
797 Graph2->InsertGlobalIndices(1, &ROW[i], count, &COL[i]);
801 Graph2->FillComplete(DomainMap, RangeMap);
813 if (!Matrix.Filled())
815 "input Epetra_RowMatrix is not FillComplete()'d"));
817 int MySize = Matrix.NumMyNonzeros();
818 int GlobalSize = Matrix.NumGlobalNonzeros();
819 std::vector<int> ROW(MySize);
820 std::vector<int> COL(MySize);
821 std::vector<double> VAL(MySize);
824 int Length = Matrix.MaxNumEntries();
825 std::vector<int> Indices(Length);
826 std::vector<double> Values(Length);
829 for (
int i = 0; i < Matrix.NumMyRows(); ++i)
831 Matrix.ExtractMyRowCopy(i, Length, NumEntries, &Values[0], &Indices[0]);
832 for (
int j = 0; j < NumEntries; ++j)
834 ROW[count] = Matrix.RowMatrixRowMap().GID(i);
835 COL[count] = Matrix.RowMatrixColMap().GID(Indices[j]);
836 VAL[count] = Values[j];
841 Write(Name,
"ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
842 Write(Name,
"COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
843 Write(Name,
"VAL", MySize, GlobalSize, H5T_NATIVE_DOUBLE, &VAL[0]);
844 Write(Name,
"NumGlobalRows", Matrix.NumGlobalRows());
845 Write(Name,
"NumGlobalCols", Matrix.NumGlobalCols());
846 Write(Name,
"NumGlobalNonzeros", Matrix.NumGlobalNonzeros());
847 Write(Name,
"NumGlobalDiagonals", Matrix.NumGlobalDiagonals());
848 Write(Name,
"MaxNumEntries", Matrix.MaxNumEntries());
849 Write(Name,
"NormOne", Matrix.NormOne());
850 Write(Name,
"NormInf", Matrix.NormInf());
851 Write(Name,
"__type__",
"Epetra_RowMatrix");
857 int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
858 int NumGlobalDiagonals, MaxNumEntries;
859 double NormOne, NormInf;
862 NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
866 Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
867 Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
869 Read(GroupName, DomainMap, RangeMap, A);
874 const Epetra_Map& RangeMap, Epetra_CrsMatrix*& A)
876 int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
877 int NumGlobalDiagonals, MaxNumEntries;
878 double NormOne, NormInf;
881 NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
884 int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
885 if (Comm_.MyPID() == 0)
886 NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
888 std::vector<int> ROW(NumMyNonzeros);
889 Read(GroupName,
"ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
891 std::vector<int> COL(NumMyNonzeros);
892 Read(GroupName,
"COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
894 std::vector<double> VAL(NumMyNonzeros);
895 Read(GroupName,
"VAL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_DOUBLE, &VAL[0]);
897 Epetra_FECrsMatrix* A2 =
new Epetra_FECrsMatrix(Copy, DomainMap, 0);
899 for (
int i = 0; i < NumMyNonzeros; )
902 while (ROW[i + count] == ROW[i])
905 A2->InsertGlobalValues(1, &ROW[i], count, &COL[i], &VAL[i]);
909 A2->GlobalAssemble(DomainMap, RangeMap);
918 int& NumGlobalNonzeros,
919 int& NumGlobalDiagonals,
926 "requested group " + GroupName +
" not found"));
929 Read(GroupName,
"__type__", Label);
931 if (Label !=
"Epetra_RowMatrix")
933 "requested group " + GroupName +
" is not an Epetra_RowMatrix",
934 "__type__ = " + Label));
936 Read(GroupName,
"NumGlobalRows", NumGlobalRows);
937 Read(GroupName,
"NumGlobalCols", NumGlobalCols);
938 Read(GroupName,
"NumGlobalNonzeros", NumGlobalNonzeros);
939 Read(GroupName,
"NumGlobalDiagonals", NumGlobalDiagonals);
940 Read(GroupName,
"MaxNumEntries", MaxNumEntries);
941 Read(GroupName,
"NormOne", NormOne);
942 Read(GroupName,
"NormInf", NormInf);
953 throw(
Exception(__FILE__, __LINE__,
"no file open yet"));
955 hid_t group_id = H5Gcreate(file_id_, GroupName.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
961 status = H5Gclose(group_id);
964 Write(GroupName,
"__type__",
"Teuchos::ParameterList");
971 throw(
Exception(__FILE__, __LINE__,
"no file open yet"));
974 Read(GroupName,
"__type__", Label);
976 if (Label !=
"Teuchos::ParameterList")
978 "requested group " + GroupName +
" is not a Teuchos::ParameterList",
979 "__type__ = " + Label));
986 group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
990 std::string xxx =
"/" + GroupName;
991 status = H5Giterate(group_id, xxx.c_str() , NULL,
f_operator, ¶ms);
995 status = H5Gclose(group_id);
1008 throw(
Exception(__FILE__, __LINE__,
"no file open yet"));
1010 hid_t group_id, dset_id;
1011 hid_t filespace_id, memspace_id;
1015 Teuchos::RCP<Epetra_MultiVector> LinearX;
1017 if (X.Map().LinearMap())
1018 LinearX = Teuchos::rcp(const_cast<Epetra_MultiVector*>(&X),
false);
1021 Epetra_Map LinearMap(X.GlobalLength(), X.Map().IndexBase(), Comm_);
1022 LinearX = Teuchos::rcp(
new Epetra_MultiVector(LinearMap, X.NumVectors()));
1023 Epetra_Import Importer(LinearMap, X.Map());
1024 LinearX->Import(X, Importer, Insert);
1027 int NumVectors = X.NumVectors();
1028 int GlobalLength = X.GlobalLength();
1034 if (writeTranspose) indexT = 1;
1036 hsize_t q_dimsf[] = {
static_cast<hsize_t
>(GlobalLength), static_cast<hsize_t>(GlobalLength)};
1037 q_dimsf[indexT] = NumVectors;
1039 filespace_id = H5Screate_simple(2, q_dimsf, NULL);
1044 group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1047 dset_id = H5Dcreate(group_id,
"Values", H5T_NATIVE_DOUBLE, filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1050 plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1052 H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1057 hsize_t offset[] = {
static_cast<hsize_t
>(LinearX->Map().GID(0)-X.Map().IndexBase()),
1058 static_cast<hsize_t>(LinearX->Map().GID(0)-X.Map().IndexBase())};
1059 hsize_t stride[] = {1, 1};
1060 hsize_t count[] = {
static_cast<hsize_t
>(LinearX->MyLength()),
1061 static_cast<hsize_t>(LinearX->MyLength())};
1062 hsize_t block[] = {1, 1};
1065 for (
int n(0); n < NumVectors; ++n)
1071 filespace_id = H5Dget_space(dset_id);
1072 H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1076 hsize_t dimsm[] = {
static_cast<hsize_t
>(LinearX->MyLength())};
1077 memspace_id = H5Screate_simple(1, dimsm, NULL);
1080 status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1081 plist_id_, LinearX->operator[](n));
1085 H5Sclose(memspace_id);
1086 H5Sclose(filespace_id);
1088 H5Pclose(plist_id_);
1090 Write(GroupName,
"GlobalLength", GlobalLength);
1091 Write(GroupName,
"NumVectors", NumVectors);
1092 Write(GroupName,
"__type__",
"Epetra_MultiVector");
1097 Epetra_MultiVector*& X,
bool writeTranspose)
1100 Epetra_MultiVector* LinearX;
1101 Read(GroupName, LinearX, writeTranspose, Map.IndexBase());
1104 Epetra_Import Importer(Map, LinearX->Map());
1105 X =
new Epetra_MultiVector(Map, LinearX->NumVectors());
1106 X->Import(*LinearX, Importer, Insert);
1114 bool readTranspose,
const int& indexBase)
1116 int GlobalLength, NumVectors;
1127 if (readTranspose) indexT = 1;
1129 hsize_t q_dimsf[] = {
static_cast<hsize_t
>(GlobalLength), static_cast<hsize_t>(GlobalLength)};
1130 q_dimsf[indexT] = NumVectors;
1132 hid_t filespace_id = H5Screate_simple(2, q_dimsf, NULL);
1136 "requested group " + GroupName +
" not found"));
1138 group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1141 hid_t dset_id = H5Dopen(group_id,
"Values", H5P_DEFAULT);
1144 plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1146 H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1148 H5Pclose(plist_id_);
1150 Epetra_Map LinearMap(GlobalLength, indexBase, Comm());
1151 LinearX =
new Epetra_MultiVector(LinearMap, NumVectors);
1154 hsize_t offset[] = {
static_cast<hsize_t
>(LinearMap.GID(0) - indexBase), static_cast<hsize_t>(LinearMap.GID(0) - indexBase)};
1155 hsize_t stride[] = {1, 1};
1162 hsize_t count[] = {1, 1};
1163 hsize_t block[] = {
static_cast<hsize_t
>(LinearX->MyLength()), static_cast<hsize_t>(LinearX->MyLength())};
1166 count [indexT] = NumVectors;
1169 filespace_id = H5Dget_space(dset_id);
1170 H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1174 hsize_t dimsm[] = {NumVectors *
static_cast<hsize_t
>(LinearX->MyLength())};
1175 memspace_id = H5Screate_simple(1, dimsm, NULL);
1178 CHECK_STATUS(H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1179 H5P_DEFAULT, LinearX->Values()));
1185 hsize_t count[] = {
static_cast<hsize_t
>(LinearX->MyLength()),
1186 static_cast<hsize_t>(LinearX->MyLength())};
1187 hsize_t block[] = {1, 1};
1190 for (
int n(0); n < NumVectors; ++n)
1196 filespace_id = H5Dget_space(dset_id);
1197 H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride,
1201 hsize_t dimsm[] = {
static_cast<hsize_t
>(LinearX->MyLength())};
1202 memspace_id = H5Screate_simple(1, dimsm, NULL);
1205 CHECK_STATUS(H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id,
1206 H5P_DEFAULT, LinearX->operator[](n)));
1224 "requested group " + GroupName +
" not found"));
1227 Read(GroupName,
"__type__", Label);
1229 if (Label !=
"Epetra_MultiVector")
1231 "requested group " + GroupName +
" is not an Epetra_MultiVector",
1232 "__type__ = " + Label));
1234 Read(GroupName,
"GlobalLength", GlobalLength);
1235 Read(GroupName,
"NumVectors", NumVectors);
1245 if (x.Map().LinearMap())
1247 Write(GroupName,
"Values", x.Map().NumMyElements() * x.
RowSize(),
1248 x.Map().NumGlobalElements() * x.
RowSize(),
1249 H5T_NATIVE_INT, x.
Values());
1255 const Epetra_BlockMap& OriginalMap = x.Map();
1256 Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
1257 Epetra_Import Importer(LinearMap, OriginalMap);
1260 LinearX.Import(x, Importer, Insert);
1262 Write(GroupName,
"Values", LinearMap.NumMyElements() * x.
RowSize(),
1263 LinearMap.NumGlobalElements() * x.
RowSize(),
1264 H5T_NATIVE_INT, LinearX.Values());
1267 Write(GroupName,
"__type__",
"EpetraExt::DistArray<int>");
1277 int GlobalLength, RowSize;
1280 if (Map.LinearMap())
1284 Read(GroupName,
"Values", Map.NumMyElements() * RowSize,
1285 Map.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->
Values());
1291 Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
1294 Read(GroupName,
"Values", LinearMap.NumMyElements() * RowSize,
1295 LinearMap.NumGlobalElements() * RowSize,
1296 H5T_NATIVE_INT, LinearX.Values());
1298 Epetra_Import Importer(Map, LinearMap);
1300 X->Import(LinearX, Importer, Insert);
1308 int GlobalLength, RowSize;
1312 Epetra_Map LinearMap(GlobalLength, 0, Comm_);
1315 Read(GroupName,
"Values", LinearMap.NumMyElements() * RowSize,
1316 LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->
Values());
1326 "requested group " + GroupName +
" not found"));
1329 Read(GroupName,
"__type__", Label);
1331 if (Label !=
"EpetraExt::DistArray<int>")
1333 "requested group " + GroupName +
" is not an EpetraExt::DistArray<int>",
1334 "__type__ = " + Label));
1336 Read(GroupName,
"GlobalLength", GlobalLength);
1337 Read(GroupName,
"RowSize", RowSize);
1347 if (x.Map().LinearMap())
1349 Write(GroupName,
"Values", x.Map().NumMyElements() * x.
RowSize(),
1350 x.Map().NumGlobalElements() * x.
RowSize(),
1351 H5T_NATIVE_DOUBLE, x.
Values());
1357 const Epetra_BlockMap& OriginalMap = x.Map();
1358 Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
1359 Epetra_Import Importer(LinearMap, OriginalMap);
1362 LinearX.Import(x, Importer, Insert);
1364 Write(GroupName,
"Values", LinearMap.NumMyElements() * x.
RowSize(),
1365 LinearMap.NumGlobalElements() * x.
RowSize(),
1366 H5T_NATIVE_DOUBLE, LinearX.Values());
1369 Write(GroupName,
"__type__",
"EpetraExt::DistArray<double>");
1379 int GlobalLength, RowSize;
1382 if (Map.LinearMap())
1386 Read(GroupName,
"Values", Map.NumMyElements() * RowSize,
1387 Map.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->
Values());
1393 Epetra_Map LinearMap(GlobalLength, Map.IndexBase(), Comm_);
1396 Read(GroupName,
"Values", LinearMap.NumMyElements() * RowSize,
1397 LinearMap.NumGlobalElements() * RowSize,
1398 H5T_NATIVE_DOUBLE, LinearX.Values());
1400 Epetra_Import Importer(Map, LinearMap);
1402 X->Import(LinearX, Importer, Insert);
1410 int GlobalLength, RowSize;
1414 Epetra_Map LinearMap(GlobalLength, 0, Comm_);
1417 Read(GroupName,
"Values", LinearMap.NumMyElements() * RowSize,
1418 LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->
Values());
1428 "requested group " + GroupName +
" not found"));
1431 Read(GroupName,
"__type__", Label);
1433 if (Label !=
"EpetraExt::DistArray<double>")
1435 "requested group " + GroupName +
" is not an EpetraExt::DistArray<double>",
1436 "__type__ = " + Label));
1438 Read(GroupName,
"GlobalLength", GlobalLength);
1439 Read(GroupName,
"RowSize", RowSize);
1453 hid_t filespace_id = H5Screate(H5S_SCALAR);
1454 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1455 hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_INT,
1456 filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1458 herr_t status = H5Dwrite(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
1459 H5P_DEFAULT, &what);
1465 H5Sclose(filespace_id);
1475 hid_t filespace_id = H5Screate(H5S_SCALAR);
1476 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1477 hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_DOUBLE,
1478 filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1480 herr_t status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL,
1481 filespace_id, H5P_DEFAULT, &what);
1487 H5Sclose(filespace_id);
1495 "requested group " + GroupName +
" not found"));
1498 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1500 hid_t filespace_id = H5Screate(H5S_SCALAR);
1501 hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1503 herr_t status = H5Dread(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
1504 H5P_DEFAULT, &data);
1507 H5Sclose(filespace_id);
1517 "requested group " + GroupName +
" not found"));
1520 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1522 hid_t filespace_id = H5Screate(H5S_SCALAR);
1523 hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1525 herr_t status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, filespace_id,
1526 H5P_DEFAULT, &data);
1529 H5Sclose(filespace_id);
1536 const std::string& DataSetName,
1537 const std::string& data)
1544 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1546 hid_t dataspace_id = H5Screate_simple(1, &len, NULL);
1548 hid_t atype = H5Tcopy(H5T_C_S1);
1549 H5Tset_size(atype, data.size() + 1);
1551 hid_t dataset_id = H5Dcreate(group_id, DataSetName.c_str(), atype,
1552 dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1553 CHECK_STATUS(H5Dwrite(dataset_id, atype, H5S_ALL, H5S_ALL,
1554 H5P_DEFAULT, data.c_str()));
1567 const std::string& DataSetName,
1572 "requested group " + GroupName +
" not found"));
1574 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1576 hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1578 hid_t datatype_id = H5Dget_type(dataset_id);
1580 H5T_class_t typeclass_id = H5Tget_class(datatype_id);
1582 if(typeclass_id != H5T_STRING)
1584 "requested group " + GroupName +
" is not a std::string"));
1586 CHECK_STATUS(H5Dread(dataset_id, datatype_id, H5S_ALL, H5S_ALL,
1587 H5P_DEFAULT, data2));
1590 H5Dclose(dataset_id);
1600 const int type,
const int Length,
1606 hsize_t dimsf = Length;
1608 hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
1610 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1612 hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type,
1613 filespace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1615 herr_t status = H5Dwrite(dset_id, type, H5S_ALL, H5S_ALL,
1622 H5Sclose(filespace_id);
1627 const int type,
const int Length,
1632 "requested group " + GroupName +
" not found"));
1634 hsize_t dimsf = Length;
1636 hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
1638 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1640 hid_t dset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1642 herr_t status = H5Dread(dset_id, type, H5S_ALL, filespace_id,
1649 H5Sclose(filespace_id);
1658 int MySize,
int GlobalSize,
int type,
const void* data)
1661 Comm_.ScanSum(&MySize, &Offset, 1);
1664 hsize_t MySize_t = MySize;
1665 hsize_t GlobalSize_t = GlobalSize;
1666 hsize_t Offset_t = Offset;
1668 hid_t filespace_id = H5Screate_simple(1, &GlobalSize_t, NULL);
1674 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1675 hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type, filespace_id,
1676 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1678 H5Sclose(filespace_id);
1683 hid_t memspace_id = H5Screate_simple(1, &MySize_t, NULL);
1686 filespace_id = H5Dget_space(dset_id);
1687 H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL, &MySize_t, NULL);
1689 plist_id_ = H5Pcreate(H5P_DATASET_XFER);
1691 H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
1694 status = H5Dwrite(dset_id, type, memspace_id, filespace_id,
1701 H5Sclose(filespace_id);
1702 H5Sclose(memspace_id);
1703 H5Pclose(plist_id_);
1708 int MySize,
int GlobalSize,
1709 const int type,
void* data)
1712 throw(
Exception(__FILE__, __LINE__,
"no file open yet"));
1714 hsize_t MySize_t = MySize;
1718 Comm_.ScanSum(&MySize, &itmp, 1);
1719 hsize_t Offset_t = itmp - MySize;
1721 hid_t group_id = H5Gopen(file_id_, GroupName.c_str(), H5P_DEFAULT);
1722 hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str(), H5P_DEFAULT);
1726 hid_t filespace_id = H5Dget_space(dataset_id);
1727 H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL,
1730 hid_t mem_dataspace = H5Screate_simple (1, &MySize_t, NULL);
1732 herr_t status = H5Dread(dataset_id, type, mem_dataspace, filespace_id,
1736 H5Sclose(mem_dataspace);
1739 H5Dclose(dataset_id);
void ReadDoubleDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
void Write(const std::string &GroupName, const std::string &DataSetName, int data)
Write an integer in group GroupName using the given DataSetName.
int GlobalLength() const
Returns the global length of the array.
void ReadMapProperties(const std::string &GroupName, int &NumGlobalElements, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_Map.
void CreateGroup(const std::string &GroupName)
Create group GroupName.
static void WriteParameterListRecursive(const Teuchos::ParameterList ¶ms, hid_t group_id)
bool IsOpen() const
Return true if a file has already been opened using Open()/Create()
void Read(const std::string &GroupName, const std::string &DataSetName, int &data)
Read an integer from group /GroupName/DataSetName.
T * Values()
Returns a pointer to the internally stored data (non-const version).
void ReadBlockMapProperties(const std::string &GroupName, int &NumGlobalElements, int &NumGlobalPoints, int &IndexBase, int &NumProc)
Read basic properties of specified Epetra_BlockMap.
void ReadIntVectorProperties(const std::string &GroupName, int &GlobalLength)
Read basic properties of specified Epetra_IntVector.
void ReadCrsGraphProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumGlobalNonzeros, int &NumGlobalDiagonals, int &MaxNumIndices)
Read basic properties of specified Epetra_CrsGraph.
bool IsContained(const std::string Name)
Return true if Name is contained in the database.
int RowSize() const
Returns the row size, that is, the amount of data associated with each element.
void Open(const std::string FileName, int AccessType=H5F_ACC_RDWR)
Open specified file with given access type.
#define CHECK_STATUS(status)
static herr_t f_operator(hid_t loc_id, const char *name, void *opdata)
HDF5(const Epetra_Comm &Comm)
Constructor.
void Create(const std::string FileName)
Create a new file.
void ReadIntDistArrayProperties(const std::string &GroupName, int &GlobalLength, int &RowSize)
Read the global number of elements and type for a generic handle object.
void ReadCrsMatrixProperties(const std::string &GroupName, int &NumGlobalRows, int &NumGlobalCols, int &NumNonzeros, int &NumGlobalDiagonals, int &MaxNumEntries, double &NormOne, double &NormInf)
Read basic properties of specified Epetra_CrsMatrix.
static herr_t FindDataset(hid_t loc_id, const char *name, void *opdata)
void ReadMultiVectorProperties(const std::string &GroupName, int &GlobalLength, int &NumVectors)
Read basic properties of specified Epetra_MultiVector.
std::string toString(const int &x)
DistArray<T>: A class to store row-oriented multivectors of type T.