From 8be5a8d9f8f6c45cff06e5ffa21bc82e7b110a87 Mon Sep 17 00:00:00 2001
From: Benjamin Brands <benjamin.brands@fau.de>
Date: Tue, 6 Feb 2018 16:24:53 +0100
Subject: [PATCH] routines to scale rows/columns of ScaLAPACKMatrix + tests

---
 include/deal.II/lac/scalapack.h               |  11 ++
 source/lac/scalapack.cc                       |  41 ++++++
 tests/scalapack/scalapack_14.cc               | 119 ++++++++++++++++++
 tests/scalapack/scalapack_14.mpirun=1.output  |  81 ++++++++++++
 tests/scalapack/scalapack_14.mpirun=11.output |  81 ++++++++++++
 tests/scalapack/scalapack_14.mpirun=9.output  |  81 ++++++++++++
 6 files changed, 414 insertions(+)
 create mode 100644 tests/scalapack/scalapack_14.cc
 create mode 100644 tests/scalapack/scalapack_14.mpirun=1.output
 create mode 100644 tests/scalapack/scalapack_14.mpirun=11.output
 create mode 100644 tests/scalapack/scalapack_14.mpirun=9.output

diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h
index 42468c81fc..609e481faf 100644
--- a/include/deal.II/lac/scalapack.h
+++ b/include/deal.II/lac/scalapack.h
@@ -26,6 +26,7 @@
 #include <deal.II/lac/lapack_support.h>
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/thread_management.h>
+#include <deal.II/base/array_view.h>
 #include <mpi.h>
 
 #include <memory>
@@ -389,6 +390,16 @@ public:
    */
   NumberType &local_el(const unsigned int loc_row, const unsigned int loc_column);
 
