EpetraExt Package Browser (Single Doxygen Collection) Development
HDF5_IO.cpp
Go to the documentation of this file.
00001 #include "EpetraExt_ConfigDefs.h"
00002 #ifdef HAVE_MPI
00003 #include "mpi.h"
00004 #include "Epetra_MpiComm.h"
00005 #else
00006 #include "Epetra_SerialComm.h"
00007 #endif
00008 #include <vector>
00009 #include "Epetra_Map.h"
00010 #include "Epetra_CrsMatrix.h"
00011 #include "Epetra_MultiVector.h"
00012 #include "EpetraExt_HDF5.h"
00013 #include "EpetraExt_Utils.h"
00014 #include "EpetraExt_Exception.h"
00015 #include "Teuchos_ParameterList.hpp"
00016 
00017 // Showing the usage of HDF5 I/O.
00018 // This example can be run with any number of processors.
00019 //
00020 // \author Marzio Sala, D-INFK/ETHZ.
00021 //
00022 // \date Last modified on 09-Mar-06.
00023 
00024 int main (int argc, char **argv)
00025 {
00026 #ifdef HAVE_MPI
00027   MPI_Init(&argc, &argv);
00028   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00029 #else
00030   Epetra_SerialComm Comm;
00031 #endif
00032 
00033   try 
00034   {
00035     int n = Comm.NumProc() * 4;
00036     Epetra_Map Map(n, 0, Comm);
00037     Epetra_MultiVector x(Map, 2); x.Random();
00038     Epetra_MultiVector b(Map, 2); x.Random();
00039     Epetra_CrsMatrix Matrix(Copy, Map, 0);
00040     // diagonal matrix
00041     for (int i = 0; i < Map.NumMyElements(); ++i)
00042     {
00043       int ii = Map.GID(i);
00044       double one = 1.0;
00045       Matrix.InsertGlobalValues(ii, 1, &one, &ii);
00046     }
00047     Matrix.FillComplete();
00048 
00049     // create a Teuchos::ParameterList and populate it
00050     Teuchos::ParameterList List;
00051     List.set("bool type", true);
00052     List.set("int type", 2);
00053     List.set("double type", 3.0);
00054     List.set("std::string type", "a std::string");
00055 
00056     // This is the HDF5 file manager
00057     EpetraExt::HDF5 HDF5(Comm);
00058 
00059     // creates a new file. To open an existing file, use Open("myfile.h5")
00060     HDF5.Create("myfile.h5");
00061 
00062     // =========================== //
00063     // P A R T   I:  W R I T I N G //
00064     // =========================== //
00065 
00066     if (Comm.MyPID() == 0)
00067       std::cout << "Writing objects to HDF5 file myfile.h5..." << std::endl << std::endl;
00068 
00069     // We first write the map, whose name contains the number of processors
00070     HDF5.Write("map-" + EpetraExt::toString(Comm.NumProc()), Map);
00071     // Then we write the matrix...
00072     HDF5.Write("matrix", Matrix);
00073     // ... and the x, b vectors
00074     HDF5.Write("x", x);
00075     HDF5.Write("b", b);
00076     // we can associate basic data types with a given group, for example:
00077     HDF5.Write("matrix", "integration order", 1);
00078     HDF5.Write("matrix", "numerical drop", 0.1);
00079     HDF5.Write("matrix", "package name", "EpetraExt");
00080     HDF5.Write("matrix", "author", "Marzio Sala");
00081     HDF5.Write("matrix", "institution", "ETHZ/D-INFK");
00082     // or put them in a new group
00083     HDF5.Write("my parameters", "latitude", 12);
00084     HDF5.Write("my parameters", "longitude", 67);
00085 
00086     // Now we write the parameter list we have used to create the matrix
00087     HDF5.Write("List", List);
00088     // We can also write integers/doubles/arrays. All these quantities are
00089     // supposed to have the same value on all processors.
00090     std::vector<int> iarray(3); 
00091     iarray[0] = 0, iarray[1] = 1; iarray[2] = 2;
00092     HDF5.Write("my parameters", "int array", H5T_NATIVE_INT, 3, &iarray[0]);
00093 
00094     std::vector<double> darray(3); 
00095     darray[0] = 0.1, darray[1] = 1.1; darray[2] = 2.1;
00096     HDF5.Write("my parameters", "double array", H5T_NATIVE_DOUBLE, 3, &darray[0]);
00097 
00098     // To analyze the content of the file, one can use 
00099     // "h5dump filename.h5" or "h5dump filename.h5 -H" o
00100 
00101     // ============================ //
00102     // P A R T   II:  R E A D I N G //
00103     // ============================ //
00104 
00105     if (Comm.MyPID() == 0)
00106       std::cout << "Reading objects from HDF5 file myfile.h5..." << std::endl << std::endl;
00107 
00108     Epetra_Map* NewMap = 0;
00109     Epetra_CrsMatrix* NewMatrix = 0;
00110     Epetra_MultiVector* NewX,* NewB;
00111 
00112     // Check if the map is there (in this case it is). If it is, read the
00113     // matrix using the map, otherwise read the matrix with a linear map.
00114     if (HDF5.IsContained("map-" + EpetraExt::toString(Comm.NumProc())))
00115     {
00116       HDF5.Read("map-" + EpetraExt::toString(Comm.NumProc()), NewMap);
00117       HDF5.Read("matrix", *NewMap, *NewMap, NewMatrix);
00118     }
00119     else
00120     {
00121       int NumGlobalRows;
00122       HDF5.Read("matrix", "NumGlobalRows", NumGlobalRows);
00123       NewMap = new Epetra_Map(NumGlobalRows, 0, Comm);
00124       HDF5.Read("matrix", *NewMap, *NewMap, NewMatrix);
00125     }
00126 
00127     // read the number of nonzeros from file, compare them with the
00128     // actual number of NewMatrix
00129     int NewNumGlobalNonzeros;
00130     HDF5.Read("matrix", "NumGlobalNonzeros", NewNumGlobalNonzeros);
00131 
00132     assert (Matrix.NumGlobalNonzeros() == NewNumGlobalNonzeros);
00133 
00134     // Now read the MultiVector's
00135     HDF5.Read("x", NewX);
00136     HDF5.Read("b", NewB);
00137 
00138     // the int/double values, and the int/double arrays
00139     int new_latitude, new_longitude;
00140     std::vector<int>    new_iarray(3);
00141     std::vector<double> new_darray(3);
00142     HDF5.Read("my parameters", "latitude", new_latitude);
00143     HDF5.Read("my parameters", "longitude", new_longitude);
00144     HDF5.Read("my parameters", "int array", H5T_NATIVE_INT, 3, &new_iarray[0]);
00145     HDF5.Read("my parameters", "double array", H5T_NATIVE_DOUBLE, 3, &new_darray[0]);
00146     // and the std::string values associated with group "matrix", recordered
00147     // with dataset "package name".
00148     std::string PackageName;
00149     HDF5.Read("matrix", "package name", PackageName);
00150 
00151     Teuchos::ParameterList newList;
00152     HDF5.Read("List", newList);
00153     if (Comm.MyPID() == 0)
00154     {
00155       std::cout << "New list as read from file is:" << std::endl;
00156       std::cout << "bool type = " << newList.get("bool type", false) << std::endl;
00157       std::cout << "int type = " << newList.get("int type", -1) << std::endl;
00158       std::cout << "double type = " << newList.get("double type", -1.0) << std::endl;
00159       std::cout << "std::string type = " << newList.get("std::string type", "not-set") << std::endl;
00160       std::cout << std::endl;
00161 
00162       std::cout << "Checking some read and written data..." << std::endl;
00163       for (int i = 0; i < 3; ++i)
00164         std::cout << "iarray[" << i << "] = " << iarray[i] << " should be " << new_iarray[i] << std::endl;
00165       for (int i = 0; i < 3; ++i)
00166         std::cout << "darray[" << i << "] = " << darray[i] << " should be " << new_darray[i] << std::endl;
00167       std::cout << std::endl;
00168       std::cout << "Try to print out the content of myfile.h5 using the command" << std::endl;
00169       std::cout << "    h5dump myfile.h5" << std::endl;
00170     }
00171 
00172     // We finally close the file. Better to close it before calling
00173     // MPI_Finalize() to avoid MPI-related errors, since Close() might call MPI
00174     // functions.
00175     HDF5.Close();
00176 
00177     // delete memory
00178     if (NewMap   ) delete NewMap;
00179     if (NewMatrix) delete NewMatrix;
00180     if (NewX     ) delete NewX;
00181     if (NewB     ) delete NewB;
00182   }
00183   catch(EpetraExt::Exception& rhs) 
00184   {
00185     rhs.Print();
00186   }
00187   catch (...) 
00188   {
00189     std::cerr << "Caught generic std::exception" << std::endl;
00190   }
00191 
00192 
00193 #ifdef HAVE_MPI
00194   MPI_Finalize();
00195 #endif
00196 
00197   return(EXIT_SUCCESS);
00198 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines