From: Daniel Arndt Date: Mon, 19 Feb 2024 22:10:07 +0000 (-0500) Subject: Tpetra: Add SparseMatrix::clear_row[s] X-Git-Tag: relicensing~14^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d9a2d3c5dedbc9a048f18a652f06c1c449d3f94;p=dealii.git Tpetra: Add SparseMatrix::clear_row[s] --- diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h index ad91ed457a..a32c520d51 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h @@ -677,6 +677,59 @@ namespace LinearAlgebra const Number *values, const bool elide_zero_values = false); + /** + * Remove all elements from this row by setting them to zero. The + * function does not modify the number of allocated nonzero entries, it + * only sets the entries to zero. + * + * This operation is used in eliminating constraints (e.g. due to hanging + * nodes) and makes sure that we can write this modification to the matrix + * without having to read entries (such as the locations of non-zero + * elements) from it — without this operation, removing constraints + * on %parallel matrices is a rather complicated procedure. + * + * The second parameter can be used to set the diagonal entry of this row + * to a value different from zero. The default is to set it to zero. + * + * @note If the matrix is stored in parallel across multiple processors + * using MPI, this function only touches rows that are locally stored and + * simply ignores all other row indices. Further, in the context of + * parallel computations, you will get into trouble if you clear a row + * while other processors still have pending writes or additions into the + * same row. In other words, if another processor still wants to add + * something to an element of a row and you call this function to zero out + * the row, then the next time you call compress() may add the remote + * value to the zero you just created. Consequently, you will want to call + * compress() after you made the last modifications to a matrix and before + * starting to clear rows. + */ + void + clear_row(const size_type row, const Number new_diag_value = 0); + + /** + * Same as clear_row(), except that it works on a number of rows at once. + * + * The second parameter can be used to set the diagonal entries of all + * cleared rows to something different from zero. Note that all of these + * diagonal entries get the same value -- if you want different values for + * the diagonal entries, you have to set them by hand. + * + * @note If the matrix is stored in parallel across multiple processors + * using MPI, this function only touches rows that are locally stored and + * simply ignores all other row indices. Further, in the context of + * parallel computations, you will get into trouble if you clear a row + * while other processors still have pending writes or additions into the + * same row. In other words, if another processor still wants to add + * something to an element of a row and you call this function to zero out + * the row, then the next time you call compress() may add the remote + * value to the zero you just created. Consequently, you will want to call + * compress() after you made the last modifications to a matrix and before + * starting to clear rows. + */ + void + clear_rows(const std::vector &rows, + const Number new_diag_value = 0); + /** * Release all memory and return to a state just like after having called * the default constructor. diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h index a8c1497817..2c3d9ab827 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h @@ -1005,6 +1005,68 @@ namespace LinearAlgebra + template + void + SparseMatrix::clear_row(const size_type row, + const Number new_diag_value) + { + clear_rows(std::vector(1, row), new_diag_value); + } + + + + template + void + SparseMatrix::clear_rows( + const std::vector &rows, + const Number new_diag_value) + { + // If the matrix is marked as compressed, we need to + // call resumeFill() first. + if (compressed || matrix->isFillComplete()) + { + matrix->resumeFill(); + compressed = false; + } + + std::vector col_indices_vector; + std::vector values_vector; + + for (size_type row : rows) + { + // Only do this on the rows owned locally on this processor. + int local_row = matrix->getRowMap()->getLocalElement(row); + if (local_row != Teuchos::OrdinalTraits::invalid()) + { + size_t nnz = matrix->getNumEntriesInLocalRow(local_row); + col_indices_vector.resize(nnz); + Teuchos::ArrayView col_indices(col_indices_vector); + values_vector.resize(nnz); + Teuchos::ArrayView values(values_vector); + + matrix->getLocalRowCopy(local_row, col_indices, values, nnz); + + const size_t diag_index = + std::find(col_indices.begin(), col_indices.end(), local_row) - + col_indices.begin(); + + for (size_t j = 0; j < nnz; ++j) + if (diag_index != j || new_diag_value == 0) + values[j] = 0.; + + if (diag_index != nnz) + values[diag_index] = new_diag_value; + + [[maybe_unused]] int n_replacements = + matrix->replaceLocalValues(local_row, col_indices, values); + AssertDimension(n_replacements, nnz); + } + } + compress(VectorOperation::insert); + } + + + template void SparseMatrix::copy_from( diff --git a/tests/trilinos_tpetra/62.cc b/tests/trilinos_tpetra/62.cc new file mode 100644 index 0000000000..68d316cdd8 --- /dev/null +++ b/tests/trilinos_tpetra/62.cc @@ -0,0 +1,89 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check LinearAlgebra::TpetraWrappers::SparseMatrix::clear () + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + AssertThrow(m.m() != 0, ExcInternalError()); + AssertThrow(m.n() != 0, ExcInternalError()); + + m.clear(); + + AssertThrow(m.m() == 0, ExcInternalError()); + AssertThrow(m.n() == 0, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix v(100U, 100U, 5U); + test(v); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/62.output b/tests/trilinos_tpetra/62.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/62.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos_tpetra/64.cc b/tests/trilinos_tpetra/64.cc new file mode 100644 index 0000000000..37ef14bfb5 --- /dev/null +++ b/tests/trilinos_tpetra/64.cc @@ -0,0 +1,113 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// This test should be run on multiple processors. + + +#include + +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +template +void +test(MatrixType &m) +{ + m.set(0, 0, 1.); + m.compress(VectorOperation::insert); + m = 0; + m.compress(VectorOperation::insert); + + Assert(fabs(m.frobenius_norm()) < 1e-15, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + const unsigned int n_dofs = 420; + // check + // LinearAlgebra::TpetraWrappers::SparseMatrix + LinearAlgebra::TpetraWrappers::SparseMatrix v1(n_dofs, + n_dofs, + 5U); + test(v1); + + // check + // LinearAlgebra::TpetraWrappers::SparseMatrix + const unsigned int n_jobs = + Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + const unsigned int my_id = + Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + Assert(n_dofs % n_jobs == 0, ExcInternalError()); + const unsigned int n_local_dofs = n_dofs / n_jobs; + IndexSet local_rows(n_dofs); + local_rows.add_range(n_local_dofs * my_id, n_local_dofs * (my_id + 1)); + LinearAlgebra::TpetraWrappers::SparseMatrix v2(local_rows, + MPI_COMM_WORLD, + 5); + test(v2); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/64.output b/tests/trilinos_tpetra/64.output new file mode 100644 index 0000000000..8b3b075900 --- /dev/null +++ b/tests/trilinos_tpetra/64.output @@ -0,0 +1,3 @@ + +DEAL::OK +DEAL::OK diff --git a/tests/trilinos_tpetra/66.cc b/tests/trilinos_tpetra/66.cc new file mode 100644 index 0000000000..4f6d653c09 --- /dev/null +++ b/tests/trilinos_tpetra/66.cc @@ -0,0 +1,141 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check LinearAlgebra::TpetraWrappers::MatrixBase::clear_row () + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + Assert(m.m() != 0, ExcInternalError()); + Assert(m.n() != 0, ExcInternalError()); + + // build a tri-diagonal pattern + double norm_sqr = 0; + unsigned int nnz = 0; + const unsigned int N = m.m(); + for (unsigned int i = 0; i < N; ++i) + { + if (i >= 5) + { + const double s = Testing::rand(); + m.set(i, i - 5, s); + norm_sqr += s * s; + ++nnz; + } + + if (i < N - 5) + { + const double s = Testing::rand(); + m.set(i, i + 5, s); + norm_sqr += s * s; + ++nnz; + } + + const double s = Testing::rand(); + m.set(i, i, s); + norm_sqr += s * s; + ++nnz; + } + m.compress(VectorOperation::insert); + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + Assert(m.n_nonzero_elements() - nnz == 0, ExcInternalError()); + + // now remove the entries of row N/2 + for (unsigned int i = 0; i < N; ++i) + { + const double s = m.el(N / 2, i); + norm_sqr -= s * s; + } + m.clear_row(N / 2); + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + + // make sure that zeroing out rows does at + // least not add new nonzero entries (it + // may remove some, though) + Assert(m.n_nonzero_elements() <= nnz, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix v(14U, 14U, 3U); + test(v); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/66.output b/tests/trilinos_tpetra/66.output new file mode 100644 index 0000000000..4572eca4cf --- /dev/null +++ b/tests/trilinos_tpetra/66.output @@ -0,0 +1,6 @@ + +DEAL::7.15267e+09 7.15267e+09 +DEAL::32 32 +DEAL::6.84337e+09 6.84337e+09 +DEAL::32 32 +DEAL::OK diff --git a/tests/trilinos_tpetra/67.cc b/tests/trilinos_tpetra/67.cc new file mode 100644 index 0000000000..16fa53189e --- /dev/null +++ b/tests/trilinos_tpetra/67.cc @@ -0,0 +1,147 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check LinearAlgebra::TpetraWrappers::MatrixBase::clear_rows () + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + Assert(m.m() != 0, ExcInternalError()); + Assert(m.n() != 0, ExcInternalError()); + + // build a tri-diagonal pattern + double norm_sqr = 0; + unsigned int nnz = 0; + const unsigned int N = m.m(); + for (unsigned int i = 0; i < N; ++i) + { + if (i >= 5) + { + const double s = Testing::rand(); + m.set(i, i - 5, s); + norm_sqr += s * s; + ++nnz; + } + + if (i < N - 5) + { + const double s = Testing::rand(); + m.set(i, i + 5, s); + norm_sqr += s * s; + ++nnz; + } + + const double s = Testing::rand(); + m.set(i, i, s); + norm_sqr += s * s; + ++nnz; + } + m.compress(VectorOperation::insert); + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + Assert(m.n_nonzero_elements() - nnz == 0, ExcInternalError()); + + // now remove the entries of rows N/2 and N/3 + for (unsigned int i = 0; i < N; ++i) + { + const double s = m.el(N / 2, i); + norm_sqr -= s * s; + } + for (unsigned int i = 0; i < N; ++i) + { + const double s = m.el(N / 3, i); + norm_sqr -= s * s; + } + const types::global_dof_index rows[2] = {N / 3, N / 2}; + m.clear_rows(std::vector(&rows[0], &rows[2])); + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + + // make sure that zeroing out rows does at + // least not add new nonzero entries (it + // may remove some, though) + AssertThrow(m.n_nonzero_elements() <= nnz, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix v(14U, 14U, 3U); + test(v); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/67.output b/tests/trilinos_tpetra/67.output new file mode 100644 index 0000000000..18be8791f8 --- /dev/null +++ b/tests/trilinos_tpetra/67.output @@ -0,0 +1,6 @@ + +DEAL::7.15267e+09 7.15267e+09 +DEAL::32 32 +DEAL::6.71272e+09 6.71272e+09 +DEAL::32 32 +DEAL::OK diff --git a/tests/trilinos_tpetra/68.cc b/tests/trilinos_tpetra/68.cc new file mode 100644 index 0000000000..7d09773906 --- /dev/null +++ b/tests/trilinos_tpetra/68.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check LinearAlgebra::TpetraWrappers::MatrixBase::clear_row () with used +// second argument + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + Assert(m.m() != 0, ExcInternalError()); + Assert(m.n() != 0, ExcInternalError()); + + // build a tri-diagonal pattern + double norm_sqr = 0; + unsigned int nnz = 0; + const unsigned int N = m.m(); + for (unsigned int i = 0; i < N; ++i) + { + if (i >= 5) + { + const double s = Testing::rand(); + m.set(i, i - 5, s); + norm_sqr += s * s; + ++nnz; + } + + if (i < N - 5) + { + const double s = Testing::rand(); + m.set(i, i + 5, s); + norm_sqr += s * s; + ++nnz; + } + + const double s = Testing::rand(); + m.set(i, i, s); + norm_sqr += s * s; + ++nnz; + } + m.compress(VectorOperation::insert); + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + Assert(m.n_nonzero_elements() - nnz == 0, ExcInternalError()); + + // now remove the entries of row N/2. set + // diagonal entries to rnd + const double rnd = Testing::rand(); + for (unsigned int i = 0; i < N; ++i) + { + const double s = m.el(N / 2, i); + norm_sqr -= s * s; + } + norm_sqr += rnd * rnd; + + m.clear_row(N / 2, rnd); + + Assert(m.el(N / 2, N / 2) == rnd, ExcInternalError()); + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + + // make sure that zeroing out rows does at + // least not add new nonzero entries (it + // may remove some, though) + Assert(m.n_nonzero_elements() <= nnz, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix v(14U, 14U, 3U); + test(v); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/68.output b/tests/trilinos_tpetra/68.output new file mode 100644 index 0000000000..bf77d2a368 --- /dev/null +++ b/tests/trilinos_tpetra/68.output @@ -0,0 +1,6 @@ + +DEAL::7.15267e+09 7.15267e+09 +DEAL::32 32 +DEAL::6.96869e+09 6.96869e+09 +DEAL::32 32 +DEAL::OK diff --git a/tests/trilinos_tpetra/69.cc b/tests/trilinos_tpetra/69.cc new file mode 100644 index 0000000000..0d74e9d3bb --- /dev/null +++ b/tests/trilinos_tpetra/69.cc @@ -0,0 +1,156 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check LinearAlgebra::TpetraWrappers::MatrixBase::clear_rows () with used +// second argument + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + Assert(m.m() != 0, ExcInternalError()); + Assert(m.n() != 0, ExcInternalError()); + + // build a tri-diagonal pattern + double norm_sqr = 0; + unsigned int nnz = 0; + const unsigned int N = m.m(); + for (unsigned int i = 0; i < N; ++i) + { + if (i >= 5) + { + const double s = Testing::rand(); + m.set(i, i - 5, s); + norm_sqr += s * s; + ++nnz; + } + + if (i < N - 5) + { + const double s = Testing::rand(); + m.set(i, i + 5, s); + norm_sqr += s * s; + ++nnz; + } + + const double s = Testing::rand(); + m.set(i, i, s); + norm_sqr += s * s; + ++nnz; + } + m.compress(VectorOperation::insert); + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + Assert(m.n_nonzero_elements() - nnz == 0, ExcInternalError()); + + // now remove the entries of rows N/2 and + // N/3. set diagonal entries to diag + const double diag = Testing::rand(); + for (unsigned int i = 0; i < N; ++i) + { + const double s = m.el(N / 2, i); + norm_sqr -= s * s; + } + for (unsigned int i = 0; i < N; ++i) + { + const double s = m.el(N / 3, i); + norm_sqr -= s * s; + } + norm_sqr += 2 * diag * diag; + + const std::vector rows = {N / 3, N / 2}; + m.clear_rows(rows, diag); + for (const auto row : rows) + { + Assert(m.el(row, row) == diag, ExcInternalError()); + } + + deallog << m.frobenius_norm() << ' ' << std::sqrt(norm_sqr) << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(std::fabs(m.frobenius_norm() - std::sqrt(norm_sqr)) < + std::fabs(std::sqrt(norm_sqr)), + ExcInternalError()); + + // make sure that zeroing out rows does at + // least not add new nonzero entries (it + // may remove some, though) + Assert(m.n_nonzero_elements() <= nnz, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix v(14U, 14U, 3U); + test(v); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/69.output b/tests/trilinos_tpetra/69.output new file mode 100644 index 0000000000..615e171a2c --- /dev/null +++ b/tests/trilinos_tpetra/69.output @@ -0,0 +1,6 @@ + +DEAL::7.15267e+09 7.15267e+09 +DEAL::32 32 +DEAL::6.96580e+09 6.96580e+09 +DEAL::32 32 +DEAL::OK