]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add tests for ScaLAPACKMatrix::pseudoinverse() for square and rectangular matrices 6119/head
authorBenjamin Brands <benjamin.brands@fau.de>
Fri, 6 Apr 2018 08:14:28 +0000 (10:14 +0200)
committerBenjamin Brands <benjamin.brands@fau.de>
Fri, 6 Apr 2018 08:39:13 +0000 (10:39 +0200)
tests/scalapack/scalapack_16.cc [new file with mode: 0644]
tests/scalapack/scalapack_16.mpirun=1.output [new file with mode: 0644]
tests/scalapack/scalapack_16.mpirun=11.output [new file with mode: 0644]
tests/scalapack/scalapack_16.mpirun=4.output [new file with mode: 0644]
tests/scalapack/scalapack_16_a.cc [new file with mode: 0644]
tests/scalapack/scalapack_16_a.mpirun=1.output [new file with mode: 0644]
tests/scalapack/scalapack_16_a.mpirun=11.output [new file with mode: 0644]
tests/scalapack/scalapack_16_a.mpirun=4.output [new file with mode: 0644]

diff --git a/tests/scalapack/scalapack_16.cc b/tests/scalapack/scalapack_16.cc
new file mode 100644 (file)
index 0000000..7bcf293
--- /dev/null
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+//
+// 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 pseudo_inverse(const NumberType ratio) for symmetric 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/base/process_grid.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/lapack_templates.h>
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include <memory>
+
+
+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));
+
+  std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
+
+  pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
+
+  // Create SPD matrices of requested size:
+  FullMatrix<NumberType> full_A(size), singular_A(size);
+  create_spd(full_A);
+  singular_A = full_A;
+  // Setting last row and column of the matrix to zero to make the matrix singular.
+  for (unsigned int i=0; i<size; ++i)
+    {
+      singular_A(i,size-1) = 0;
+      singular_A(size-1,i) = 0;
+    }
+  const NumberType ratio = 1e-8;
+
+  ScaLAPACKMatrix<NumberType> scalapack_A_1(size, grid, block_size);
+  scalapack_A_1 = singular_A;
+  scalapack_A_1.set_property(LAPACKSupport::Property::general);
+  const unsigned int rank_1 = scalapack_A_1.pseudoinverse(ratio);
+
+  ScaLAPACKMatrix<NumberType> inverse_A(size, grid, block_size);
+  inverse_A = full_A;
+  inverse_A.set_property(LAPACKSupport::Property::symmetric);
+  inverse_A.invert();
+
+  ScaLAPACKMatrix<NumberType> scalapack_A_2(size, grid, block_size);
+  scalapack_A_2 = full_A;
+  const unsigned int rank_2 = scalapack_A_2.pseudoinverse(ratio);
+
+  FullMatrix<NumberType> full_inverse(size);
+  inverse_A.copy_to(full_inverse);
+  FullMatrix<NumberType> full_pseudo_inverse(size);
+  scalapack_A_2.copy_to(full_pseudo_inverse);
+
+  if (this_mpi_process==0)
+    {
+      std::cout << "ranks: " << rank_1 << "/" << size << "  &  ";
+      std::cout << rank_2 << "/" << size << std::endl;
+      double norm = full_pseudo_inverse.frobenius_norm();
+      full_inverse.add(-1,full_pseudo_inverse);
+      norm = full_inverse.frobenius_norm() / norm;
+      AssertThrow(norm<tol,ExcInternalError());
+    }
+  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> sizes = {{200,400,600}};
+  const std::vector<unsigned int> blocks = {{32,64}};
+
+  const double tol = 1e-10;
+
+  for (const auto &s : sizes)
+    for (const auto &b : blocks)
+      if (b <= s)
+        test<double>(s,b,tol);
+}
diff --git a/tests/scalapack/scalapack_16.mpirun=1.output b/tests/scalapack/scalapack_16.mpirun=1.output
new file mode 100644 (file)
index 0000000..2c50041
--- /dev/null
@@ -0,0 +1,18 @@
+200 32 1 1
+ranks: 199/200  &  200/200
+
+200 64 1 1
+ranks: 199/200  &  200/200
+
+400 32 1 1
+ranks: 399/400  &  400/400
+
+400 64 1 1
+ranks: 399/400  &  400/400
+
+600 32 1 1
+ranks: 599/600  &  600/600
+
+600 64 1 1
+ranks: 599/600  &  600/600
+
diff --git a/tests/scalapack/scalapack_16.mpirun=11.output b/tests/scalapack/scalapack_16.mpirun=11.output
new file mode 100644 (file)
index 0000000..d58d105
--- /dev/null
@@ -0,0 +1,18 @@
+200 32 3 3
+ranks: 199/200  &  200/200
+
+200 64 3 3
+ranks: 199/200  &  200/200
+
+400 32 3 3
+ranks: 399/400  &  400/400
+
+400 64 3 3
+ranks: 399/400  &  400/400
+
+600 32 3 3
+ranks: 599/600  &  600/600
+
+600 64 3 3
+ranks: 599/600  &  600/600
+
diff --git a/tests/scalapack/scalapack_16.mpirun=4.output b/tests/scalapack/scalapack_16.mpirun=4.output
new file mode 100644 (file)
index 0000000..33bb38e
--- /dev/null
@@ -0,0 +1,18 @@
+200 32 2 2
+ranks: 199/200  &  200/200
+
+200 64 2 2
+ranks: 199/200  &  200/200
+
+400 32 2 2
+ranks: 399/400  &  400/400
+
+400 64 2 2
+ranks: 399/400  &  400/400
+
+600 32 2 2
+ranks: 599/600  &  600/600
+
+600 64 2 2
+ranks: 599/600  &  600/600
+
diff --git a/tests/scalapack/scalapack_16_a.cc b/tests/scalapack/scalapack_16_a.cc
new file mode 100644 (file)
index 0000000..cec06a2
--- /dev/null
@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// 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 pseudo_inverse(const NumberType ratio) for rectangular 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/base/process_grid.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/lapack_templates.h>
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include <memory>
+
+
+template <typename NumberType>
+void test(const unsigned int size_1, const unsigned int size_2, 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));
+
+  std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size_1,size_2,block_size,block_size);
+
+  pcout << size_1 << "x" << size_2 << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
+
+  // Create randomly filled matrix of requested dimension size_1xsize_2:
+  FullMatrix<NumberType> full_A(size_1,size_2);
+  create_random(full_A);
+  const NumberType ratio = 1e-8;
+
+  ScaLAPACKMatrix<NumberType> pseudoinverse_A(size_1,size_2,grid,block_size,block_size,LAPACKSupport::Property::general);
+  pseudoinverse_A = full_A;
+  const unsigned int rank = pseudoinverse_A.pseudoinverse(ratio);
+
+  // The pseudoinverse_A has to fulfill: A * A+ * A = A
+
+  ScaLAPACKMatrix<NumberType> scalapack_A(size_1,size_2,grid,block_size,block_size,LAPACKSupport::Property::general);
+  scalapack_A = full_A;
+  ScaLAPACKMatrix<NumberType> scalapack_tmp(size_2,size_2,grid,block_size,block_size,LAPACKSupport::Property::general);
+  // compute tmp = A+ * A
+  pseudoinverse_A.mmult(scalapack_tmp,scalapack_A);
+  ScaLAPACKMatrix<NumberType> scalapack_result(size_1,size_2,grid,block_size,block_size,LAPACKSupport::Property::general);
+  // compute result = A * tmp
+  scalapack_A.mmult(scalapack_result,scalapack_tmp);
+
+  // compute difference
+  scalapack_result.add(-1,scalapack_A);
+  const double norm = scalapack_result.frobenius_norm();
+
+  if (this_mpi_process==0)
+    {
+      std::cout << "rank=" << rank << "/" << std::min(size_1,size_2) << std::endl;
+      AssertThrow(norm<tol,ExcInternalError());
+    }
+  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> sizes_1 = {{200,400}};
+  const std::vector<unsigned int> sizes_2 = {{250,350}};
+  const std::vector<unsigned int> blocks = {{16,32}};
+
+  const double tol = 1e-10;
+
+  for (const auto &s_1 : sizes_1)
+    for (const auto &s_2 : sizes_2)
+      for (const auto &b : blocks)
+        if (b <= s_1 && b <= s_2)
+          test<double>(s_1,s_2,b,tol);
+}
diff --git a/tests/scalapack/scalapack_16_a.mpirun=1.output b/tests/scalapack/scalapack_16_a.mpirun=1.output
new file mode 100644 (file)
index 0000000..ed93062
--- /dev/null
@@ -0,0 +1,24 @@
+200x250 16 1 1
+rank=200/200
+
+200x250 32 1 1
+rank=200/200
+
+200x350 16 1 1
+rank=200/200
+
+200x350 32 1 1
+rank=200/200
+
+400x250 16 1 1
+rank=250/250
+
+400x250 32 1 1
+rank=250/250
+
+400x350 16 1 1
+rank=350/350
+
+400x350 32 1 1
+rank=350/350
+
diff --git a/tests/scalapack/scalapack_16_a.mpirun=11.output b/tests/scalapack/scalapack_16_a.mpirun=11.output
new file mode 100644 (file)
index 0000000..94c5f66
--- /dev/null
@@ -0,0 +1,24 @@
+200x250 16 3 3
+rank=200/200
+
+200x250 32 3 3
+rank=200/200
+
+200x350 16 2 4
+rank=200/200
+
+200x350 32 2 4
+rank=200/200
+
+400x250 16 5 2
+rank=250/250
+
+400x250 32 5 2
+rank=250/250
+
+400x350 16 3 3
+rank=350/350
+
+400x350 32 3 3
+rank=350/350
+
diff --git a/tests/scalapack/scalapack_16_a.mpirun=4.output b/tests/scalapack/scalapack_16_a.mpirun=4.output
new file mode 100644 (file)
index 0000000..c18a21d
--- /dev/null
@@ -0,0 +1,24 @@
+200x250 16 2 2
+rank=200/200
+
+200x250 32 2 2
+rank=200/200
+
+200x350 16 2 2
+rank=200/200
+
+200x350 32 2 2
+rank=200/200
+
+400x250 16 2 2
+rank=250/250
+
+400x250 32 2 2
+rank=250/250
+
+400x350 16 2 2
+rank=350/350
+
+400x350 32 2 2
+rank=350/350
+

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.