From: Benjamin Brands Date: Fri, 6 Apr 2018 08:14:28 +0000 (+0200) Subject: add tests for ScaLAPACKMatrix::pseudoinverse() for square and rectangular matrices X-Git-Tag: v9.0.0-rc1~214^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F6119%2Fhead;p=dealii.git add tests for ScaLAPACKMatrix::pseudoinverse() for square and rectangular matrices --- diff --git a/tests/scalapack/scalapack_16.cc b/tests/scalapack/scalapack_16.cc new file mode 100644 index 0000000000..7bcf2934c5 --- /dev/null +++ b/tests/scalapack/scalapack_16.cc @@ -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 +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + + +template +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 grid = std::make_shared(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 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 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 inverse_A(size, grid, block_size); + inverse_A = full_A; + inverse_A.set_property(LAPACKSupport::Property::symmetric); + inverse_A.invert(); + + ScaLAPACKMatrix scalapack_A_2(size, grid, block_size); + scalapack_A_2 = full_A; + const unsigned int rank_2 = scalapack_A_2.pseudoinverse(ratio); + + FullMatrix full_inverse(size); + inverse_A.copy_to(full_inverse); + FullMatrix 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 sizes = {{200,400,600}}; + const std::vector blocks = {{32,64}}; + + const double tol = 1e-10; + + for (const auto &s : sizes) + for (const auto &b : blocks) + if (b <= s) + test(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 index 0000000000..2c50041b4d --- /dev/null +++ b/tests/scalapack/scalapack_16.mpirun=1.output @@ -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 index 0000000000..d58d105860 --- /dev/null +++ b/tests/scalapack/scalapack_16.mpirun=11.output @@ -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 index 0000000000..33bb38e95b --- /dev/null +++ b/tests/scalapack/scalapack_16.mpirun=4.output @@ -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 index 0000000000..cec06a21da --- /dev/null +++ b/tests/scalapack/scalapack_16_a.cc @@ -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 +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + + +template +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 grid = std::make_shared(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 full_A(size_1,size_2); + create_random(full_A); + const NumberType ratio = 1e-8; + + ScaLAPACKMatrix 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 scalapack_A(size_1,size_2,grid,block_size,block_size,LAPACKSupport::Property::general); + scalapack_A = full_A; + ScaLAPACKMatrix 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 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 sizes_1 = {{200,400}}; + const std::vector sizes_2 = {{250,350}}; + const std::vector 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(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 index 0000000000..ed93062d25 --- /dev/null +++ b/tests/scalapack/scalapack_16_a.mpirun=1.output @@ -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 index 0000000000..94c5f664dc --- /dev/null +++ b/tests/scalapack/scalapack_16_a.mpirun=11.output @@ -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 index 0000000000..c18a21d67d --- /dev/null +++ b/tests/scalapack/scalapack_16_a.mpirun=4.output @@ -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 +