improve documentation and minor clean-ups

diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h
index 584e714b1c155312350f1f544e42cb8441483a4d..5abfedd1069271b47cb8fb2648e722c0465bf0ac 100644 (file)
--- a/include/deal.II/lac/scalapack.h
+++ b/include/deal.II/lac/scalapack.h
@@ -177,6 +177,9 @@ public:
    * - The dimension of the submatrix to be copied is given by @p submatrix_size
    *   with number of rows=<code>submatrix_size.first</code> and number of columns=<code>submatrix_size.second</code>
    * .
+   *
+   * If it is necessary to copy complete matrices with an identical block-cyclic distribution,
+   * use copy_to(ScaLAPACKMatrix<NumberType> &dest) with only one argument to avoid communication
    */
   void copy_to(ScaLAPACKMatrix<NumberType> &B,
                const std::pair<unsigned int,unsigned int> &offset_A,
@@ -185,13 +188,23 @@ public:
 
   /**
    * Stores the distributed matrix in @p filename using HDF5
+   *
+   * If HDF5 was build with MPI, parallel I/O is used to save the matrix.
+   * Otherwise, just one process will do the output.
    */
   void save(const char *filename) const;
 
   /**
-   * Loads the distributed matrix from file @p filename using HDF5
+   * Loads the distributed matrix from file @p filename using HDF5.
+   *
+   * The matrix must have the same dimensions as the matrix in stored in the file.
+   *
+   * If HDF5 was build with MPI, parallel I/O is used to load the matrix.
+   * Otherwise, just one process will load the matrix from storage
+   * and distribute the content to the other processes subsequently.
    */
   void load(const char *filename);
+
   /**
    * Compute the Cholesky factorization of the matrix using ScaLAPACK
    * function <code>pXpotrf</code>. The result of the factorization is stored in this object.

diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc
index 6879cec96ce4d4d0562ea66e1412ae1766df957d..ffafdf5d76012786cb25b20419ddb1db46d37c91 100644 (file)
--- a/source/lac/scalapack.cc
+++ b/source/lac/scalapack.cc
@@ -23,9 +23,6 @@
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/mpi.templates.h>
 #include <deal.II/lac/scalapack.templates.h>
-#include <deal.II/base/conditional_ostream.h>
-
-#include <array>
 
 #  ifdef DEAL_II_WITH_HDF5
 #include <hdf5.h>
@@ -935,7 +932,7 @@ void ScaLAPACKMatrix<NumberType>::save_serial(const char *filename) const
       // create the dataset
       hid_t type_id = hdf5_type_id(&tmp.values[0]);
       hid_t dataset_id = H5Dcreate2(file_id, "/matrix",
-                                                               type_id, dataspace_id,
+                                    type_id, dataspace_id,
                                     H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
       // write the dataset
@@ -983,14 +980,14 @@ void ScaLAPACKMatrix<NumberType>::save_parallel(const char *filename) const
   copy_to(tmp);
 
   // get pointer to data held by the process
-  NumberType* data = (tmp.values.size()>0) ? &tmp.values[0] : nullptr;
+  NumberType *data = (tmp.values.size()>0) ? &tmp.values[0] : nullptr;
 
   herr_t status;
   // dataset dimensions
   hsize_t dims[2];
 
   // set up file access property list with parallel I/O access
-  hid_t        plist_id = H5Pcreate(H5P_FILE_ACCESS);
+  hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
   status = H5Pset_fapl_mpio(plist_id, tmp.grid->mpi_communicator, info);
   AssertThrow(status >= 0, ExcIO());
 
@@ -1011,7 +1008,7 @@ void ScaLAPACKMatrix<NumberType>::save_parallel(const char *filename) const
   // create the dataset with default properties and close filespace
   hid_t type_id = hdf5_type_id(data);
   hid_t dset_id = H5Dcreate2(file_id, "/matrix", type_id,
-                                                    filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+                             filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
   status = H5Sclose(filespace);
   AssertThrow(status >= 0, ExcIO());
 
@@ -1031,7 +1028,7 @@ void ScaLAPACKMatrix<NumberType>::save_parallel(const char *filename) const
 
   hsize_t offset[2] = {0};
   for (unsigned int i=0; i<my_rank; ++i)
-         offset[0] += proc_n_local_columns[i];
+    offset[0] += proc_n_local_columns[i];
 
   // select hyperslab in the file.
   filespace = H5Dget_space(dset_id);
@@ -1045,11 +1042,11 @@ void ScaLAPACKMatrix<NumberType>::save_parallel(const char *filename) const
 
   // process with no data will not participate in writing to the file
   if (tmp.values.size()>0)
-  {
-         status = H5Dwrite(dset_id, type_id, memspace, filespace,
-                                       plist_id, data);
-         AssertThrow(status >= 0, ExcIO());
-  }
+    {
+      status = H5Dwrite(dset_id, type_id, memspace, filespace,
+                        plist_id, data);
+      AssertThrow(status >= 0, ExcIO());
+    }
 
   // close/release sources
   status = H5Dclose(dset_id);
@@ -1191,88 +1188,88 @@ void ScaLAPACKMatrix<NumberType>::load_parallel(const char *filename)
   ScaLAPACKMatrix<NumberType> tmp(n_rows,n_columns,column_grid,MB,NB);
 
   // get pointer to data held by the process
-  NumberType* data = (tmp.values.size()>0) ? &tmp.values[0] : nullptr;
+  NumberType *data = (tmp.values.size()>0) ? &tmp.values[0] : nullptr;
 
   herr_t status;
 
-   // set up file access property list with parallel I/O access
-   hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
-   status = H5Pset_fapl_mpio(plist_id, tmp.grid->mpi_communicator, info);
-   AssertThrow(status >= 0, ExcIO());
-
-   // open file collectively in read-only mode and release property list identifier
-   hid_t file_id = H5Fopen(filename, H5F_ACC_RDONLY, plist_id);
-   status = H5Pclose(plist_id);
-   AssertThrow(status >= 0, ExcIO());
-
-   // open the dataset in the file collectively
-   hid_t dataset_id = H5Dopen2(file_id, "/matrix", H5P_DEFAULT);
-
-   // check the datatype of the dataset in the file
-   // if the classes of type of the dataset and the matrix do not match abort
-   // see HDF User's Guide: 6.10. Data Transfer: Datatype Conversion and Selection
-   hid_t datatype = hdf5_type_id(data);
-   hid_t datatype_inp  = H5Dget_type(dataset_id);
-   H5T_class_t t_class_inp = H5Tget_class(datatype_inp);
-   H5T_class_t t_class = H5Tget_class(datatype);
-   AssertThrow(t_class_inp == t_class,
-               ExcMessage("The data type of the matrix to be read does not match the archive"));
-
-   // get the dimensions of the matrix stored in the file
-   // get dataspace handle
-   hid_t dataspace_id = H5Dget_space(dataset_id);
-   // get number of dimensions
-   const int ndims = H5Sget_simple_extent_ndims(dataspace_id);
-   AssertThrow(ndims==2, ExcIO());
-   // get every dimension
-   hsize_t dims[2];
-   status = H5Sget_simple_extent_dims(dataspace_id, dims, nullptr);
-   AssertThrow(status >= 0, ExcIO());
-   AssertThrow((int)dims[0]==n_columns,
-               ExcMessage("The number of columns of the matrix does not match the content of the archive"));
-   AssertThrow((int)dims[1]==n_rows,
-               ExcMessage("The number of rows of the matrix does not match the content of the archive"));
-
-   // gather the number of local rows and columns from all processes
-   std::vector<int> proc_n_local_rows(n_mpi_processes), proc_n_local_columns(n_mpi_processes);
-   MPI_Allgather(&tmp.n_local_rows,1,MPI_INT,proc_n_local_rows.data(),1,MPI_INT,tmp.grid->mpi_communicator);
-   MPI_Allgather(&tmp.n_local_columns,1,MPI_INT,proc_n_local_columns.data(),1,MPI_INT,tmp.grid->mpi_communicator);
-
-   const unsigned int my_rank(Utilities::MPI::this_mpi_process(tmp.grid->mpi_communicator));
-
-   // hyperslab selection parameters
-   // each process defines dataset in memory and writes it to the hyperslab in the file
-   hsize_t count[2];
-   count[0] = tmp.n_local_columns;
-   count[1] = tmp.n_local_rows;
-
-   hsize_t offset[2] = {0};
-   for (unsigned int i=0; i<my_rank; ++i)
-         offset[0] += proc_n_local_columns[i];
-
-   // select hyperslab in the file
-   status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, nullptr, count, nullptr);
-   AssertThrow(status >= 0, ExcIO());
-
-   // create a memory dataspace independently
-   hid_t memspace = H5Screate_simple(2, count, nullptr);
-
-   // read data independently
-   status = H5Dread(dataset_id, datatype, memspace, dataspace_id, H5P_DEFAULT, data);
-   AssertThrow(status >= 0, ExcIO());
-
-   // close/release sources
-   status = H5Sclose(memspace);
-   AssertThrow(status >= 0, ExcIO());
-   status = H5Dclose(dataset_id);
-   AssertThrow(status >= 0, ExcIO());
-   status = H5Sclose(dataspace_id);
-   AssertThrow(status >= 0, ExcIO());
-   status = H5Fclose(file_id);
-   AssertThrow(status >= 0, ExcIO());
-
-   // copying the distributed matrices
-   tmp.copy_to(*this);
+  // set up file access property list with parallel I/O access
+  hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
+  status = H5Pset_fapl_mpio(plist_id, tmp.grid->mpi_communicator, info);
+  AssertThrow(status >= 0, ExcIO());
+
+  // open file collectively in read-only mode and release property list identifier
+  hid_t file_id = H5Fopen(filename, H5F_ACC_RDONLY, plist_id);
+  status = H5Pclose(plist_id);
+  AssertThrow(status >= 0, ExcIO());
+
+  // open the dataset in the file collectively
+  hid_t dataset_id = H5Dopen2(file_id, "/matrix", H5P_DEFAULT);
+
+  // check the datatype of the dataset in the file
+  // if the classes of type of the dataset and the matrix do not match abort
+  // see HDF User's Guide: 6.10. Data Transfer: Datatype Conversion and Selection
+  hid_t datatype = hdf5_type_id(data);
+  hid_t datatype_inp  = H5Dget_type(dataset_id);
+  H5T_class_t t_class_inp = H5Tget_class(datatype_inp);
+  H5T_class_t t_class = H5Tget_class(datatype);
+  AssertThrow(t_class_inp == t_class,
+              ExcMessage("The data type of the matrix to be read does not match the archive"));
+
+  // get the dimensions of the matrix stored in the file
+  // get dataspace handle
+  hid_t dataspace_id = H5Dget_space(dataset_id);
+  // get number of dimensions
+  const int ndims = H5Sget_simple_extent_ndims(dataspace_id);
+  AssertThrow(ndims==2, ExcIO());
+  // get every dimension
+  hsize_t dims[2];
+  status = H5Sget_simple_extent_dims(dataspace_id, dims, nullptr);
+  AssertThrow(status >= 0, ExcIO());
+  AssertThrow((int)dims[0]==n_columns,
+              ExcMessage("The number of columns of the matrix does not match the content of the archive"));
+  AssertThrow((int)dims[1]==n_rows,
+              ExcMessage("The number of rows of the matrix does not match the content of the archive"));
+
+  // gather the number of local rows and columns from all processes
+  std::vector<int> proc_n_local_rows(n_mpi_processes), proc_n_local_columns(n_mpi_processes);
+  MPI_Allgather(&tmp.n_local_rows,1,MPI_INT,proc_n_local_rows.data(),1,MPI_INT,tmp.grid->mpi_communicator);
+  MPI_Allgather(&tmp.n_local_columns,1,MPI_INT,proc_n_local_columns.data(),1,MPI_INT,tmp.grid->mpi_communicator);
+
+  const unsigned int my_rank(Utilities::MPI::this_mpi_process(tmp.grid->mpi_communicator));
+
+  // hyperslab selection parameters
+  // each process defines dataset in memory and writes it to the hyperslab in the file
+  hsize_t count[2];
+  count[0] = tmp.n_local_columns;
+  count[1] = tmp.n_local_rows;
+
+  hsize_t offset[2] = {0};
+  for (unsigned int i=0; i<my_rank; ++i)
+    offset[0] += proc_n_local_columns[i];
+
+  // select hyperslab in the file
+  status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, nullptr, count, nullptr);
+  AssertThrow(status >= 0, ExcIO());
+
+  // create a memory dataspace independently
+  hid_t memspace = H5Screate_simple(2, count, nullptr);
+
+  // read data independently
+  status = H5Dread(dataset_id, datatype, memspace, dataspace_id, H5P_DEFAULT, data);
+  AssertThrow(status >= 0, ExcIO());
+
+  // close/release sources
+  status = H5Sclose(memspace);
+  AssertThrow(status >= 0, ExcIO());
+  status = H5Dclose(dataset_id);
+  AssertThrow(status >= 0, ExcIO());
+  status = H5Sclose(dataspace_id);
+  AssertThrow(status >= 0, ExcIO());
+  status = H5Fclose(file_id);
+  AssertThrow(status >= 0, ExcIO());
+
+  // copying the distributed matrices
+  tmp.copy_to(*this);
 
 #    endif // H5_HAVE_PARALLEL
 #  endif // DEAL_II_WITH_HDF5