From: Benjamin Brands Date: Tue, 23 Jan 2018 09:37:01 +0000 (+0100) Subject: improve documentation and minor clean-ups X-Git-Tag: v9.0.0-rc1~483^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d93f320be90a7266affabd04ab8b81857f9a8a61;p=dealii.git improve documentation and minor clean-ups --- diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index 584e714b1c..5abfedd106 100644 --- 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=submatrix_size.first and number of columns=submatrix_size.second * . + * + * If it is necessary to copy complete matrices with an identical block-cyclic distribution, + * use copy_to(ScaLAPACKMatrix &dest) with only one argument to avoid communication */ void copy_to(ScaLAPACKMatrix &B, const std::pair &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 pXpotrf. The result of the factorization is stored in this object. diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index 6879cec96c..ffafdf5d76 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -23,9 +23,6 @@ #include #include #include -#include - -#include # ifdef DEAL_II_WITH_HDF5 #include @@ -935,7 +932,7 @@ void ScaLAPACKMatrix::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::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::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::save_parallel(const char *filename) const hsize_t offset[2] = {0}; for (unsigned int i=0; i::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::load_parallel(const char *filename) ScaLAPACKMatrix 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 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= 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 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= 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