EpetraExt_HDF5.cpp

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

Generated on Tue Jul 13 09:23:06 2010 for EpetraExt by  doxygen 1.4.7