EpetraExt_HDF5.cpp

Go to the documentation of this file.
00001 #include "EpetraExt_ConfigDefs.h"
00002 #ifdef HAVE_EPETRAEXT_HDF5
00003 #ifdef HAVE_MPI
00004 #include "Epetra_MpiComm.h"
00005 #include "mpi.h"
00006 #else
00007 #include "Epetra_SerialComm.h"
00008 #endif
00009 #include "Teuchos_ParameterList.hpp"
00010 #include "Teuchos_RefCountPtr.hpp"
00011 #include "Epetra_Map.h"
00012 #include "Epetra_BlockMap.h"
00013 #include "Epetra_CrsGraph.h"
00014 #include "Epetra_FECrsGraph.h"
00015 #include "Epetra_RowMatrix.h"
00016 #include "Epetra_CrsMatrix.h"
00017 #include "Epetra_FECrsMatrix.h"
00018 #include "Epetra_IntVector.h"
00019 #include "Epetra_MultiVector.h"
00020 #include "Epetra_Import.h"
00021 #include "EpetraExt_Exception.h"
00022 #include "EpetraExt_Utils.h"
00023 #include "EpetraExt_HDF5.h"
00024 #include "EpetraExt_DistArray.h"
00025 
00026 #define CHECK_HID(hid_t) \
00027   { if (hid_t < 0) \
00028     throw(EpetraExt::Exception(__FILE__, __LINE__, \
00029                     "hid_t is negative")); }
00030 
00031 #define CHECK_STATUS(status) \
00032   { if (status < 0) \
00033     throw(EpetraExt::Exception(__FILE__, __LINE__, \
00034                     "function H5Giterater returned a negative value")); }
00035 
00036 // ==========================================================================
00037 // data container and iterators to find a dataset with a given name
00038 struct FindDataset_t
00039 {
00040   std::string name;
00041   bool found;
00042 };
00043 
00044 static herr_t FindDataset(hid_t loc_id, const char *name, void *opdata)
00045 {
00046   std::string& token = ((FindDataset_t*)opdata)->name;
00047   if (token == name)
00048     ((FindDataset_t*)opdata)->found = true;
00049 
00050   return(0);
00051 }
00052 
00053 // ==========================================================================
00054 // This function copied from Roman Geus' FEMAXX code
00055 static void WriteParameterListRecursive(const Teuchos::ParameterList& params, 
00056                                         hid_t group_id) 
00057 {
00058   Teuchos::ParameterList::ConstIterator it = params.begin();
00059   for (; it != params.end(); ++ it) 
00060   {
00061     std::string key(params.name(it));
00062     if (params.isSublist(key)) 
00063     {
00064       // Sublist
00065 
00066       // Create subgroup for sublist.
00067       hid_t child_group_id = H5Gcreate(group_id, key.c_str(), 0);
00068       WriteParameterListRecursive(params.sublist(key), child_group_id);
00069       H5Gclose(child_group_id);
00070     } 
00071     else 
00072     {
00073       //
00074       // Regular parameter
00075       //
00076 
00077       // Create dataspace/dataset.
00078       herr_t status;
00079       hsize_t one = 1;
00080       hid_t dataspace_id, dataset_id;
00081       bool found = false; // to avoid a compiler error on MAC OS X GCC 4.0.0
00082 
00083       // Write the dataset.
00084       if (params.isType<std::string>(key)) 
00085       {
00086         std::string value = params.get<std::string>(key);
00087         hsize_t len = value.size() + 1;
00088         dataspace_id = H5Screate_simple(1, &len, NULL);
00089         dataset_id = H5Dcreate(group_id, key.c_str(), H5T_C_S1, dataspace_id, H5P_DEFAULT);
00090         status = H5Dwrite(dataset_id, H5T_C_S1, H5S_ALL, H5S_ALL, 
00091                           H5P_DEFAULT, value.c_str());
00092         CHECK_STATUS(status);
00093         status = H5Dclose(dataset_id);
00094         CHECK_STATUS(status);
00095         status = H5Sclose(dataspace_id);
00096         CHECK_STATUS(status);
00097         found = true;
00098       } 
00099 
00100       if (params.isType<bool>(key)) 
00101       {
00102         // Use H5T_NATIVE_USHORT to store a bool value
00103         unsigned short value = params.get<bool>(key) ? 1 : 0;
00104         dataspace_id = H5Screate_simple(1, &one, NULL);
00105         dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_USHORT, dataspace_id, H5P_DEFAULT);
00106         status = H5Dwrite(dataset_id, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL, 
00107                           H5P_DEFAULT, &value);
00108         CHECK_STATUS(status);
00109         status = H5Dclose(dataset_id);
00110         CHECK_STATUS(status);
00111         status = H5Sclose(dataspace_id);
00112         CHECK_STATUS(status);
00113         found = true;
00114       } 
00115       
00116       if (params.isType<int>(key)) 
00117       {
00118         int value = params.get<int>(key);
00119         dataspace_id = H5Screate_simple(1, &one, NULL);
00120         dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_INT, dataspace_id, H5P_DEFAULT);
00121         status = H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, 
00122                           H5P_DEFAULT, &value);
00123         CHECK_STATUS(status);
00124         status = H5Dclose(dataset_id);
00125         CHECK_STATUS(status);
00126         status = H5Sclose(dataspace_id);
00127         CHECK_STATUS(status);
00128         found = true;
00129       } 
00130 
00131       if (params.isType<double>(key)) 
00132       {
00133         double value = params.get<double>(key);
00134         dataspace_id = H5Screate_simple(1, &one, NULL);
00135         dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT);
00136         status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, 
00137                           H5P_DEFAULT, &value);
00138         CHECK_STATUS(status);
00139         status = H5Dclose(dataset_id);
00140         CHECK_STATUS(status);
00141         status = H5Sclose(dataspace_id);
00142         CHECK_STATUS(status);
00143         found = true;
00144       } 
00145 
00146       if (!found)
00147       {
00148         throw(EpetraExt::Exception(__FILE__, __LINE__, 
00149                                 "type for parameter " + key + " not supported"));
00150       }
00151     }
00152   }
00153 }
00154 
00155 // ==========================================================================
00156 // Recursive Operator function called by H5Giterate for each entity in group.
00157 // This function copied from Roman Geus' FEMAXX code
00158 static herr_t f_operator(hid_t loc_id, const char *name, void *opdata) 
00159 {
00160   H5G_stat_t statbuf;
00161   hid_t dataset_id, space_id, type_id;
00162   Teuchos::ParameterList* sublist;
00163   Teuchos::ParameterList* params = (Teuchos::ParameterList*)opdata;
00164   /*
00165    * Get type of the object and display its name and type.
00166    * The name of the object is passed to this function by 
00167    * the Library. Some magic :-)
00168    */
00169   H5Gget_objinfo(loc_id, name, 0, &statbuf);
00170   if (strcmp(name, "__type__") == 0)
00171     return(0); // skip list identifier
00172 
00173   switch (statbuf.type) {
00174   case H5G_GROUP:
00175     sublist = &params->sublist(name);
00176     H5Giterate(loc_id, name , NULL, f_operator, sublist);
00177     break;
00178   case H5G_DATASET:
00179     hsize_t len;
00180     dataset_id = H5Dopen(loc_id, name);
00181     space_id = H5Dget_space(dataset_id);
00182     if (H5Sget_simple_extent_ndims(space_id) != 1)
00183       throw(EpetraExt::Exception(__FILE__, __LINE__, 
00184                               "dimensionality of parameters must be 1."));
00185     H5Sget_simple_extent_dims(space_id, &len, NULL);
00186     type_id = H5Dget_type(dataset_id);
00187     if (H5Tequal(type_id, H5T_NATIVE_DOUBLE) > 0) {
00188       double value;
00189       H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
00190       params->set(name, value);
00191     } else if (H5Tequal(type_id, H5T_NATIVE_INT) > 0) {
00192       int value;
00193       H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
00194       params->set(name, value);
00195     } else if (H5Tequal(type_id, H5T_C_S1) > 0) {
00196       char* buf = new char[len];
00197       H5Dread(dataset_id, H5T_C_S1, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
00198       params->set(name, std::string(buf));
00199       delete[] buf;
00200     } else if (H5Tequal(type_id, H5T_NATIVE_USHORT) > 0) {
00201       unsigned short value;
00202       H5Dread(dataset_id, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
00203       params->set(name, value != 0 ? true : false);
00204     } else {
00205       throw(EpetraExt::Exception(__FILE__, __LINE__,
00206                               "unsupported datatype")); // FIXME
00207     }
00208     H5Tclose(type_id);
00209     H5Sclose(space_id);  
00210     H5Dclose(dataset_id);  
00211     break;
00212   default:
00213     throw(EpetraExt::Exception(__FILE__, __LINE__,
00214                             "unsupported datatype")); // FIXME
00215   }
00216   return 0;
00217 }
00218 
00219 // ==========================================================================
00220 EpetraExt::HDF5::HDF5(const Epetra_Comm& Comm) :
00221   Comm_(Comm),
00222   IsOpen_(false)
00223 {}
00224 
00225 // ==========================================================================
00226 void EpetraExt::HDF5::Create(const std::string FileName)
00227 {
00228   if (IsOpen())
00229     throw(Exception(__FILE__, __LINE__,
00230                     "an HDF5 is already open, first close the current one",
00231                     "using method Close(), then open/create a new one"));
00232 
00233   FileName_ = FileName;
00234 
00235   // Set up file access property list with parallel I/O access
00236   plist_id_ = H5Pcreate(H5P_FILE_ACCESS);
00237 #ifdef HAVE_MPI
00238   // Create property list for collective dataset write.
00239   H5Pset_fapl_mpio(plist_id_, MPI_COMM_WORLD, MPI_INFO_NULL);
00240 #endif
00241 
00242 #if 0
00243   unsigned int boh = H5Z_FILTER_MAX;
00244   H5Pset_filter(plist_id_, H5Z_FILTER_DEFLATE, H5Z_FILTER_MAX, 0, &boh);
00245 #endif
00246 
00247   // create the file collectively and release property list identifier.
00248   file_id_ = H5Fcreate(FileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, 
00249                       plist_id_);
00250   H5Pclose(plist_id_);
00251 
00252   IsOpen_ = true;
00253 }
00254 
00255 // ==========================================================================
00256 void EpetraExt::HDF5::Open(const std::string FileName, int AccessType)
00257 {
00258   if (IsOpen())
00259     throw(Exception(__FILE__, __LINE__,
00260                     "an HDF5 is already open, first close the current one",
00261                     "using method Close(), then open/create a new one"));
00262 
00263   FileName_ = FileName;
00264 
00265   // create the file collectively and release property list identifier.
00266   file_id_ = H5Fopen(FileName.c_str(), AccessType, H5P_DEFAULT);
00267 
00268 #ifdef HAVE_MPI
00269 // FIXME: DO I NEED THE MPIO_COLLECTIVE??
00270 //  plist_id_ = H5Pcreate(H5P_DATASET_XFER);
00271 //  H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
00272 #endif
00273 
00274   IsOpen_ = true;
00275 }
00276 
00277 // ==========================================================================
00278 bool EpetraExt::HDF5::IsContained(std::string Name)
00279 {
00280   if (!IsOpen())
00281     throw(Exception(__FILE__, __LINE__, "no file open yet"));
00282 
00283   FindDataset_t data;
00284   data.name = Name;
00285   data.found = false;
00286 
00287   int idx_f = H5Giterate(file_id_, "/", NULL, FindDataset, (void*)&data);
00288 
00289   return(data.found);
00290 }
00291 
00292 // ============================ //
00293 // Epetra_BlockMap / Epetra_Map //
00294 // ============================ //
00295 
00296 // ==========================================================================
00297 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_BlockMap& BlockMap)
00298 {
00299   int NumMyPoints       = BlockMap.NumMyPoints();
00300   int NumMyElements     = BlockMap.NumMyElements();
00301   int NumGlobalElements = BlockMap.NumGlobalElements();
00302   int NumGlobalPoints   = BlockMap.NumGlobalPoints();
00303   int* MyGlobalElements = BlockMap.MyGlobalElements();
00304   int* ElementSizeList  = BlockMap.ElementSizeList();
00305 
00306   std::vector<int> NumMyElements_v(Comm_.NumProc());
00307   Comm_.GatherAll(&NumMyElements, &NumMyElements_v[0], 1);
00308 
00309   std::vector<int> NumMyPoints_v(Comm_.NumProc());
00310   Comm_.GatherAll(&NumMyPoints, &NumMyPoints_v[0], 1);
00311 
00312   Write(Name, "MyGlobalElements", NumMyElements, NumGlobalElements, 
00313         H5T_NATIVE_INT, MyGlobalElements);
00314   Write(Name, "ElementSizeList", NumMyElements, NumGlobalElements, 
00315         H5T_NATIVE_INT, ElementSizeList);
00316   Write(Name, "NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
00317 
00318   // need to know how many processors currently host this map
00319   Write(Name, "NumProc", Comm_.NumProc());
00320   // write few more data about this map
00321   Write(Name, "NumGlobalPoints", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalPoints);
00322   Write(Name, "NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &NumGlobalElements);
00323   Write(Name, "IndexBase", BlockMap.IndexBase());
00324   Write(Name, "__type__", "Epetra_BlockMap");
00325 }
00326 
00327 // ==========================================================================
00328 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_BlockMap*& BlockMap)
00329 {
00330   int NumGlobalElements, NumGlobalPoints, IndexBase, NumProc;
00331 
00332   ReadBlockMapProperties(GroupName, NumGlobalElements, NumGlobalPoints,
00333                         IndexBase, NumProc);
00334 
00335   std::vector<int> NumMyPoints_v(Comm_.NumProc());
00336   std::vector<int> NumMyElements_v(Comm_.NumProc());
00337 
00338   Read(GroupName, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
00339   Read(GroupName, "NumMyPoints", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyPoints_v[0]);
00340   int NumMyElements = NumMyElements_v[Comm_.MyPID()];
00341   int NumMyPoints   = NumMyPoints_v[Comm_.MyPID()];
00342 
00343   if (NumProc != Comm_.NumProc())
00344     throw(Exception(__FILE__, __LINE__,
00345                     "requested map not compatible with current number of",
00346                     "processors, " + toString(Comm_.NumProc()) + 
00347                     " vs. " + toString(NumProc)));
00348 
00349   std::vector<int> MyGlobalElements(NumMyElements);
00350   std::vector<int> ElementSizeList(NumMyElements);
00351 
00352   Read(GroupName, "MyGlobalElements", NumMyElements, NumGlobalElements, 
00353        H5T_NATIVE_INT, &MyGlobalElements[0]);
00354 
00355   Read(GroupName, "ElementSizeList", NumMyElements, NumGlobalElements, 
00356        H5T_NATIVE_INT, &ElementSizeList[0]);
00357 
00358   BlockMap = new Epetra_BlockMap(NumGlobalElements, NumMyElements, 
00359                                  &MyGlobalElements[0], &ElementSizeList[0],
00360                                  IndexBase, Comm_);
00361 }
00362 
00363 // ==========================================================================
00364 void EpetraExt::HDF5::ReadBlockMapProperties(const std::string& GroupName, 
00365                                           int& NumGlobalElements,
00366                                           int& NumGlobalPoints,
00367                                           int& IndexBase,
00368                                           int& NumProc)
00369 {
00370   if (!IsContained(GroupName))
00371     throw(Exception(__FILE__, __LINE__,
00372                     "requested group " + GroupName + " not found"));
00373 
00374   std::string Label;
00375   Read(GroupName, "__type__", Label);
00376 
00377   if (Label != "Epetra_BlockMap")
00378     throw(Exception(__FILE__, __LINE__,
00379                     "requested group " + GroupName + " is not an Epetra_BlockMap",
00380                     "__type__ = " + Label));
00381 
00382   Read(GroupName, "NumGlobalElements", NumGlobalElements);
00383   Read(GroupName, "NumGlobalPoints", NumGlobalPoints);
00384   Read(GroupName, "IndexBase", IndexBase);
00385   Read(GroupName, "NumProc", NumProc);
00386 }
00387 
00388 // ==========================================================================
00389 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_Map& Map)
00390 {
00391   int MySize = Map.NumMyElements();
00392   int GlobalSize = Map.NumGlobalElements();
00393   int* MyGlobalElements = Map.MyGlobalElements();
00394 
00395   std::vector<int> NumMyElements(Comm_.NumProc());
00396   Comm_.GatherAll(&MySize, &NumMyElements[0], 1);
00397 
00398   Write(Name, "MyGlobalElements", MySize, GlobalSize, 
00399         H5T_NATIVE_INT, MyGlobalElements);
00400   Write(Name, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements[0]);
00401   Write(Name, "NumGlobalElements", 1, Comm_.NumProc(), H5T_NATIVE_INT, &GlobalSize);
00402   Write(Name, "NumProc", Comm_.NumProc());
00403   Write(Name, "IndexBase", Map.IndexBase());
00404   Write(Name, "__type__", "Epetra_Map");
00405 }
00406 
00407 // ==========================================================================
00408 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_Map*& Map)
00409 {
00410   int NumGlobalElements, IndexBase, NumProc;
00411 
00412   ReadMapProperties(GroupName, NumGlobalElements, IndexBase, NumProc);
00413 
00414   std::vector<int> NumMyElements_v(Comm_.NumProc());
00415 
00416   Read(GroupName, "NumMyElements", H5T_NATIVE_INT, Comm_.NumProc(), &NumMyElements_v[0]);
00417   int NumMyElements = NumMyElements_v[Comm_.MyPID()];
00418 
00419   if (NumProc != Comm_.NumProc())
00420     throw(Exception(__FILE__, __LINE__,
00421                     "requested map not compatible with current number of",
00422                     "processors, " + toString(Comm_.NumProc()) + 
00423                     " vs. " + toString(NumProc)));
00424 
00425   std::vector<int> MyGlobalElements(NumMyElements);
00426 
00427   Read(GroupName, "MyGlobalElements", NumMyElements, NumGlobalElements, 
00428        H5T_NATIVE_INT, &MyGlobalElements[0]);
00429 
00430   Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], IndexBase, Comm_);
00431 }
00432 
00433 // ==========================================================================
00434 void EpetraExt::HDF5::ReadMapProperties(const std::string& GroupName, 
00435                                      int& NumGlobalElements,
00436                                      int& IndexBase,
00437                                      int& NumProc)
00438 {
00439   if (!IsContained(GroupName))
00440     throw(Exception(__FILE__, __LINE__,
00441                     "requested group " + GroupName + " not found"));
00442 
00443   std::string Label;
00444   Read(GroupName, "__type__", Label);
00445 
00446   if (Label != "Epetra_Map")
00447     throw(Exception(__FILE__, __LINE__,
00448                     "requested group " + GroupName + " is not an Epetra_Map",
00449                     "__type__ = " + Label));
00450 
00451   Read(GroupName, "NumGlobalElements", NumGlobalElements);
00452   Read(GroupName, "IndexBase", IndexBase);
00453   Read(GroupName, "NumProc", NumProc);
00454 }
00455 
00456 // ================ //
00457 // Epetra_IntVector //
00458 // ================ //
00459 
00460 // ==========================================================================
00461 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_IntVector& x)
00462 {
00463   if (x.Map().LinearMap())
00464   {
00465     Write(Name, "GlobalLength", x.GlobalLength());
00466     Write(Name, "Values", x.Map().NumMyElements(), x.Map().NumGlobalElements(),
00467           H5T_NATIVE_INT, x.Values());
00468   }
00469   else
00470   {
00471     // need to build a linear map first, the import data, then
00472     // finally write them 
00473     const Epetra_BlockMap& OriginalMap = x.Map();
00474     Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
00475     Epetra_Import Importer(LinearMap, OriginalMap);
00476 
00477     Epetra_IntVector LinearX(LinearMap);
00478     LinearX.Import(x, Importer, Insert);
00479 
00480     Write(Name, "GlobalLength", x.GlobalLength());
00481     Write(Name, "Values", LinearMap.NumMyElements(), LinearMap.NumGlobalElements(),
00482           H5T_NATIVE_INT, LinearX.Values());
00483   }
00484   Write(Name, "__type__", "Epetra_IntVector");
00485 }
00486 
00487 // ==========================================================================
00488 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_IntVector*& X)
00489 {
00490   // gets the length of the std::vector
00491   int GlobalLength;
00492   ReadIntVectorProperties(GroupName, GlobalLength);
00493 
00494   // creates a new linear map and use it to read data
00495   Epetra_Map LinearMap(GlobalLength, 0, Comm_);
00496   X = new Epetra_IntVector(LinearMap);
00497 
00498   Read(GroupName, "Values", LinearMap.NumMyElements(), 
00499        LinearMap.NumGlobalElements(), H5T_NATIVE_INT, X->Values());
00500 }
00501 
00502 // ==========================================================================
00503 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
00504                         Epetra_IntVector*& X)
00505 {
00506   // gets the length of the std::vector
00507   int GlobalLength;
00508   ReadIntVectorProperties(GroupName, GlobalLength);
00509 
00510   if (Map.LinearMap())
00511   {
00512     X = new Epetra_IntVector(Map);
00513     // simply read stuff and go home
00514     Read(GroupName, "Values", Map.NumMyElements(), Map.NumGlobalElements(),
00515          H5T_NATIVE_INT, X->Values());
00516 
00517   }
00518   else
00519   {
00520     // we need to first create a linear map, read the std::vector,
00521     // then import it to the actual nonlinear map
00522     Epetra_Map LinearMap(GlobalLength, 0, Comm_);
00523     Epetra_IntVector LinearX(LinearMap);
00524 
00525     Read(GroupName, "Values", LinearMap.NumMyElements(), 
00526          LinearMap.NumGlobalElements(),
00527          H5T_NATIVE_INT, LinearX.Values());
00528 
00529     Epetra_Import Importer(Map, LinearMap);
00530     X = new Epetra_IntVector(Map);
00531     X->Import(LinearX, Importer, Insert);
00532   }
00533 }
00534 
00535 // ==========================================================================
00536 void EpetraExt::HDF5::ReadIntVectorProperties(const std::string& GroupName, 
00537                                            int& GlobalLength)
00538 {
00539   if (!IsContained(GroupName))
00540     throw(Exception(__FILE__, __LINE__,
00541                     "requested group " + GroupName + " not found"));
00542 
00543   std::string Label;
00544   Read(GroupName, "__type__", Label);
00545 
00546   if (Label != "Epetra_IntVector")
00547     throw(Exception(__FILE__, __LINE__,
00548                     "requested group " + GroupName + " is not an Epetra_IntVector",
00549                     "__type__ = " + Label));
00550 
00551   Read(GroupName, "GlobalLength", GlobalLength);
00552 }
00553 
00554 // =============== //
00555 // Epetra_CrsGraph //
00556 // =============== //
00557 
00558 // ==========================================================================
00559 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_CrsGraph& Graph)
00560 {
00561   if (!Graph.Filled())
00562     throw(Exception(__FILE__, __LINE__,
00563                     "input Epetra_CrsGraph is not FillComplete()'d"));
00564 
00565   // like RowMatrix, only without values
00566   int MySize = Graph.NumMyNonzeros();
00567   int GlobalSize = Graph.NumGlobalNonzeros();
00568   std::vector<int> ROW(MySize);
00569   std::vector<int> COL(MySize);
00570 
00571   int count = 0;
00572   int* RowIndices;
00573   int NumEntries;
00574 
00575   for (int i = 0; i < Graph.NumMyRows(); ++i)
00576   {
00577     Graph.ExtractMyRowView(i, NumEntries, RowIndices);
00578     for (int j = 0; j < NumEntries; ++j)
00579     {
00580       ROW[count] = Graph.GRID(i);
00581       COL[count] = Graph.GCID(RowIndices[j]);
00582       ++count;
00583     }
00584   }
00585 
00586   Write(Name, "ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
00587   Write(Name, "COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
00588   Write(Name, "MaxNumIndices", Graph.MaxNumIndices());
00589   Write(Name, "NumGlobalRows", Graph.NumGlobalRows());
00590   Write(Name, "NumGlobalCols", Graph.NumGlobalCols());
00591   Write(Name, "NumGlobalNonzeros", Graph.NumGlobalNonzeros());
00592   Write(Name, "NumGlobalDiagonals", Graph.NumGlobalDiagonals());
00593   Write(Name, "__type__", "Epetra_CrsGraph");
00594 }
00595 
00596 // ==========================================================================
00597 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_CrsGraph*& Graph)
00598 {
00599   int NumGlobalRows, NumGlobalCols;
00600   int NumGlobalNonzeros, NumGlobalDiagonals, MaxNumIndices;
00601 
00602   ReadCrsGraphProperties(GroupName, NumGlobalRows,
00603                          NumGlobalCols, NumGlobalNonzeros,
00604                          NumGlobalDiagonals, MaxNumIndices);
00605 
00606   Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
00607   Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
00608 
00609   Read(GroupName, DomainMap, RangeMap, Graph);
00610 }
00611 
00612 // ==========================================================================
00613 void EpetraExt::HDF5::ReadCrsGraphProperties(const std::string& GroupName, 
00614                                           int& NumGlobalRows,
00615                                           int& NumGlobalCols,
00616                                           int& NumGlobalNonzeros,
00617                                           int& NumGlobalDiagonals,
00618                                           int& MaxNumIndices)
00619 {
00620   if (!IsContained(GroupName))
00621     throw(Exception(__FILE__, __LINE__,
00622                     "requested group " + GroupName + " not found"));
00623 
00624   std::string Label;
00625   Read(GroupName, "__type__", Label);
00626 
00627   if (Label != "Epetra_CrsGraph")
00628     throw(Exception(__FILE__, __LINE__,
00629                     "requested group " + GroupName + " is not an Epetra_CrsGraph",
00630                     "__type__ = " + Label));
00631 
00632   Read(GroupName, "NumGlobalRows",      NumGlobalRows);
00633   Read(GroupName, "NumGlobalCols",      NumGlobalCols);
00634   Read(GroupName, "NumGlobalNonzeros",  NumGlobalNonzeros);
00635   Read(GroupName, "NumGlobalDiagonals", NumGlobalDiagonals);
00636   Read(GroupName, "MaxNumIndices",      MaxNumIndices);
00637 }
00638 
00639 // ==========================================================================
00640 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& DomainMap, 
00641                         const Epetra_Map& RangeMap, Epetra_CrsGraph*& Graph)
00642 {
00643   if (!IsContained(GroupName))
00644     throw(Exception(__FILE__, __LINE__,
00645                     "requested group " + GroupName + " not found in database"));
00646 
00647   std::string Label;
00648   Read(GroupName, "__type__", Label);
00649 
00650   if (Label != "Epetra_CrsGraph")
00651     throw(Exception(__FILE__, __LINE__,
00652                     "requested group " + GroupName + " is not an Epetra_CrsGraph",
00653                     "__type__ = " + Label));
00654 
00655   int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
00656   Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
00657   Read(GroupName, "NumGlobalRows", NumGlobalRows);
00658   Read(GroupName, "NumGlobalCols", NumGlobalCols);
00659 
00660   // linear distribution for nonzeros
00661   int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
00662   if (Comm_.MyPID() == 0)
00663     NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
00664 
00665   std::vector<int> ROW(NumMyNonzeros);
00666   Read(GroupName, "ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
00667 
00668   std::vector<int> COL(NumMyNonzeros);
00669   Read(GroupName, "COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
00670 
00671   Epetra_FECrsGraph* Graph2 = new Epetra_FECrsGraph(Copy, DomainMap, 0);
00672 
00673   for (int i = 0; i < NumMyNonzeros; )
00674   {
00675     int count = 1;
00676     while (ROW[i + count] == ROW[i])
00677       ++count;
00678 
00679     Graph2->InsertGlobalIndices(1, &ROW[i], count, &COL[i]);
00680     i += count;
00681   }
00682 
00683   Graph2->FillComplete(DomainMap, RangeMap);
00684 
00685   Graph = Graph2;
00686 }
00687 
00688 // ================ //
00689 // Epetra_RowMatrix //
00690 // ================ //
00691 
00692 // ==========================================================================
00693 void EpetraExt::HDF5::Write(const std::string& Name, const Epetra_RowMatrix& Matrix)
00694 {
00695   if (!Matrix.Filled())
00696     throw(Exception(__FILE__, __LINE__,
00697                     "input Epetra_RowMatrix is not FillComplete()'d"));
00698 
00699   int MySize = Matrix.NumMyNonzeros();
00700   int GlobalSize = Matrix.NumGlobalNonzeros();
00701   std::vector<int> ROW(MySize);
00702   std::vector<int> COL(MySize);
00703   std::vector<double> VAL(MySize);
00704 
00705   int count = 0;
00706   int Length = Matrix.MaxNumEntries();
00707   std::vector<int> Indices(Length);
00708   std::vector<double> Values(Length);
00709   int NumEntries;
00710 
00711   for (int i = 0; i < Matrix.NumMyRows(); ++i)
00712   {
00713     Matrix.ExtractMyRowCopy(i, Length, NumEntries, &Values[0], &Indices[0]);
00714     for (int j = 0; j < NumEntries; ++j)
00715     {
00716       ROW[count] = Matrix.RowMatrixRowMap().GID(i);
00717       COL[count] = Matrix.RowMatrixColMap().GID(Indices[j]);
00718       VAL[count] = Values[j];
00719       ++count;
00720     }
00721   }
00722 
00723   Write(Name, "ROW", MySize, GlobalSize, H5T_NATIVE_INT, &ROW[0]);
00724   Write(Name, "COL", MySize, GlobalSize, H5T_NATIVE_INT, &COL[0]);
00725   Write(Name, "VAL", MySize, GlobalSize, H5T_NATIVE_DOUBLE, &VAL[0]);
00726   Write(Name, "NumGlobalRows", Matrix.NumGlobalRows());
00727   Write(Name, "NumGlobalCols", Matrix.NumGlobalCols());
00728   Write(Name, "NumGlobalNonzeros", Matrix.NumGlobalNonzeros());
00729   Write(Name, "NumGlobalDiagonals", Matrix.NumGlobalDiagonals());
00730   Write(Name, "MaxNumEntries", Matrix.MaxNumEntries());
00731   Write(Name, "NormOne", Matrix.NormOne());
00732   Write(Name, "NormInf", Matrix.NormInf());
00733   Write(Name, "__type__", "Epetra_RowMatrix");
00734 }
00735 
00736 // ==========================================================================
00737 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_CrsMatrix*& A)
00738 {
00739   int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
00740   int NumGlobalDiagonals, MaxNumEntries;
00741   double NormOne, NormInf;
00742 
00743   ReadCrsMatrixProperties(GroupName, NumGlobalRows, NumGlobalCols,
00744                           NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
00745                           NormOne, NormInf);
00746 
00747   // build simple linear maps for domain and range space
00748   Epetra_Map RangeMap(NumGlobalRows, 0, Comm_);
00749   Epetra_Map DomainMap(NumGlobalCols, 0, Comm_);
00750 
00751   Read(GroupName, DomainMap, RangeMap, A);
00752 }
00753 
00754 // ==========================================================================
00755 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& DomainMap, 
00756                         const Epetra_Map& RangeMap, Epetra_CrsMatrix*& A)
00757 {
00758   int NumGlobalRows, NumGlobalCols, NumGlobalNonzeros;
00759   int NumGlobalDiagonals, MaxNumEntries;
00760   double NormOne, NormInf;
00761 
00762   ReadCrsMatrixProperties(GroupName, NumGlobalRows, NumGlobalCols,
00763                           NumGlobalNonzeros, NumGlobalDiagonals, MaxNumEntries,
00764                           NormOne, NormInf);
00765 
00766   int NumMyNonzeros = NumGlobalNonzeros / Comm_.NumProc();
00767   if (Comm_.MyPID() == 0)
00768     NumMyNonzeros += NumGlobalNonzeros % Comm_.NumProc();
00769 
00770   std::vector<int> ROW(NumMyNonzeros);
00771   Read(GroupName, "ROW", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &ROW[0]);
00772 
00773   std::vector<int> COL(NumMyNonzeros);
00774   Read(GroupName, "COL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_INT, &COL[0]);
00775 
00776   std::vector<double> VAL(NumMyNonzeros);
00777   Read(GroupName, "VAL", NumMyNonzeros, NumGlobalNonzeros, H5T_NATIVE_DOUBLE, &VAL[0]);
00778 
00779   Epetra_FECrsMatrix* A2 = new Epetra_FECrsMatrix(Copy, DomainMap, 0);
00780 
00781   for (int i = 0; i < NumMyNonzeros; )
00782   {
00783     int count = 1;
00784     while (ROW[i + count] == ROW[i])
00785       ++count;
00786 
00787     A2->InsertGlobalValues(1, &ROW[i], count, &COL[i], &VAL[i]);
00788     i += count;
00789   }
00790 
00791   A2->FillComplete(DomainMap, RangeMap);
00792 
00793   A = A2;
00794 }
00795 
00796 // ==========================================================================
00797 void EpetraExt::HDF5::ReadCrsMatrixProperties(const std::string& GroupName, 
00798                                            int& NumGlobalRows,
00799                                            int& NumGlobalCols,
00800                                            int& NumGlobalNonzeros,
00801                                            int& NumGlobalDiagonals,
00802                                            int& MaxNumEntries,
00803                                            double& NormOne,
00804                                            double& NormInf)
00805 {
00806   if (!IsContained(GroupName))
00807     throw(Exception(__FILE__, __LINE__,
00808                     "requested group " + GroupName + " not found"));
00809 
00810   std::string Label;
00811   Read(GroupName, "__type__", Label);
00812 
00813   if (Label != "Epetra_RowMatrix")
00814     throw(Exception(__FILE__, __LINE__,
00815                     "requested group " + GroupName + " is not an Epetra_RowMatrix",
00816                     "__type__ = " + Label));
00817 
00818   Read(GroupName, "NumGlobalRows", NumGlobalRows);
00819   Read(GroupName, "NumGlobalCols", NumGlobalCols);
00820   Read(GroupName, "NumGlobalNonzeros", NumGlobalNonzeros);
00821   Read(GroupName, "NumGlobalDiagonals", NumGlobalDiagonals);
00822   Read(GroupName, "MaxNumEntries", MaxNumEntries);
00823   Read(GroupName, "NormOne", NormOne);
00824   Read(GroupName, "NormInf", NormInf);
00825 }
00826 
00827 // ============= //
00828 // ParameterList //
00829 // ============= //
00830 
00831 // ==========================================================================
00832 void EpetraExt::HDF5::Write(const std::string& GroupName, const Teuchos::ParameterList& params)
00833 {
00834   if (!IsOpen())
00835     throw(Exception(__FILE__, __LINE__, "no file open yet"));
00836 
00837   hid_t group_id = H5Gcreate(file_id_, GroupName.c_str(), 0);
00838 
00839   // Iterate through parameter list 
00840   WriteParameterListRecursive(params, group_id);
00841 
00842   // Finalize hdf5 file
00843   status = H5Gclose(group_id);
00844   CHECK_STATUS(status);
00845 
00846   Write(GroupName, "__type__", "Teuchos::ParameterList");
00847 }
00848 
00849 // ==========================================================================
00850 void EpetraExt::HDF5::Read(const std::string& GroupName, Teuchos::ParameterList& params) 
00851 {
00852   if (!IsOpen())
00853     throw(Exception(__FILE__, __LINE__, "no file open yet"));
00854 
00855   std::string Label;
00856   Read(GroupName, "__type__", Label);
00857 
00858   if (Label != "Teuchos::ParameterList")
00859     throw(Exception(__FILE__, __LINE__,
00860                     "requested group " + GroupName + " is not a Teuchos::ParameterList",
00861                     "__type__ = " + Label));
00862 
00863   // Open hdf5 file
00864   hid_t       group_id;  /* identifiers */
00865   herr_t      status;
00866 
00867   // open group in the root group using absolute name.
00868   group_id = H5Gopen(file_id_, GroupName.c_str());
00869   CHECK_HID(group_id);
00870 
00871   // Iterate through parameter list 
00872   std::string xxx = "/" + GroupName;
00873   status = H5Giterate(group_id, xxx.c_str() , NULL, f_operator, &params);
00874   CHECK_STATUS(status);
00875 
00876   // Finalize hdf5 file
00877   status = H5Gclose(group_id);
00878   CHECK_STATUS(status);
00879 }
00880 
00881 // ================== //
00882 // Epetra_MultiVector //
00883 // ================== //
00884 
00885 // ==========================================================================
00886 void EpetraExt::HDF5::Write(const std::string& GroupName, const Epetra_MultiVector& X)
00887 {
00888   if (!IsOpen())
00889     throw(Exception(__FILE__, __LINE__, "no file open yet"));
00890 
00891   hid_t       group_id, dset_id;
00892   hid_t       filespace_id, memspace_id;
00893   herr_t      status;
00894 
00895   // need a linear distribution to use hyperslabs
00896   Teuchos::RefCountPtr<Epetra_MultiVector> LinearX;
00897 
00898   if (Comm().NumProc() == 1 || X.Map().LinearMap())
00899     LinearX = Teuchos::rcp(const_cast<Epetra_MultiVector*>(&X), false);
00900   else
00901   {
00902     Epetra_Map LinearMap(X.GlobalLength(), 0, Comm_);
00903     LinearX = Teuchos::rcp(new Epetra_MultiVector(LinearMap, X.NumVectors()));
00904     Epetra_Import Importer(LinearMap, X.Map());
00905     LinearX->Import(X, Importer, Insert);
00906   }
00907 
00908   int NumVectors = X.NumVectors();
00909   int GlobalLength = X.GlobalLength();
00910 
00911   // Create the dataspace for the dataset.
00912   hsize_t q_dimsf[] = {NumVectors, GlobalLength};
00913   filespace_id = H5Screate_simple(2, q_dimsf, NULL);
00914 
00915   if (!IsContained(GroupName))
00916     CreateGroup(GroupName);
00917 
00918   group_id = H5Gopen(file_id_, GroupName.c_str());
00919 
00920   // Create the dataset with default properties and close filespace_id.
00921   dset_id = H5Dcreate(group_id, "Values", H5T_NATIVE_DOUBLE, filespace_id, H5P_DEFAULT);
00922 
00923   // Create property list for collective dataset write.
00924   plist_id_ = H5Pcreate(H5P_DATASET_XFER);
00925 #ifdef HAVE_MPMI
00926   H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
00927 #endif
00928 
00929   // Select hyperslab in the file.
00930   hsize_t offset[] = {0, LinearX->Map().GID(0)};
00931   hsize_t stride[] = {1, 1};
00932   hsize_t count[] = {NumVectors, 1};
00933   hsize_t block[] = {1, LinearX->MyLength()};
00934   filespace_id = H5Dget_space(dset_id);
00935   H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride, 
00936                       count, block);
00937 
00938   // Each process defines dataset in memory and writes it to the hyperslab in the file.
00939   hsize_t dimsm[] = {NumVectors * LinearX->MyLength()};
00940   memspace_id = H5Screate_simple(1, dimsm, NULL);
00941 
00942   // Write hyperslab
00943   status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id, 
00944                     plist_id_, LinearX->Values());
00945   CHECK_STATUS(status);
00946   H5Gclose(group_id);
00947   H5Sclose(memspace_id);
00948   H5Sclose(filespace_id);
00949   H5Dclose(dset_id);
00950   H5Pclose(plist_id_);
00951 
00952   Write(GroupName, "GlobalLength", GlobalLength);
00953   Write(GroupName, "NumVectors", NumVectors);
00954   Write(GroupName, "__type__", "Epetra_MultiVector");
00955 }
00956 
00957 // ==========================================================================
00958 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
00959                         Epetra_MultiVector*& X)
00960 {
00961   // first read it with linear distribution
00962   Epetra_MultiVector* LinearX;
00963   Read(GroupName, LinearX);
00964   
00965   // now build the importer to the actual one
00966   Epetra_Import Importer(Map, X->Map());
00967   X = new Epetra_MultiVector(Map, LinearX->NumVectors());
00968   X->Import(*LinearX, Importer, Insert);
00969 
00970   // delete memory
00971   delete LinearX;
00972 }
00973 
00974 // ==========================================================================
00975 void EpetraExt::HDF5::Read(const std::string& GroupName, Epetra_MultiVector*& LinearX)
00976 {
00977   int GlobalLength, NumVectors;
00978 
00979   ReadMultiVectorProperties(GroupName, GlobalLength, NumVectors);
00980 
00981   hid_t group_id;
00982 
00983   // Create the dataspace for the dataset.
00984   hsize_t q_dimsf[] = {NumVectors, GlobalLength};
00985   hid_t filespace_id = H5Screate_simple(2, q_dimsf, NULL);
00986 
00987   if (!IsContained(GroupName))
00988     throw(Exception(__FILE__, __LINE__,
00989                     "requested group " + GroupName + " not found"));
00990 
00991   group_id = H5Gopen(file_id_, GroupName.c_str());
00992 
00993   // Create the dataset with default properties and close filespace_id.
00994   hid_t dset_id = H5Dopen(group_id, "Values");
00995 
00996   // Create property list for collective dataset write.
00997   plist_id_ = H5Pcreate(H5P_DATASET_XFER);
00998 #ifdef HAVE_MPI
00999   H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
01000 #endif
01001   H5Pclose(plist_id_);
01002 
01003   Epetra_Map LinearMap(GlobalLength, 0, Comm());
01004   LinearX = new Epetra_MultiVector(LinearMap, NumVectors);
01005 
01006   // Select hyperslab in the file.
01007   hsize_t offset[] = {0, LinearMap.GID(0)};
01008   hsize_t stride[] = {1, 1};
01009   hsize_t count[] = {NumVectors, 1};
01010   hsize_t block[] = {1, LinearX->MyLength()};
01011 
01012   filespace_id = H5Dget_space(dset_id);
01013   H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, offset, stride, 
01014                       count, block);
01015 
01016   // Each process defines dataset in memory and writes it to the hyperslab in the file.
01017   hsize_t dimsm[] = {NumVectors * LinearX->MyLength()};
01018   hid_t memspace_id = H5Screate_simple(1, dimsm, NULL);
01019 
01020   // Write hyperslab
01021   CHECK_STATUS(H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace_id, filespace_id, 
01022                        H5P_DEFAULT, LinearX->Values()));
01023 
01024   CHECK_STATUS(H5Gclose(group_id));
01025   CHECK_STATUS(H5Sclose(memspace_id));
01026   CHECK_STATUS(H5Sclose(filespace_id));
01027   CHECK_STATUS(H5Dclose(dset_id));
01028 }
01029 
01030 // ==========================================================================
01031 void EpetraExt::HDF5::ReadMultiVectorProperties(const std::string& GroupName, 
01032                                              int& GlobalLength,
01033                                              int& NumVectors)
01034 {
01035   if (!IsContained(GroupName))
01036     throw(Exception(__FILE__, __LINE__,
01037                     "requested group " + GroupName + " not found"));
01038 
01039   std::string Label;
01040   Read(GroupName, "__type__", Label);
01041 
01042   if (Label != "Epetra_MultiVector")
01043     throw(Exception(__FILE__, __LINE__,
01044                     "requested group " + GroupName + " is not an Epetra_MultiVector",
01045                     "__type__ = " + Label));
01046 
01047   Read(GroupName, "GlobalLength", GlobalLength);
01048   Read(GroupName, "NumVectors", NumVectors);
01049 }
01050 
01051 // ========================= //
01052 // EpetraExt::DistArray<int> //
01053 // ========================= //
01054 
01055 // ==========================================================================
01056 void EpetraExt::HDF5::Write(const std::string& GroupName, const DistArray<int>& x)
01057 {
01058   if (x.Map().LinearMap())
01059   {
01060     Write(GroupName, "Values", x.Map().NumMyElements() * x.RowSize(), 
01061           x.Map().NumGlobalElements() * x.RowSize(),
01062           H5T_NATIVE_INT, x.Values());
01063   }
01064   else
01065   {
01066     // need to build a linear map first, the import data, then
01067     // finally write them 
01068     const Epetra_BlockMap& OriginalMap = x.Map();
01069     Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
01070     Epetra_Import Importer(LinearMap, OriginalMap);
01071 
01072     EpetraExt::DistArray<int> LinearX(LinearMap, x.RowSize());
01073     LinearX.Import(x, Importer, Insert);
01074 
01075     Write(GroupName, "Values", LinearMap.NumMyElements() * x.RowSize(), 
01076           LinearMap.NumGlobalElements() * x.RowSize(),
01077           H5T_NATIVE_INT, LinearX.Values());
01078   }
01079 
01080   Write(GroupName, "__type__", "EpetraExt::DistArray<int>");
01081   Write(GroupName, "GlobalLength", x.GlobalLength());
01082   Write(GroupName, "RowSize", x.RowSize());
01083 }
01084 
01085 // ==========================================================================
01086 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
01087                            DistArray<int>*& X)
01088 {
01089   // gets the length of the std::vector
01090   int GlobalLength, RowSize;
01091   ReadIntDistArrayProperties(GroupName, GlobalLength, RowSize);
01092 
01093   if (Map.LinearMap())
01094   {
01095     X = new EpetraExt::DistArray<int>(Map, RowSize);
01096     // simply read stuff and go home
01097     Read(GroupName, "Values", Map.NumMyElements() * RowSize, 
01098          Map.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->Values());
01099   }
01100   else
01101   {
01102     // we need to first create a linear map, read the std::vector,
01103     // then import it to the actual nonlinear map
01104     Epetra_Map LinearMap(GlobalLength, 0, Comm_);
01105     EpetraExt::DistArray<int> LinearX(LinearMap, RowSize);
01106 
01107     Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize, 
01108          LinearMap.NumGlobalElements() * RowSize,
01109          H5T_NATIVE_INT, LinearX.Values());
01110 
01111     Epetra_Import Importer(Map, LinearMap);
01112     X = new EpetraExt::DistArray<int>(Map, RowSize);
01113     X->Import(LinearX, Importer, Insert);
01114   }
01115 }
01116 
01117 // ==========================================================================
01118 void EpetraExt::HDF5::Read(const std::string& GroupName, DistArray<int>*& X)
01119 {
01120   // gets the length of the std::vector
01121   int GlobalLength, RowSize;
01122   ReadIntDistArrayProperties(GroupName, GlobalLength, RowSize);
01123 
01124   // creates a new linear map and use it to read data
01125   Epetra_Map LinearMap(GlobalLength, 0, Comm_);
01126   X = new EpetraExt::DistArray<int>(LinearMap, RowSize);
01127 
01128   Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize, 
01129        LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_INT, X->Values());
01130 }
01131 
01132 // ==========================================================================
01133 void EpetraExt::HDF5::ReadIntDistArrayProperties(const std::string& GroupName, 
01134                                                  int& GlobalLength,
01135                                                  int& RowSize)
01136 {
01137   if (!IsContained(GroupName))
01138     throw(Exception(__FILE__, __LINE__,
01139                     "requested group " + GroupName + " not found"));
01140 
01141   std::string Label;
01142   Read(GroupName, "__type__", Label);
01143 
01144   if (Label != "EpetraExt::DistArray<int>")
01145     throw(Exception(__FILE__, __LINE__,
01146                     "requested group " + GroupName + " is not an EpetraExt::DistArray<int>",
01147                     "__type__ = " + Label));
01148 
01149   Read(GroupName, "GlobalLength", GlobalLength);
01150   Read(GroupName, "RowSize", RowSize);
01151 }
01152 
01153 // ============================ //
01154 // EpetraExt::DistArray<double> //
01155 // ============================ //
01156 
01157 // ==========================================================================
01158 void EpetraExt::HDF5::Write(const std::string& GroupName, const DistArray<double>& x)
01159 {
01160   if (x.Map().LinearMap())
01161   {
01162     Write(GroupName, "Values", x.Map().NumMyElements() * x.RowSize(), 
01163           x.Map().NumGlobalElements() * x.RowSize(),
01164           H5T_NATIVE_DOUBLE, x.Values());
01165   }
01166   else
01167   {
01168     // need to build a linear map first, the import data, then
01169     // finally write them 
01170     const Epetra_BlockMap& OriginalMap = x.Map();
01171     Epetra_Map LinearMap(OriginalMap.NumGlobalElements(), OriginalMap.IndexBase(), Comm_);
01172     Epetra_Import Importer(LinearMap, OriginalMap);
01173 
01174     EpetraExt::DistArray<double> LinearX(LinearMap, x.RowSize());
01175     LinearX.Import(x, Importer, Insert);
01176 
01177     Write(GroupName, "Values", LinearMap.NumMyElements() * x.RowSize(), 
01178           LinearMap.NumGlobalElements() * x.RowSize(),
01179           H5T_NATIVE_DOUBLE, LinearX.Values());
01180   }
01181 
01182   Write(GroupName, "__type__", "EpetraExt::DistArray<double>");
01183   Write(GroupName, "GlobalLength", x.GlobalLength());
01184   Write(GroupName, "RowSize", x.RowSize());
01185 }
01186 
01187 // ==========================================================================
01188 void EpetraExt::HDF5::Read(const std::string& GroupName, const Epetra_Map& Map,
01189                            DistArray<double>*& X)
01190 {
01191   // gets the length of the std::vector
01192   int GlobalLength, RowSize;
01193   ReadDoubleDistArrayProperties(GroupName, GlobalLength, RowSize);
01194 
01195   if (Map.LinearMap())
01196   {
01197     X = new EpetraExt::DistArray<double>(Map, RowSize);
01198     // simply read stuff and go home
01199     Read(GroupName, "Values", Map.NumMyElements() * RowSize, 
01200          Map.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->Values());
01201   }
01202   else
01203   {
01204     // we need to first create a linear map, read the std::vector,
01205     // then import it to the actual nonlinear map
01206     Epetra_Map LinearMap(GlobalLength, 0, Comm_);
01207     EpetraExt::DistArray<double> LinearX(LinearMap, RowSize);
01208 
01209     Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize, 
01210          LinearMap.NumGlobalElements() * RowSize,
01211          H5T_NATIVE_DOUBLE, LinearX.Values());
01212 
01213     Epetra_Import Importer(Map, LinearMap);
01214     X = new EpetraExt::DistArray<double>(Map, RowSize);
01215     X->Import(LinearX, Importer, Insert);
01216   }
01217 }
01218 
01219 // ==========================================================================
01220 void EpetraExt::HDF5::Read(const std::string& GroupName, DistArray<double>*& X)
01221 {
01222   // gets the length of the std::vector
01223   int GlobalLength, RowSize;
01224   ReadDoubleDistArrayProperties(GroupName, GlobalLength, RowSize);
01225 
01226   // creates a new linear map and use it to read data
01227   Epetra_Map LinearMap(GlobalLength, 0, Comm_);
01228   X = new EpetraExt::DistArray<double>(LinearMap, RowSize);
01229 
01230   Read(GroupName, "Values", LinearMap.NumMyElements() * RowSize, 
01231        LinearMap.NumGlobalElements() * RowSize, H5T_NATIVE_DOUBLE, X->Values());
01232 }
01233 //
01234 // ==========================================================================
01235 void EpetraExt::HDF5::ReadDoubleDistArrayProperties(const std::string& GroupName, 
01236                                                     int& GlobalLength,
01237                                                     int& RowSize)
01238 {
01239   if (!IsContained(GroupName))
01240     throw(Exception(__FILE__, __LINE__,
01241                     "requested group " + GroupName + " not found"));
01242 
01243   std::string Label;
01244   Read(GroupName, "__type__", Label);
01245 
01246   if (Label != "EpetraExt::DistArray<double>")
01247     throw(Exception(__FILE__, __LINE__,
01248                     "requested group " + GroupName + " is not an EpetraExt::DistArray<double>",
01249                     "__type__ = " + Label));
01250 
01251   Read(GroupName, "GlobalLength", GlobalLength);
01252   Read(GroupName, "RowSize", RowSize);
01253 }
01254 
01255 // ======================= //
01256 // basic serial data types //
01257 // ======================= //
01258 
01259 // ==========================================================================
01260 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
01261                             int what)
01262 {
01263   if (!IsContained(GroupName))
01264     CreateGroup(GroupName);
01265 
01266   hid_t filespace_id = H5Screate(H5S_SCALAR);
01267   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01268   hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_INT, 
01269                       filespace_id, H5P_DEFAULT);
01270 
01271   herr_t status = H5Dwrite(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
01272                            H5P_DEFAULT, &what);
01273   CHECK_STATUS(status);
01274 
01275   // Close/release resources.
01276   H5Dclose(dset_id);
01277   H5Gclose(group_id);
01278   H5Sclose(filespace_id);
01279 }
01280 
01281 // ==========================================================================
01282 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
01283                             double what)
01284 {
01285   if (!IsContained(GroupName))
01286     CreateGroup(GroupName);
01287 
01288   hid_t filespace_id = H5Screate(H5S_SCALAR);
01289   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01290   hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), H5T_NATIVE_DOUBLE, 
01291                             filespace_id, H5P_DEFAULT);
01292 
01293   herr_t status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, 
01294                            filespace_id, H5P_DEFAULT, &what);
01295   CHECK_STATUS(status);
01296 
01297   // Close/release resources.
01298   H5Dclose(dset_id);
01299   H5Gclose(group_id);
01300   H5Sclose(filespace_id);
01301 }
01302 
01303 // ==========================================================================
01304 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName, int& data)
01305 {
01306   if (!IsContained(GroupName))
01307     throw(Exception(__FILE__, __LINE__,
01308                     "requested group " + GroupName + " not found"));
01309 
01310   // Create group in the root group using absolute name.
01311   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01312 
01313   hid_t filespace_id = H5Screate(H5S_SCALAR);
01314   hid_t dset_id = H5Dopen(group_id, DataSetName.c_str());
01315 
01316   herr_t status = H5Dread(dset_id, H5T_NATIVE_INT, H5S_ALL, filespace_id,
01317                     H5P_DEFAULT, &data);
01318   CHECK_STATUS(status);
01319 
01320   H5Sclose(filespace_id);
01321   H5Dclose(dset_id);
01322   H5Gclose(group_id);
01323 }
01324 
01325 // ==========================================================================
01326 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName, double& data)
01327 {
01328   if (!IsContained(GroupName))
01329     throw(Exception(__FILE__, __LINE__,
01330                     "requested group " + GroupName + " not found"));
01331 
01332   // Create group in the root group using absolute name.
01333   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01334 
01335   hid_t filespace_id = H5Screate(H5S_SCALAR);
01336   hid_t dset_id = H5Dopen(group_id, DataSetName.c_str());
01337 
01338   herr_t status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, filespace_id,
01339                     H5P_DEFAULT, &data);
01340   CHECK_STATUS(status);
01341 
01342   H5Sclose(filespace_id);
01343   H5Dclose(dset_id);
01344   H5Gclose(group_id);
01345 }
01346 
01347 // ==========================================================================
01348 void EpetraExt::HDF5::Write(const std::string& GroupName, 
01349                             const std::string& DataSetName, 
01350                             const std::string& data)
01351 {
01352   if (!IsContained(GroupName))
01353     CreateGroup(GroupName);
01354 
01355   hsize_t len = 1;
01356 
01357   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01358 
01359   hid_t dataspace_id = H5Screate_simple(1, &len, NULL);
01360 
01361   hid_t atype = H5Tcopy(H5T_C_S1);
01362   H5Tset_size(atype, data.size() + 1);
01363 
01364   hid_t dataset_id = H5Dcreate(group_id, DataSetName.c_str(), atype, 
01365                                dataspace_id, H5P_DEFAULT);
01366   CHECK_STATUS(H5Dwrite(dataset_id, atype, H5S_ALL, H5S_ALL, 
01367                         H5P_DEFAULT, data.c_str()));
01368 
01369   CHECK_STATUS(H5Tclose(atype));
01370 
01371   CHECK_STATUS(H5Dclose(dataset_id));
01372 
01373   CHECK_STATUS(H5Sclose(dataspace_id));
01374 
01375   CHECK_STATUS(H5Gclose(group_id));
01376 }
01377 
01378 // ==========================================================================
01379 void EpetraExt::HDF5::Read(const std::string& GroupName, 
01380                            const std::string& DataSetName, 
01381                            std::string& data)
01382 {
01383   if (!IsContained(GroupName))
01384     throw(Exception(__FILE__, __LINE__,
01385                     "requested group " + GroupName + " not found"));
01386 
01387   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01388 
01389   hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str());
01390 
01391   hid_t datatype_id = H5Dget_type(dataset_id);
01392   size_t typesize_id = H5Tget_size(datatype_id);
01393   H5T_class_t typeclass_id = H5Tget_class(datatype_id);
01394 
01395   if(typeclass_id != H5T_STRING) 
01396     throw(Exception(__FILE__, __LINE__,
01397                     "requested group " + GroupName + " is not a std::string"));
01398   char data2[160];
01399   CHECK_STATUS(H5Dread(dataset_id, datatype_id, H5S_ALL, H5S_ALL,
01400                        H5P_DEFAULT, data2));
01401   data = data2;
01402 
01403   H5Dclose(dataset_id);  
01404   H5Gclose(group_id);  
01405 }
01406 
01407 // ============= //
01408 // serial arrays //
01409 // ============= //
01410 
01411 // ==========================================================================
01412 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
01413                          const int type, const int Length, 
01414                          void* data)
01415 {
01416   if (!IsContained(GroupName))
01417     CreateGroup(GroupName);
01418 
01419   hsize_t dimsf = Length;
01420 
01421   hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
01422 
01423   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01424 
01425   hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type, 
01426                             filespace_id, H5P_DEFAULT);
01427 
01428   herr_t status = H5Dwrite(dset_id, type, H5S_ALL, H5S_ALL,
01429                            H5P_DEFAULT, data);
01430   CHECK_STATUS(status);
01431 
01432   // Close/release resources.
01433   H5Dclose(dset_id);
01434   H5Gclose(group_id);
01435   H5Sclose(filespace_id);
01436 }
01437 
01438 // ==========================================================================
01439 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName,
01440                         const int type, const int Length, 
01441                         void* data)
01442 {
01443   if (!IsContained(GroupName))
01444     throw(Exception(__FILE__, __LINE__,
01445                     "requested group " + GroupName + " not found"));
01446 
01447   hsize_t dimsf = Length;
01448 
01449   hid_t filespace_id = H5Screate_simple(1, &dimsf, NULL);
01450 
01451   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01452 
01453   hid_t dset_id = H5Dopen(group_id, DataSetName.c_str());
01454 
01455   herr_t status = H5Dread(dset_id, type, H5S_ALL, filespace_id,
01456                           H5P_DEFAULT, data);
01457   CHECK_STATUS(status);
01458 
01459   // Close/release resources.
01460   H5Dclose(dset_id);
01461   H5Gclose(group_id);
01462   H5Sclose(filespace_id);
01463 }
01464 
01465 // ================== //
01466 // distributed arrays //
01467 // ================== //
01468 
01469 // ==========================================================================
01470 void EpetraExt::HDF5::Write(const std::string& GroupName, const std::string& DataSetName,
01471                          int MySize, int GlobalSize, int type, const void* data)
01472 {
01473   int Offset;
01474   Comm_.ScanSum(&MySize, &Offset, 1);
01475   Offset -= MySize;
01476 
01477   hsize_t MySize_t = MySize;
01478   hsize_t GlobalSize_t = GlobalSize;
01479   hsize_t Offset_t = Offset;
01480 
01481   hid_t filespace_id = H5Screate_simple(1, &GlobalSize_t, NULL); 
01482 
01483   // Create the dataset with default properties and close filespace.
01484   if (!IsContained(GroupName))
01485     CreateGroup(GroupName);
01486 
01487   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01488   hid_t dset_id = H5Dcreate(group_id, DataSetName.c_str(), type, filespace_id, 
01489                             H5P_DEFAULT);
01490 
01491   H5Sclose(filespace_id);
01492 
01493   // Each process defines dataset in memory and writes it to the hyperslab
01494   // in the file.
01495 
01496   hid_t memspace_id = H5Screate_simple(1, &MySize_t, NULL);
01497 
01498   // Select hyperslab in the file.
01499   filespace_id = H5Dget_space(dset_id);
01500   H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL, &MySize_t, NULL);
01501 
01502   plist_id_ = H5Pcreate(H5P_DATASET_XFER);
01503 #ifdef HAVE_MPI
01504   H5Pset_dxpl_mpio(plist_id_, H5FD_MPIO_COLLECTIVE);
01505 #endif
01506 
01507   status = H5Dwrite(dset_id, type, memspace_id, filespace_id,
01508                     plist_id_, data);
01509   CHECK_STATUS(status);
01510 
01511   // Close/release resources.
01512   H5Dclose(dset_id);
01513   H5Gclose(group_id);
01514   H5Sclose(filespace_id);
01515   H5Sclose(memspace_id);
01516   H5Pclose(plist_id_);
01517 }
01518 
01519 // ==========================================================================
01520 void EpetraExt::HDF5::Read(const std::string& GroupName, const std::string& DataSetName,
01521                         int MySize, int GlobalSize,
01522                         const int type, void* data)
01523 {
01524   if (!IsOpen())
01525     throw(Exception(__FILE__, __LINE__, "no file open yet"));
01526 
01527   hsize_t MySize_t = MySize;
01528 
01529   // offset
01530   int itmp;
01531   Comm_.ScanSum(&MySize, &itmp, 1);
01532   hsize_t Offset_t = itmp - MySize;
01533 
01534   hid_t group_id = H5Gopen(file_id_, GroupName.c_str());
01535   hid_t dataset_id = H5Dopen(group_id, DataSetName.c_str());
01536   //hid_t space_id = H5Screate_simple(1, &Offset_t, 0);
01537 
01538   // Select hyperslab in the file.
01539   hid_t filespace_id = H5Dget_space(dataset_id);
01540   H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET, &Offset_t, NULL, 
01541                       &MySize_t, NULL);
01542 
01543   hid_t mem_dataspace = H5Screate_simple (1, &MySize_t, NULL);
01544 
01545   herr_t status = H5Dread(dataset_id, type, mem_dataspace, filespace_id, 
01546                           H5P_DEFAULT, data);
01547   CHECK_STATUS(status);
01548 
01549   H5Sclose(mem_dataspace);
01550   H5Gclose(group_id);  
01551   //H5Sclose(space_id);  
01552   H5Dclose(dataset_id);  
01553 //  H5Dclose(filespace_id);  
01554 }
01555 #endif

Generated on Tue Oct 20 12:45:29 2009 for EpetraExt by doxygen 1.4.7