]> https://gitweb.dealii.org/ - dealii.git/commitdiff
fix bug in computation of norms for ScaLAPACKMatrix 5868/head
authorBenjamin Brands <benjamin.brands@fau.de>
Tue, 6 Feb 2018 10:26:23 +0000 (11:26 +0100)
committerbenjamin <benjamin.brands@fau.de>
Fri, 9 Feb 2018 10:37:46 +0000 (11:37 +0100)
include/deal.II/lac/scalapack.h
include/deal.II/lac/scalapack.templates.h
source/lac/scalapack.cc
tests/scalapack/scalapack_04_b.cc [new file with mode: 0644]
tests/scalapack/scalapack_04_b.mpirun=1.output [new file with mode: 0644]
tests/scalapack/scalapack_04_b.mpirun=11.output [new file with mode: 0644]
tests/scalapack/scalapack_04_b.mpirun=4.output [new file with mode: 0644]

index 42468c81fc958e2c1b71e37253ee8be366564aac..27e5a2b8ee778b90ae7f52edad932fb10407062b 100644 (file)
@@ -391,11 +391,17 @@ public:
 
 private:
 
+  /**
+   * Calculate the norm of a distributed symmetric dense matrix using ScaLAPACK's
+   * internal function.
+   */
+  NumberType norm_symmetric(const char type) const;
+
   /**
    * Calculate the norm of a distributed dense matrix using ScaLAPACK's
    * internal function.
    */
-  NumberType norm(const char type) const;
+  NumberType norm_general(const char type) const;
 
   /**
    * Computing selected eigenvalues and, optionally, the eigenvectors.
index 6ca40dfbe08eb0ee112bb7049f54895d60916631..7479394c8e8e8505f00ba39298a20bafddf17c04 100644 (file)
@@ -403,20 +403,20 @@ extern "C"
    * or the element of largest absolute value of a distributed matrix
    */
   double pdlange_(char const *norm,
-                  int const &m,
-                  int const &n,
-                  double *A,
-                  int const &ia,
-                  int const &ja,
-                  int *desca,
+                  const int *m,
+                  const int *n,
+                  const double *A,
+                  const int *ia,
+                  const int *ja,
+                  const int *desca,
                   double *work);
