--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+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
+
--- /dev/null
+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
+
--- /dev/null
+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
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+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
+
--- /dev/null
+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
+
--- /dev/null
+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
+
projects / dealii.git / commitdiff