From 93f3c810548e7f58b30715d833f74ab70dce0a8f Mon Sep 17 00:00:00 2001
From: Benjamin Brands <benjamin.brands@fau.de>
Date: Tue, 6 Feb 2018 11:26:23 +0100
Subject: [PATCH] fix bug in computation of norms for ScaLAPACKMatrix

---
 include/deal.II/lac/scalapack.h               |   8 +-
 include/deal.II/lac/scalapack.templates.h     |  68 ++++++------
 source/lac/scalapack.cc                       |  61 +++++++++--
 tests/scalapack/scalapack_04_b.cc             | 100 ++++++++++++++++++
 .../scalapack/scalapack_04_b.mpirun=1.output  |  36 +++++++
 .../scalapack/scalapack_04_b.mpirun=11.output |  36 +++++++
 .../scalapack/scalapack_04_b.mpirun=4.output  |  36 +++++++
 7 files changed, 304 insertions(+), 41 deletions(-)
 create mode 100644 tests/scalapack/scalapack_04_b.cc
 create mode 100644 tests/scalapack/scalapack_04_b.mpirun=1.output
 create mode 100644 tests/scalapack/scalapack_04_b.mpirun=11.output
 create mode 100644 tests/scalapack/scalapack_04_b.mpirun=4.output

diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h
index 42468c81fc..27e5a2b8ee 100644
--- a/include/deal.II/lac/scalapack.h
+++ b/include/deal.II/lac/scalapack.h
@@ -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.
diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h
index 6ca40dfbe0..7479394c8e 100644
--- a/include/deal.II/lac/scalapack.templates.h
+++ b/include/deal.II/lac/scalapack.templates.h
@@ -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);
diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc
index 55d592b31a..236123895b 100644
--- a/source/lac/scalapack.cc
+++ b/source/lac/scalapack.cc
@@ -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
index 0000000000..f472085f24
--- /dev/null
+++ b/tests/scalapack/scalapack_04_b.cc
@@ -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
index 0000000000..820f498747
--- /dev/null
+++ b/tests/scalapack/scalapack_04_b.mpirun=1.output
@@ -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
index 0000000000..9f198a0d7d
--- /dev/null
+++ b/tests/scalapack/scalapack_04_b.mpirun=11.output
@@ -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
index 0000000000..603d0a4cb7
--- /dev/null
+++ b/tests/scalapack/scalapack_04_b.mpirun=4.output
@@ -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
-- 
2.39.5