-  float pslange_(char const *norm,
-                 int const &m,
-                 int const &n,
-                 float *A,
-                 int const &ia,
-                 int const &ja,
-                 int *desca,
+  float pslange_(const char *norm,
+                 const int *m,
+                 const int *n,
+                 const float *A,
+                 const int *ia,
+                 const int *ja,
+                 const int *desca,
                  float *work);
 
   /**
@@ -1043,37 +1043,37 @@ inline void pgemm(const char *transa,
 
 
 template <typename number>
-inline number plange(char const * /*norm*/,
-                     int const &/*m*/,
-                     int const &/*n*/,
-                     number * /*A*/,
-                     int const &/*ia*/,
-                     int const &/*ja*/,
-                     int * /*desca*/,
+inline number plange(const char * /*norm*/,
+                     const int * /*m*/,
+                     const int * /*n*/,
+                     const number * /*A*/,
+                     const int * /*ia*/,
+                     const int * /*ja*/,
+                     const int * /*desca*/,
                      number * /*work*/)
 {
   Assert (false, dealii::ExcNotImplemented());
 }
 
-inline double plange(char const *norm,
-                     int const &m,
-                     int const &n,
-                     double *A,
-                     int const &ia,
-                     int const &ja,
-                     int *desca,
+inline double plange(const char *norm,
+                     const int *m,
+                     const int *n,
+                     const double *A,
+                     const int *ia,
+                     const int *ja,
+                     const int *desca,
                      double *work)
 {
   return pdlange_(norm, m, n, A, ia, ja,desca, work);
 }
 
-inline float plange(char const *norm,
-                    int const &m,
-                    int const &n,
-                    float *A,
-                    int const &ia,
-                    int const &ja,
-                    int *desca,
+inline float plange(const char *norm,
+                    const int *m,
+                    const int *n,
+                    const float *A,
+                    const int *ia,
+                    const int *ja,
+                    const int *desca,
                     float *work)
 {
   return pslange_(norm, m, n, A, ia, ja,desca, work);
index 55d592b31a9e505d00a7156320de694d930dcecd..236123895bd211048c87221481286dc1a38146dc 100644 (file)
@@ -851,7 +851,11 @@ template <typename NumberType>
 NumberType ScaLAPACKMatrix<NumberType>::l1_norm() const
 {
   const char type('O');
-  return norm(type);
+
+  if (property == LAPACKSupport::symmetric)
+    return norm_symmetric(type);
+  else
+    return norm_general(type);
 }
 
 
@@ -860,7 +864,11 @@ template <typename NumberType>
 NumberType ScaLAPACKMatrix<NumberType>::linfty_norm() const
 {
   const char type('I');
-  return norm(type);
+
+  if (property == LAPACKSupport::symmetric)
+    return norm_symmetric(type);
+  else
+    return norm_general(type);
 }
 
 
@@ -869,17 +877,58 @@ template <typename NumberType>
 NumberType ScaLAPACKMatrix<NumberType>::frobenius_norm() const
 {
   const char type('F');
-  return norm(type);
+
+  if (property == LAPACKSupport::symmetric)
+    return norm_symmetric(type);
+  else
+    return norm_general(type);
+}
+
+
+
+template <typename NumberType>
+NumberType ScaLAPACKMatrix<NumberType>::norm_general(const char type) const
+{
+  Assert (state == LAPACKSupport::matrix ||
+          state == LAPACKSupport::inverse_matrix,
+          ExcMessage("norms can be called in matrix state only."));
+  Threads::Mutex::ScopedLock lock (mutex);
+  NumberType res = 0.;
+
+  if (grid->mpi_process_is_active)
+    {
+      const int iarow = indxg2p_(&submatrix_row, &row_block_size, &(grid->this_process_row), &first_process_row, &(grid->n_process_rows));
+      const int iacol = indxg2p_(&submatrix_column, &column_block_size, &(grid->this_process_column), &first_process_column, &(grid->n_process_columns));
+      const int mp0   = numroc_(&n_rows, &row_block_size, &(grid->this_process_row), &iarow, &(grid->n_process_rows));
+      const int nq0   = numroc_(&n_columns, &column_block_size, &(grid->this_process_column), &iacol, &(grid->n_process_columns));
+
+      // type='M': compute largest absolute value
+      // type='F' || type='E': compute Frobenius norm
+      // type='0' || type='1': compute infinity norm
+      int lwork=0; // for type == 'M' || type == 'F' || type == 'E'
+      if (type=='O' || type=='1')
+        lwork = nq0;
+      else if (type=='I')
+        lwork = mp0;
+
+      work.resize(lwork);
+      const NumberType *A_loc = this->values.begin();
+      res = plange(&type, &n_rows, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, work.data());
+    }
+  grid->send_to_inactive(&res);
+  return res;
 }
 
 
 
 template <typename NumberType>
-NumberType ScaLAPACKMatrix<NumberType>::norm(const char type) const
+NumberType ScaLAPACKMatrix<NumberType>::norm_symmetric(const char type) const
 {
   Assert (state == LAPACKSupport::matrix ||
           state == LAPACKSupport::inverse_matrix,
           ExcMessage("norms can be called in matrix state only."));
+  Assert (property == LAPACKSupport::symmetric,
+          ExcMessage("Matrix has to be symmetric for this operation."));
   Threads::Mutex::ScopedLock lock (mutex);
   NumberType res = 0.;
 
@@ -904,8 +953,8 @@ NumberType ScaLAPACKMatrix<NumberType>::norm(const char type) const
                         0 :
                         2*Nq0+Np0+ldw;
       work.resize(lwork);
-      const NumberType *A_loc = &this->values[0];
-      res = plansy(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, &work[0]);
+      const NumberType *A_loc = this->values.begin();
+      res = plansy(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, work.data());
     }
   grid->send_to_inactive(&res);
   return res;
diff --git a/tests/scalapack/scalapack_04_b.cc b/tests/scalapack/scalapack_04_b.cc
new file mode 100644 (file)
index 0000000..f472085
--- /dev/null
@@ -0,0 +1,100 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 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 norms of ScaLAPACK vs FullMatrix for general matrices
+
+#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/vector.h>
+
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+
+template <typename NumberType>
+void test(const unsigned int size, const unsigned int block_size, const NumberType tol)
+{
+  MPI_Comm mpi_communicator(MPI_COMM_WORLD);
+  const unsigned int n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator));
+  const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator));
+
+  ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
+
+  // Create SPD matrices of requested size:
+  FullMatrix<NumberType> full(size);
+
+  std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
+  ScaLAPACKMatrix<NumberType> scalapack (size, size, grid, block_size, block_size,
+                                         LAPACKSupport::Property::general);
+
+  pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
+
+  create_random (full);
+  scalapack = full;
+
+  const NumberType l1 = full.l1_norm();
+  const NumberType linfty = full.linfty_norm();
+  const NumberType frobenius = full.frobenius_norm();
+
+  // local result on this core:
+  const NumberType s_l1 = scalapack.l1_norm();
+  const NumberType s_linfty = scalapack.linfty_norm();
+  const NumberType s_frobenius = scalapack.frobenius_norm();
+
+  // make sure we have the same result on all cores, do average:
+  const NumberType as_l1 = dealii::Utilities::MPI::sum(s_l1, mpi_communicator)/n_mpi_processes;
+  const NumberType as_linfty = dealii::Utilities::MPI::sum(s_linfty, mpi_communicator)/n_mpi_processes;
+  const NumberType as_frobenius = dealii::Utilities::MPI::sum(s_frobenius, mpi_communicator)/n_mpi_processes;
+
+  pcout << l1 << " " << s_l1 << " " << as_l1 << std::endl
+        << linfty << " " << s_linfty << " " << as_linfty << std::endl
+        << frobenius << " " << s_frobenius << " " << as_frobenius << std::endl;
+
+  AssertThrow (std::abs(l1 -as_l1) < tol*std::abs(l1), dealii::ExcInternalError());
+  AssertThrow (std::abs(linfty -as_linfty) < tol*std::abs(linfty), dealii::ExcInternalError());
+  AssertThrow (std::abs(frobenius -as_frobenius) < tol*std::abs(frobenius), dealii::ExcInternalError());
+}
+
+
+
+int main (int argc,char **argv)
+{
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+  const std::vector<unsigned int> sizes = {{32,64,120,320,640}};
+  const std::vector<unsigned int> blocks = {{32,64}};
+
+  const double tol_double = 1e-10;
+  const float tol_float = 1e-5;
+
+  /*for (const auto &s : sizes)
+    for (const auto &b : blocks)
+      if (b <= s)
+        test<float>(s,b,tol_float);*/
+
+  for (const auto &s : sizes)
+    for (const auto &b : blocks)
+      if (b <= s)
+        test<double>(s,b,tol_double);
+}
diff --git a/tests/scalapack/scalapack_04_b.mpirun=1.output b/tests/scalapack/scalapack_04_b.mpirun=1.output
new file mode 100644 (file)
index 0000000..820f498
--- /dev/null
@@ -0,0 +1,36 @@
+32 32 1 1
+19.2741 19.2741 19.2741
+20.6796 20.6796 20.6796
+18.5785 18.5785 18.5785
+64 32 1 1
+38.5552 38.5552 38.5552
+38.2739 38.2739 38.2739
+36.8389 36.8389 36.8389
+64 64 1 1
+37.8671 37.8671 37.8671
+37.8657 37.8657 37.8657
+36.7171 36.7171 36.7171
+120 32 1 1
+69.0663 69.0663 69.0663
+68.0208 68.0208 68.0208
+69.1394 69.1394 69.1394
+120 64 1 1
+68.2078 68.2078 68.2078
+66.8542 66.8542 66.8542
+69.0112 69.0112 69.0112
+320 32 1 1
+174.387 174.387 174.387
+172.758 172.758 172.758
+184.876 184.876 184.876
+320 64 1 1
+182.948 182.948 182.948
+174.945 174.945 174.945
+184.521 184.521 184.521
+640 32 1 1
+338.935 338.935 338.935
+342.898 342.898 342.898
+369.65 369.65 369.65
+640 64 1 1
+342.421 342.421 342.421
+341.617 341.617 341.617
+369.497 369.497 369.497
diff --git a/tests/scalapack/scalapack_04_b.mpirun=11.output b/tests/scalapack/scalapack_04_b.mpirun=11.output
new file mode 100644 (file)
index 0000000..9f198a0
--- /dev/null
@@ -0,0 +1,36 @@
+32 32 1 1
+19.2741 19.2741 19.2741
+20.6796 20.6796 20.6796
+18.5785 18.5785 18.5785
+64 32 2 2
+38.5552 38.5552 38.5552
+38.2739 38.2739 38.2739
+36.8389 36.8389 36.8389
+64 64 1 1
+37.8671 37.8671 37.8671
+37.8657 37.8657 37.8657
+36.7171 36.7171 36.7171
+120 32 3 3
+69.0663 69.0663 69.0663
+68.0208 68.0208 68.0208
+69.1394 69.1394 69.1394
+120 64 2 2
+68.2078 68.2078 68.2078
+66.8542 66.8542 66.8542
+69.0112 69.0112 69.0112
+320 32 3 3
+174.387 174.387 174.387
+172.758 172.758 172.758
+184.876 184.876 184.876
+320 64 3 3
+182.948 182.948 182.948
+174.945 174.945 174.945
+184.521 184.521 184.521
+640 32 3 3
+338.935 338.935 338.935
+342.898 342.898 342.898
+369.65 369.65 369.65
+640 64 3 3
+342.421 342.421 342.421
+341.617 341.617 341.617
+369.497 369.497 369.497
diff --git a/tests/scalapack/scalapack_04_b.mpirun=4.output b/tests/scalapack/scalapack_04_b.mpirun=4.output
new file mode 100644 (file)
index 0000000..603d0a4
--- /dev/null
@@ -0,0 +1,36 @@
+32 32 1 1
+19.2741 19.2741 19.2741
+20.6796 20.6796 20.6796
+18.5785 18.5785 18.5785
+64 32 2 2
+38.5552 38.5552 38.5552
+38.2739 38.2739 38.2739
+36.8389 36.8389 36.8389
+64 64 1 1
+37.8671 37.8671 37.8671
+37.8657 37.8657 37.8657
+36.7171 36.7171 36.7171
+120 32 2 2
+69.0663 69.0663 69.0663
+68.0208 68.0208 68.0208
+69.1394 69.1394 69.1394
+120 64 2 2
+68.2078 68.2078 68.2078
+66.8542 66.8542 66.8542
+69.0112 69.0112 69.0112
+320 32 2 2
+174.387 174.387 174.387
+172.758 172.758 172.758
+184.876 184.876 184.876
+320 64 2 2
+182.948 182.948 182.948
+174.945 174.945 174.945
+184.521 184.521 184.521
+640 32 2 2
+338.935 338.935 338.935
+342.898 342.898 342.898
+369.65 369.65 369.65
+640 64 2 2
+342.421 342.421 342.421
+341.617 341.617 341.617
+369.497 369.497 369.497

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.