save/load state and property of ScaLAPACKMatrix

diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h
index 08a75f2d78ec3a28465d14644268e4ad65fbfe91..5cff9a44d5d7faa7f500dd93b37c3517699c85d4 100644 (file)
--- a/include/deal.II/lac/scalapack.h
+++ b/include/deal.II/lac/scalapack.h
@@ -141,6 +141,16 @@ public:
    */
   void set_property(const LAPACKSupport::Property property);
 
+  /**
+   * Return current @p property of this matrix
+   */
+  LAPACKSupport::Property get_property() const;
+
+  /**
+   * Return current @p state of this matrix
+   */
+  LAPACKSupport::State get_state() const;
+
   /**
    * Assignment operator from a regular FullMatrix.
    *

diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc
index 7bd19f121e38d7484ee216c2dba3abe1aa64d614..3d40a53cf2b629116acf8ba311e36d57c8446e41 100644 (file)
--- a/source/lac/scalapack.cc
+++ b/source/lac/scalapack.cc
@@ -155,6 +155,24 @@ ScaLAPACKMatrix<NumberType>::set_property(const LAPACKSupport::Property property
 
 
 
+template <typename NumberType>
+LAPACKSupport::Property
+ScaLAPACKMatrix<NumberType>::get_property() const
+{
+  return property;
+}
+
+
+
+template <typename NumberType>
+LAPACKSupport::State
+ScaLAPACKMatrix<NumberType>::get_state() const
+{
+  return state;
+}
+
+
+
 template <typename NumberType>
 ScaLAPACKMatrix<NumberType> &
 ScaLAPACKMatrix<NumberType>::operator = (const FullMatrix<NumberType> &matrix)
@@ -1158,6 +1176,69 @@ NumberType ScaLAPACKMatrix<NumberType>::norm_symmetric(const char type) const
 
 
 
+namespace internal
+{
+  namespace
+  {
+    void create_HDF5_state_enum_id(hid_t &state_enum_id)
+    {
+      // create HDF5 enum type for LAPACKSupport::State
+      LAPACKSupport::State val;
+      state_enum_id = H5Tcreate (H5T_ENUM, sizeof(LAPACKSupport::State));
+      val = LAPACKSupport::State::cholesky;
+      herr_t status = H5Tenum_insert (state_enum_id, "cholesky", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+      val = LAPACKSupport::State::eigenvalues;
+      status = H5Tenum_insert (state_enum_id, "eigenvalues", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+      val = LAPACKSupport::State::inverse_matrix;
+      status = H5Tenum_insert (state_enum_id, "inverse_matrix", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+      val = LAPACKSupport::State::inverse_svd;
+      status = H5Tenum_insert (state_enum_id, "inverse_svd", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+      val = LAPACKSupport::State::lu;
+      status = H5Tenum_insert (state_enum_id, "lu", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+      val = LAPACKSupport::State::matrix;
+      status = H5Tenum_insert (state_enum_id, "matrix", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+      val = LAPACKSupport::State::svd;
+      status = H5Tenum_insert (state_enum_id, "svd", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+      val = LAPACKSupport::State::unusable;
+      status = H5Tenum_insert (state_enum_id, "unusable", (int *)&val);
+      AssertThrow(status >= 0, ExcInternalError());
+    }
+
+    void create_HDF5_property_enum_id(hid_t &property_enum_id)
+    {
+      // create HDF5 enum type for LAPACKSupport::Property
+      property_enum_id = H5Tcreate (H5T_ENUM, sizeof(LAPACKSupport::Property));
+      LAPACKSupport::Property prop = LAPACKSupport::Property::diagonal;
+      herr_t status = H5Tenum_insert (property_enum_id, "diagonal", (int *)&prop);
+      AssertThrow(status >= 0, ExcInternalError());
+      prop = LAPACKSupport::Property::general;
+      status = H5Tenum_insert (property_enum_id, "general", (int *)&prop);
+      AssertThrow(status >= 0, ExcInternalError());
+      prop = LAPACKSupport::Property::hessenberg;
+      status = H5Tenum_insert (property_enum_id, "hessenberg", (int *)&prop);
+      AssertThrow(status >= 0, ExcInternalError());
+      prop = LAPACKSupport::Property::lower_triangular;
+      status = H5Tenum_insert (property_enum_id, "lower_triangular", (int *)&prop);
+      AssertThrow(status >= 0, ExcInternalError());
+      prop = LAPACKSupport::Property::symmetric;
+      status = H5Tenum_insert (property_enum_id, "symmetric", (int *)&prop);
+      AssertThrow(status >= 0, ExcInternalError());
+      prop = LAPACKSupport::Property::upper_triangular;
+      status = H5Tenum_insert (property_enum_id, "upper_triangular", (int *)&prop);
+      AssertThrow(status >= 0, ExcInternalError());
+    }
+  }
+}
+
+
+
 template <typename NumberType>
 void ScaLAPACKMatrix<NumberType>::save(const char *filename,
                                        const std::pair<unsigned int,unsigned int> &chunk_size) const
@@ -1230,8 +1311,8 @@ void ScaLAPACKMatrix<NumberType>::save_serial(const char *filename,
       //revert order of rows and columns as ScaLAPACK uses column-major ordering
       chunk_dims[0] = chunk_size.second;
       chunk_dims[1] = chunk_size.first;
-      hid_t property = H5Pcreate (H5P_DATASET_CREATE);
-      status = H5Pset_chunk (property, 2, chunk_dims);
+      hid_t data_property = H5Pcreate (H5P_DATASET_CREATE);
+      status = H5Pset_chunk (data_property, 2, chunk_dims);
       AssertThrow(status >= 0, ExcIO());
 
       // create the data space for the dataset
@@ -1245,7 +1326,7 @@ void ScaLAPACKMatrix<NumberType>::save_serial(const char *filename,
       hid_t type_id = hdf5_type_id(&tmp.values[0]);
       hid_t dataset_id = H5Dcreate2(file_id, "/matrix",
                                     type_id, dataspace_id,
-                                    H5P_DEFAULT, property, H5P_DEFAULT);
+                                    H5P_DEFAULT, data_property, H5P_DEFAULT);
 
       // write the dataset
       status = H5Dwrite(dataset_id, type_id,
@@ -1253,16 +1334,61 @@ void ScaLAPACKMatrix<NumberType>::save_serial(const char *filename,
                         &tmp.values[0]);
       AssertThrow(status >= 0, ExcIO());
 
-      // end access to the dataset and release resources used by it
+      // create HDF5 enum type for LAPACKSupport::State and LAPACKSupport::Property
+      hid_t state_enum_id, property_enum_id;
+      internal::create_HDF5_state_enum_id(state_enum_id);
+      internal::create_HDF5_property_enum_id(property_enum_id);
+
+      // create the data space for the state enum
+      hsize_t dims_state[1];
+      dims_state[0]=1;
+      hid_t state_enum_dataspace = H5Screate_simple(1, dims_state, nullptr);
+      // create the dataset for the state enum
+      hid_t state_enum_dataset = H5Dcreate2(file_id, "/state", state_enum_id, state_enum_dataspace,
+                                            H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+      // write the dataset for the state enum
+      status = H5Dwrite(state_enum_dataset, state_enum_id,
+                        H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                        &state);
+      AssertThrow(status >= 0, ExcIO());
+
+      // create the data space for the property enum
+      hsize_t dims_property[1];
+      dims_property[0]=1;
+      hid_t property_enum_dataspace = H5Screate_simple(1, dims_property, nullptr);
+      // create the dataset for the property enum
+      hid_t property_enum_dataset = H5Dcreate2(file_id, "/property", property_enum_id, property_enum_dataspace,
+                                               H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+      // write the dataset for the property enum
+      status = H5Dwrite(property_enum_dataset, property_enum_id,
+                        H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                        &property);
+      AssertThrow(status >= 0, ExcIO());
+
+      // end access to the datasets and release resources used by them
       status = H5Dclose(dataset_id);
       AssertThrow(status >= 0, ExcIO());
+      status = H5Dclose(state_enum_dataset);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Dclose(property_enum_dataset);
+      AssertThrow(status >= 0, ExcIO());
 
-      // terminate access to the data space
+      // terminate access to the data spaces
       status = H5Sclose(dataspace_id);
       AssertThrow(status >= 0, ExcIO());
+      status = H5Sclose(state_enum_dataspace);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Sclose(property_enum_dataspace);
+      AssertThrow(status >= 0, ExcIO());
+
+      // release enum data types
+      status = H5Tclose(state_enum_id);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Tclose(property_enum_id);
+      AssertThrow(status >= 0, ExcIO());
 
       // release the creation property
-      status = H5Pclose (property);
+      status = H5Pclose (data_property);
       AssertThrow(status >= 0, ExcIO());
 
       // close the file.
@@ -1334,8 +1460,10 @@ void ScaLAPACKMatrix<NumberType>::save_parallel(const char *filename,
   hid_t type_id = hdf5_type_id(data);
   hid_t dset_id = H5Dcreate2(file_id, "/matrix", type_id,
                              filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT);
+
   status = H5Sclose(filespace);
   AssertThrow(status >= 0, ExcIO());
+
   status = H5Pclose(plist_id);
   AssertThrow(status >= 0, ExcIO());
 
@@ -1374,7 +1502,6 @@ void ScaLAPACKMatrix<NumberType>::save_parallel(const char *filename,
                         plist_id, data);
       AssertThrow(status >= 0, ExcIO());
     }
-
   // close/release sources
   status = H5Dclose(dset_id);
   AssertThrow(status >= 0, ExcIO());
@@ -1387,6 +1514,63 @@ void ScaLAPACKMatrix<NumberType>::save_parallel(const char *filename,
   status = H5Fclose(file_id);
   AssertThrow(status >= 0, ExcIO());
 
+  // before writing the state and property to file wait for
+  // all processes to finish writing the matrix content to the file
+  MPI_Barrier(tmp.grid->mpi_communicator);
+
+  // only root process will write state and property to the file
+  if (tmp.grid->this_mpi_process==0)
+    {
+      // open file using default properties
+      hid_t file_id_reopen = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
+
+      // create HDF5 enum type for LAPACKSupport::State and LAPACKSupport::Property
+      hid_t state_enum_id, property_enum_id;
+      internal::create_HDF5_state_enum_id(state_enum_id);
+      internal::create_HDF5_property_enum_id(property_enum_id);
+
+      // create the data space for the state enum
+      hsize_t dims_state[1];
+      dims_state[0]=1;
+      hid_t state_enum_dataspace = H5Screate_simple(1, dims_state, nullptr);
+      // create the dataset for the state enum
+      hid_t state_enum_dataset = H5Dcreate2(file_id_reopen, "/state", state_enum_id, state_enum_dataspace,
+                                            H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+      // write the dataset for the state enum
+      status = H5Dwrite(state_enum_dataset, state_enum_id,
+                        H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                        &state);
+      AssertThrow(status >= 0, ExcIO());
+
+      // create the data space for the property enum
+      hsize_t dims_property[1];
+      dims_property[0]=1;
+      hid_t property_enum_dataspace = H5Screate_simple(1, dims_property, nullptr);
+      // create the dataset for the property enum
+      hid_t property_enum_dataset = H5Dcreate2(file_id_reopen, "/property", property_enum_id, property_enum_dataspace,
+                                               H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+      // write the dataset for the property enum
+      status = H5Dwrite(property_enum_dataset, property_enum_id,
+                        H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                        &property);
+      AssertThrow(status >= 0, ExcIO());
+
+      status = H5Dclose(state_enum_dataset);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Dclose(property_enum_dataset);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Sclose(state_enum_dataspace);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Sclose(property_enum_dataspace);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Tclose(state_enum_id);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Tclose(property_enum_id);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Fclose(file_id_reopen);
+      AssertThrow(status >= 0, ExcIO());
+    }
+
 #  endif
 }
 
@@ -1402,6 +1586,7 @@ void ScaLAPACKMatrix<NumberType>::load(const char *filename)
 #  ifdef H5_HAVE_PARALLEL
   //implementation for configurations equipped with a parallel file system
   load_parallel(filename);
+
 #  else
   //implementation for configurations with no parallel file system
   load_serial(filename);
@@ -1429,6 +1614,9 @@ void ScaLAPACKMatrix<NumberType>::load_serial(const char *filename)
   const int MB=n_rows, NB=n_columns;
   ScaLAPACKMatrix<NumberType> tmp(n_rows,n_columns,one_grid,MB,NB);
 
+  int state_int = -1;
+  int property_int = -1;
+
   // the 1x1 grid has only one process and this one reads
   // the content of the matrix from the HDF5 file
   if (tmp.grid->mpi_process_is_active)
@@ -1468,22 +1656,91 @@ void ScaLAPACKMatrix<NumberType>::load_serial(const char *filename)
                        H5P_DEFAULT, &tmp.values[0]);
       AssertThrow(status >= 0, ExcIO());
 
-      // terminate access to the data space
+      // create HDF5 enum type for LAPACKSupport::State and LAPACKSupport::Property
+      hid_t state_enum_id, property_enum_id;
+      internal::create_HDF5_state_enum_id(state_enum_id);
+      internal::create_HDF5_property_enum_id(property_enum_id);
+
+      // open the datasets for the state and property enum in the file
+      hid_t dataset_state_id = H5Dopen2(file_id, "/state", H5P_DEFAULT);
+      hid_t datatype_state  = H5Dget_type(dataset_state_id);
+      H5T_class_t t_class_state = H5Tget_class(datatype_state);
+      AssertThrow(t_class_state == H5T_ENUM, ExcIO());
+
+      hid_t dataset_property_id = H5Dopen2(file_id, "/property", H5P_DEFAULT);
+      hid_t datatype_property  = H5Dget_type(dataset_property_id);
+      H5T_class_t t_class_property = H5Tget_class(datatype_property);
+      AssertThrow(t_class_property == H5T_ENUM, ExcIO());
+
+      // get dataspace handles
+      hid_t dataspace_state = H5Dget_space(dataset_state_id);
+      hid_t dataspace_property = H5Dget_space(dataset_property_id);
+      // get number of dimensions
+      const int ndims_state = H5Sget_simple_extent_ndims(dataspace_state);
+      AssertThrow(ndims_state==1, ExcIO());
+      const int ndims_property = H5Sget_simple_extent_ndims(dataspace_property);
+      AssertThrow(ndims_property==1, ExcIO());
+      // get every dimension
+      hsize_t dims_state[1];
+      H5Sget_simple_extent_dims(dataspace_state, dims_state, nullptr);
+      AssertThrow((int)dims_state[0]==1,ExcIO());
+      hsize_t dims_property[1];
+      H5Sget_simple_extent_dims(dataspace_property, dims_property, nullptr);
+      AssertThrow((int)dims_property[0]==1,ExcIO());
+
+      // read data
+      status = H5Dread(dataset_state_id, state_enum_id, H5S_ALL, H5S_ALL,
+                       H5P_DEFAULT, &tmp.state);
+      AssertThrow(status >= 0, ExcIO());
+      // To send the state from the root process to the other processes
+      // the state enum is casted to an integer, that will be broadcasted and
+      // subsequently casted back to the enum type
+      state_int = static_cast<int>(tmp.state);
+
+      status = H5Dread(dataset_property_id, property_enum_id, H5S_ALL, H5S_ALL,
+                       H5P_DEFAULT, &tmp.property);
+      AssertThrow(status >= 0, ExcIO());
+      // To send the property from the root process to the other processes
+      // the state enum is casted to an integer, that will be broadcasted and
+      // subsequently casted back to the enum type
+      property_int = static_cast<int>(tmp.property);
+
+      // terminate access to the data spaces
       status = H5Sclose(dataspace_id);
       AssertThrow(status >= 0, ExcIO());
+      status = H5Sclose(dataspace_state);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Sclose(dataspace_property);
+      AssertThrow(status >= 0, ExcIO());
 
-      // release data type handle
+      // release data type handles
       status = H5Tclose(datatype);
       AssertThrow(status >= 0, ExcIO());
+      status = H5Tclose(state_enum_id);
+      AssertThrow(status >= 0, ExcIO());
+      status = H5Tclose(property_enum_id);
+      AssertThrow(status >= 0, ExcIO());
 
-      // end access to the dataset and release resources used by it
+      // end access to the data sets and release resources used by them
+      status = H5Dclose(dataset_state_id);
+      AssertThrow(status >= 0, ExcIO());
       status = H5Dclose(dataset_id);
       AssertThrow(status >= 0, ExcIO());
+      status = H5Dclose(dataset_property_id);
+      AssertThrow(status >= 0, ExcIO());
 
       // close the file.
       status = H5Fclose(file_id);
       AssertThrow(status >= 0, ExcIO());
     }
+  // so far only the root process has the correct state integer --> broadcasting
+  tmp.grid->send_to_inactive(&state_int,1);
+  // so far only the root process has the correct property integer --> broadcasting
+  tmp.grid->send_to_inactive(&property_int,1);
+
+  tmp.state = static_cast<LAPACKSupport::State>(state_int);
+  tmp.property = static_cast<LAPACKSupport::Property>(property_int);
+
   tmp.copy_to(*this);
 
 #  endif // DEAL_II_WITH_HDF5
@@ -1585,13 +1842,68 @@ void ScaLAPACKMatrix<NumberType>::load_parallel(const char *filename)
   status = H5Dread(dataset_id, datatype, memspace, dataspace_id, H5P_DEFAULT, data);
   AssertThrow(status >= 0, ExcIO());
 
+  // create HDF5 enum type for LAPACKSupport::State and LAPACKSupport::Property
+  hid_t state_enum_id, property_enum_id;
+  internal::create_HDF5_state_enum_id(state_enum_id);
+  internal::create_HDF5_property_enum_id(property_enum_id);
+
+  // open the datasets for the state and property enum in the file
+  hid_t dataset_state_id = H5Dopen2(file_id, "/state", H5P_DEFAULT);
+  hid_t datatype_state  = H5Dget_type(dataset_state_id);
+  H5T_class_t t_class_state = H5Tget_class(datatype_state);
+  AssertThrow(t_class_state == H5T_ENUM, ExcIO());
+
+  hid_t dataset_property_id = H5Dopen2(file_id, "/property", H5P_DEFAULT);
+  hid_t datatype_property  = H5Dget_type(dataset_property_id);
+  H5T_class_t t_class_property = H5Tget_class(datatype_property);
+  AssertThrow(t_class_property == H5T_ENUM, ExcIO());
+
+  // get dataspace handles
+  hid_t dataspace_state = H5Dget_space(dataset_state_id);
+  hid_t dataspace_property = H5Dget_space(dataset_property_id);
+  // get number of dimensions
+  const int ndims_state = H5Sget_simple_extent_ndims(dataspace_state);
+  AssertThrow(ndims_state==1, ExcIO());
+  const int ndims_property = H5Sget_simple_extent_ndims(dataspace_property);
+  AssertThrow(ndims_property==1, ExcIO());
+  // get every dimension
+  hsize_t dims_state[1];
+  H5Sget_simple_extent_dims(dataspace_state, dims_state, nullptr);
+  AssertThrow((int)dims_state[0]==1,ExcIO());
+  hsize_t dims_property[1];
+  H5Sget_simple_extent_dims(dataspace_property, dims_property, nullptr);
+  AssertThrow((int)dims_property[0]==1,ExcIO());
+
+  // read data
+  status = H5Dread(dataset_state_id, state_enum_id, H5S_ALL, H5S_ALL,
+                   H5P_DEFAULT, &tmp.state);
+  AssertThrow(status >= 0, ExcIO());
+
+  status = H5Dread(dataset_property_id, property_enum_id, H5S_ALL, H5S_ALL,
+                   H5P_DEFAULT, &tmp.property);
+  AssertThrow(status >= 0, ExcIO());
+
   // close/release sources
   status = H5Sclose(memspace);
   AssertThrow(status >= 0, ExcIO());
   status = H5Dclose(dataset_id);
   AssertThrow(status >= 0, ExcIO());
+  status = H5Dclose(dataset_state_id);
+  AssertThrow(status >= 0, ExcIO());
+  status = H5Dclose(dataset_property_id);
+  AssertThrow(status >= 0, ExcIO());
   status = H5Sclose(dataspace_id);
   AssertThrow(status >= 0, ExcIO());
+  status = H5Sclose(dataspace_state);
+  AssertThrow(status >= 0, ExcIO());
+  status = H5Sclose(dataspace_property);
+  AssertThrow(status >= 0, ExcIO());
+  //status = H5Tclose(datatype);
+  //AssertThrow(status >= 0, ExcIO());
+  status = H5Tclose(state_enum_id);
+  AssertThrow(status >= 0, ExcIO());
+  status = H5Tclose(property_enum_id);
+  AssertThrow(status >= 0, ExcIO());
   status = H5Fclose(file_id);
   AssertThrow(status >= 0, ExcIO());
 

diff --git a/tests/scalapack/scalapack_10.cc b/tests/scalapack/scalapack_10.cc
index 02f0cfd94a5b2bfbdc594679a64d5244355a32d5..962fa7cad028fd376ecf50d49f6276ee2ca58d4d 100644 (file)
--- a/tests/scalapack/scalapack_10.cc
+++ b/tests/scalapack/scalapack_10.cc
@@ -16,7 +16,7 @@
 #include "../tests.h"
 #include "../lapack/create_matrix.h"
 
-// test serial saving and loading of distributed ScaLAPACKMatrices
+// test saving and loading of distributed ScaLAPACKMatrices
 
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/utilities.h>
@@ -34,7 +34,7 @@
 template <typename NumberType>
 void test(const unsigned int size, const unsigned int block_size)
 {
-  const std::string filename ("scalapck_10_test.h5");
+  const std::string filename ("scalapack_10_test.h5");
 
   MPI_Comm mpi_communicator(MPI_COMM_WORLD);
   const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator));
@@ -62,6 +62,10 @@ void test(const unsigned int size, const unsigned int block_size)
   copy.add(-1,full);
 
   pcout << size << " " << block_size << std::endl;
+
+  if (copy.frobenius_norm()>1e-12)
+    pcout << "norm of difference: " << copy.frobenius_norm() << std::endl;
+
   AssertThrow(copy.frobenius_norm()<1e-12,ExcInternalError());
   std::remove(filename.c_str());
 }

diff --git a/tests/scalapack/scalapack_10_b.cc b/tests/scalapack/scalapack_10_b.cc
index 3e09b515bb3bf1b8f8c061b8f64ce71b49271bf6..a56c9a4704ca39c2bfaf69263f6f20ae86d5d948 100644 (file)
--- a/tests/scalapack/scalapack_10_b.cc
+++ b/tests/scalapack/scalapack_10_b.cc
@@ -16,7 +16,7 @@
 #include "../tests.h"
 #include "../lapack/create_matrix.h"
 
-// test serial saving and loading of distributed ScaLAPACKMatrices with prescribed chunk sizes
+// test saving and loading of distributed ScaLAPACKMatrices with prescribed chunk sizes
 
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/utilities.h>
@@ -34,7 +34,7 @@
 template <typename NumberType>
 void test(const std::pair<unsigned int,unsigned int> &size, const unsigned int block_size, const std::pair<unsigned int,unsigned int> &chunk_size)
 {
-  const std::string filename ("scalapck_10_b_test.h5");
+  const std::string filename ("scalapack_10_b_test.h5");
 
   MPI_Comm mpi_communicator(MPI_COMM_WORLD);
   const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator));

diff --git a/tests/scalapack/scalapack_10_c.cc b/tests/scalapack/scalapack_10_c.cc
new file mode 100644 (file)
index 0000000..e6bf6f2
--- /dev/null
+++ b/tests/scalapack/scalapack_10_c.cc
@@ -0,0 +1,131 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include "../tests.h"
+#include "../lapack/create_matrix.h"
+
+// test saving and loading of the State and Property of distributed ScaLAPACKMatrices
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/multithread_info.h>
+
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+#include <cstdio>
+
+
+template <typename NumberType>
+void test()
+{
+  const std::string filename ("scalapack_10_test.h5");
+
+  MPI_Comm mpi_communicator(MPI_COMM_WORLD);
+  const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator));
+  ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
+
+  pcout << "Saving and restoring the state and property of ScaLAPACKMatrix" << std::endl;
+
+  const unsigned int size=100, block_size=8;
+
+  //create FullMatrix and fill it
+  FullMatrix<NumberType> full(100);
+  create_spd (full);
+
+  //create 2d process grid
+  std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
+
+  ScaLAPACKMatrix<NumberType> scalapack_matrix(size,size,grid,block_size,block_size);
+  ScaLAPACKMatrix<NumberType> scalapack_matrix_copy(size,size,grid,block_size,block_size);
+
+  scalapack_matrix.set_property(LAPACKSupport::Property::diagonal);
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_property()==scalapack_matrix_copy.get_property(),
+              ExcInternalError());
+
+  scalapack_matrix.set_property(LAPACKSupport::Property::general);
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_property()==scalapack_matrix_copy.get_property(),
+              ExcInternalError());
+
+  scalapack_matrix.set_property(LAPACKSupport::Property::hessenberg);
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_property()==scalapack_matrix_copy.get_property(),
+              ExcInternalError());
+
+  scalapack_matrix.set_property(LAPACKSupport::Property::lower_triangular);
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_property()==scalapack_matrix_copy.get_property(),
+              ExcInternalError());
+
+  scalapack_matrix.set_property(LAPACKSupport::Property::symmetric);
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_property()==scalapack_matrix_copy.get_property(),
+              ExcInternalError());
+
+  scalapack_matrix.set_property(LAPACKSupport::Property::upper_triangular);
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_property()==scalapack_matrix_copy.get_property(),
+              ExcInternalError());
+
+  // after construction the matrix state is LAPACKSupport::State::unusable
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_state()==scalapack_matrix_copy.get_state(),
+              ExcInternalError());
+
+  // the assignment operator changes the state to LAPACKSupport::State::matrix
+  scalapack_matrix = full;
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_state()==scalapack_matrix_copy.get_state(),
+              ExcInternalError());
+
+  // calling invert changes the state to LAPACKSupport::inverse_matrix
+  scalapack_matrix.set_property(LAPACKSupport::Property::symmetric);
+  scalapack_matrix.invert();
+  scalapack_matrix.save(filename.c_str());
+  scalapack_matrix_copy.load(filename.c_str());
+  std::remove(filename.c_str());
+  AssertThrow(scalapack_matrix.get_state()==scalapack_matrix_copy.get_state(),
+              ExcInternalError());
+}
+
+
+
+int main (int argc,char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+  test<double>();
+}

diff --git a/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=1.output b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=1.output
new file mode 100644 (file)
index 0000000..b65f4d9
--- /dev/null
+++ b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=1.output
@@ -0,0 +1 @@
+Saving and restoring the state and property of ScaLAPACKMatrix

diff --git a/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=11.output b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=11.output
new file mode 100644 (file)
index 0000000..b65f4d9
--- /dev/null
+++ b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=11.output
@@ -0,0 +1 @@
+Saving and restoring the state and property of ScaLAPACKMatrix

diff --git a/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=4.output b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=4.output
new file mode 100644 (file)
index 0000000..b65f4d9
--- /dev/null
+++ b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=4.output
@@ -0,0 +1 @@
+Saving and restoring the state and property of ScaLAPACKMatrix

diff --git a/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=5.output b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=5.output
new file mode 100644 (file)
index 0000000..b65f4d9
--- /dev/null
+++ b/tests/scalapack/scalapack_10_c.with_hdf5=on.mpirun=5.output
@@ -0,0 +1 @@
+Saving and restoring the state and property of ScaLAPACKMatrix