+  /**
+   * scaling the columns of the distributed matrix by scalars in array @p factors
+   */
+  void scale_columns(const ArrayView<const NumberType> &factors);
+
+  /**
+   * scaling the rows of the distributed matrix by scalars in array @p factors
+   */
+  void scale_rows(const ArrayView<const NumberType> &factors);
+
 private:
 
   /**
diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc
index 55d592b31a..aafc4efc8a 100644
--- a/source/lac/scalapack.cc
+++ b/source/lac/scalapack.cc
@@ -22,6 +22,7 @@
 
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/mpi.templates.h>
+#include <deal.II/base/array_view.h>
 #include <deal.II/lac/scalapack.templates.h>
 
 #ifdef DEAL_II_WITH_HDF5
@@ -1320,6 +1321,46 @@ void ScaLAPACKMatrix<NumberType>::load_parallel(const char *filename)
 
 
 
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::scale_columns(const ArrayView<const NumberType> &factors)
+{
+  Assert(n_columns==(int)factors.size(),ExcDimensionMismatch(n_columns,factors.size()));
+
+  if (this->grid->mpi_process_is_active)
+    {
+      for (int i=0; i<n_local_rows; ++i)
+        {
+          for (int j=0; j<n_local_columns; ++j)
+            {
+              const int glob_j = global_column(j);
+              local_el(i,j) *= factors[glob_j];
+            }
+        }
+    }
+}
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::scale_rows(const ArrayView<const NumberType> &factors)
+{
+  Assert(n_rows==(int)factors.size(),ExcDimensionMismatch(n_rows,factors.size()));
+
+  if (this->grid->mpi_process_is_active)
+    {
+      for (int i=0; i<n_local_rows; ++i)
+        {
+          const int glob_i = global_row(i);
+          for (int j=0; j<n_local_columns; ++j)
+            {
+              local_el(i,j) *= factors[glob_i];
+            }
+        }
+    }
+}
+
+
+
 // instantiations
 template class ScaLAPACKMatrix<double>;
 template class ScaLAPACKMatrix<float>;
diff --git a/tests/scalapack/scalapack_14.cc b/tests/scalapack/scalapack_14.cc
new file mode 100644
index 0000000000..d5501b2986
--- /dev/null
+++ b/tests/scalapack/scalapack_14.cc
@@ -0,0 +1,119 @@
+// ---------------------------------------------------------------------
+//
+// 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 scaling of rows and columns 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/base/array_view.h>
+
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+#include <typeinfo>
+
+
+template <typename NumberType>
+void test(const unsigned int block_size_i, const unsigned int block_size_j)
+{
+  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));
+
+  std::cout << std::setprecision(10);
+  ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
+
+  const unsigned int proc_rows = std::floor(std::sqrt(n_mpi_processes));
+  const unsigned int proc_columns = std::floor(n_mpi_processes/proc_rows);
+  //create 2d process grid
+  const std::vector<unsigned int> sizes = {{400,500}};
+  std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,sizes[0],sizes[1],block_size_i,block_size_i);
+  pcout << "2D process grid: " << grid->get_process_grid_rows() << "x" << grid->get_process_grid_columns() << std::endl << std::endl;
+
+  // test scaling of rows
+  {
+    FullMatrix<NumberType> full_A(sizes[0],sizes[1]);
+    create_random(full_A);
+
+    std::vector<NumberType> scaling_factors(full_A.m());
+    for (unsigned int i=0; i<scaling_factors.size(); ++i)
+      scaling_factors[i] = std::sqrt(i+1);
+
+    ScaLAPACKMatrix<NumberType> scalapack_A (full_A.m(),full_A.n(),grid,block_size_i,block_size_j);
+    scalapack_A = full_A;
+    const ArrayView<NumberType> view_rows(scaling_factors);
+    scalapack_A.scale_rows(view_rows);
+    FullMatrix<NumberType> tmp_full_A(scalapack_A.m(),scalapack_A.n());
+    scalapack_A.copy_to(tmp_full_A);
+
+    for (unsigned int i=0; i<full_A.m(); ++i)
+      for (unsigned int j=0; j<full_A.n(); ++j)
+        full_A(i,j) *= scaling_factors[i];
+
+    pcout << "   Row scaling for"
+          << " A in R^(" << scalapack_A.m() << "x" << scalapack_A.n() << ")" << std::endl;
+    pcout << "   norms: " << tmp_full_A.frobenius_norm()<< " & "
+          << full_A.frobenius_norm() << "  for "
+          << typeid(NumberType).name() << std::endl << std::endl;
+  }
+  // test scaling of columns
+  {
+    FullMatrix<NumberType> full_A(sizes[0],sizes[1]);
+    create_random(full_A);
+
+    std::vector<NumberType> scaling_factors(full_A.n());
+    for (unsigned int i=0; i<scaling_factors.size(); ++i)
+      scaling_factors[i] = std::sqrt(i+1);
+
+    ScaLAPACKMatrix<NumberType> scalapack_A (full_A.m(),full_A.n(),grid,block_size_i,block_size_j);
+    scalapack_A = full_A;
+    const ArrayView<NumberType> view_columns(scaling_factors);
+    scalapack_A.scale_columns(view_columns);
+    FullMatrix<NumberType> tmp_full_A(scalapack_A.m(),scalapack_A.n());
+    scalapack_A.copy_to(tmp_full_A);
+
+    for (unsigned int i=0; i<full_A.m(); ++i)
+      for (unsigned int j=0; j<full_A.n(); ++j)
+        full_A(i,j) *= scaling_factors[j];
+
+    pcout << "   Column scaling for"
+          << " A in R^(" << scalapack_A.m() << "x" << scalapack_A.n() << ")" << std::endl;
+    pcout << "   norms: " << tmp_full_A.frobenius_norm()<< " & "
+          << full_A.frobenius_norm() << "  for "
+          << typeid(NumberType).name() << std::endl << std::endl;
+  }
+  pcout << std::endl;
+}
+
+
+
+int main (int argc,char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+  const std::vector<unsigned int> blocks_i = {{16,32,64}};
+  const std::vector<unsigned int> blocks_j = {{16,32,64}};
+
+  for (const auto &s : blocks_i)
+    for (const auto &b : blocks_j)
+      test<double>(s,b);
+}
diff --git a/tests/scalapack/scalapack_14.mpirun=1.output b/tests/scalapack/scalapack_14.mpirun=1.output
new file mode 100644
index 0000000000..4dfac23a67
--- /dev/null
+++ b/tests/scalapack/scalapack_14.mpirun=1.output
@@ -0,0 +1,81 @@
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.285795 & 3654.285795  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4090.694598 & 4090.694598  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.738726 & 3654.738726  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4085.273405 & 4085.273405  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3660.42714 & 3660.42714  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4082.735765 & 4082.735765  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.388891 & 3654.388891  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4084.61645 & 4084.61645  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3655.261484 & 3655.261484  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4084.19574 & 4084.19574  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3655.547554 & 3655.547554  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4088.843423 & 4088.843423  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.096322 & 3654.096322  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4083.807486 & 4083.807486  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3662.0357 & 3662.0357  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4089.328538 & 4089.328538  for d
+
+
+2D process grid: 1x1
+
+   Row scaling for A in R^(400x500)
+   norms: 3663.288472 & 3663.288472  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4081.172517 & 4081.172517  for d
+
+
diff --git a/tests/scalapack/scalapack_14.mpirun=11.output b/tests/scalapack/scalapack_14.mpirun=11.output
new file mode 100644
index 0000000000..58e9b076e7
--- /dev/null
+++ b/tests/scalapack/scalapack_14.mpirun=11.output
@@ -0,0 +1,81 @@
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.285795 & 3654.285795  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4090.694598 & 4090.694598  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.738726 & 3654.738726  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4085.273405 & 4085.273405  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3660.42714 & 3660.42714  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4082.735765 & 4082.735765  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.388891 & 3654.388891  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4084.61645 & 4084.61645  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3655.261484 & 3655.261484  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4084.19574 & 4084.19574  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3655.547554 & 3655.547554  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4088.843423 & 4088.843423  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.096322 & 3654.096322  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4083.807486 & 4083.807486  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3662.0357 & 3662.0357  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4089.328538 & 4089.328538  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3663.288472 & 3663.288472  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4081.172517 & 4081.172517  for d
+
+
diff --git a/tests/scalapack/scalapack_14.mpirun=9.output b/tests/scalapack/scalapack_14.mpirun=9.output
new file mode 100644
index 0000000000..58e9b076e7
--- /dev/null
+++ b/tests/scalapack/scalapack_14.mpirun=9.output
@@ -0,0 +1,81 @@
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.285795 & 3654.285795  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4090.694598 & 4090.694598  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.738726 & 3654.738726  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4085.273405 & 4085.273405  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3660.42714 & 3660.42714  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4082.735765 & 4082.735765  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.388891 & 3654.388891  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4084.61645 & 4084.61645  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3655.261484 & 3655.261484  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4084.19574 & 4084.19574  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3655.547554 & 3655.547554  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4088.843423 & 4088.843423  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3654.096322 & 3654.096322  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4083.807486 & 4083.807486  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3662.0357 & 3662.0357  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4089.328538 & 4089.328538  for d
+
+
+2D process grid: 3x3
+
+   Row scaling for A in R^(400x500)
+   norms: 3663.288472 & 3663.288472  for d
+
+   Column scaling for A in R^(400x500)
+   norms: 4081.172517 & 4081.172517  for d
+
+
-- 
2.39.5