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