const Number *values,
const bool elide_zero_values = false);
+ /**
+ * Remove all elements from this <tt>row</tt> 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<size_type> &rows,
+ const Number new_diag_value = 0);
+
/**
* Release all memory and return to a state just like after having called
* the default constructor.
+ template <typename Number, typename MemorySpace>
+ void
+ SparseMatrix<Number, MemorySpace>::clear_row(const size_type row,
+ const Number new_diag_value)
+ {
+ clear_rows(std::vector<size_type>(1, row), new_diag_value);
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ void
+ SparseMatrix<Number, MemorySpace>::clear_rows(
+ const std::vector<size_type> &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<int> col_indices_vector;
+ std::vector<Number> 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<int>::invalid())
+ {
+ size_t nnz = matrix->getNumEntriesInLocalRow(local_row);
+ col_indices_vector.resize(nnz);
+ Teuchos::ArrayView<int> col_indices(col_indices_vector);
+ values_vector.resize(nnz);
+ Teuchos::ArrayView<Number> 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 <typename Number, typename MemorySpace>
void
SparseMatrix<Number, MemorySpace>::copy_from(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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<double>::clear ()
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &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<double> 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <Epetra_Comm.h>
+#include <Epetra_Map.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+template <typename MatrixType>
+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<double>
+ LinearAlgebra::TpetraWrappers::SparseMatrix<double> v1(n_dofs,
+ n_dofs,
+ 5U);
+ test(v1);
+
+ // check
+ // LinearAlgebra::TpetraWrappers::SparseMatrix<double>
+ 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<double> 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &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<double> 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;
+ };
+}
--- /dev/null
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.84337e+09 6.84337e+09
+DEAL::32 32
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &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<types::global_dof_index>(&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<double> 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;
+ };
+}
--- /dev/null
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.71272e+09 6.71272e+09
+DEAL::32 32
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &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<double> 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;
+ };
+}
--- /dev/null
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.96869e+09 6.96869e+09
+DEAL::32 32
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &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<types::global_dof_index> 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<double> 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;
+ };
+}
--- /dev/null
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.96580e+09 6.96580e+09
+DEAL::32 32
+DEAL::OK
projects / dealii.git / commitdiff