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