]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Tpetra: Add SparseMatrix::clear_row[s]
authorDaniel Arndt <arndtd@ornl.gov>
Mon, 19 Feb 2024 22:10:07 +0000 (17:10 -0500)
committerDaniel Arndt <arndtd@ornl.gov>
Mon, 19 Feb 2024 22:10:07 +0000 (17:10 -0500)
14 files changed:
include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
tests/trilinos_tpetra/62.cc [new file with mode: 0644]
tests/trilinos_tpetra/62.output [new file with mode: 0644]
tests/trilinos_tpetra/64.cc [new file with mode: 0644]
tests/trilinos_tpetra/64.output [new file with mode: 0644]
tests/trilinos_tpetra/66.cc [new file with mode: 0644]
tests/trilinos_tpetra/66.output [new file with mode: 0644]
tests/trilinos_tpetra/67.cc [new file with mode: 0644]
tests/trilinos_tpetra/67.output [new file with mode: 0644]
tests/trilinos_tpetra/68.cc [new file with mode: 0644]
tests/trilinos_tpetra/68.output [new file with mode: 0644]
tests/trilinos_tpetra/69.cc [new file with mode: 0644]
tests/trilinos_tpetra/69.output [new file with mode: 0644]

index ad91ed457a65947698d4e6c4b84ec94309900497..a32c520d51d716327feffc5b526d0b1d3a8a6b83 100644 (file)
@@ -677,6 +677,59 @@ namespace LinearAlgebra
           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 &mdash; 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.
index a8c14978179dc610afc7809031441e84dd40b254..2c3d9ab8270d3cc3ec47b788e29aeecc84f95cc1 100644 (file)
@@ -1005,6 +1005,68 @@ namespace LinearAlgebra
 
 
 
+    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(
diff --git a/tests/trilinos_tpetra/62.cc b/tests/trilinos_tpetra/62.cc
new file mode 100644 (file)
index 0000000..68d316c
--- /dev/null
@@ -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<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;
+    };
+}
diff --git a/tests/trilinos_tpetra/62.output b/tests/trilinos_tpetra/62.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos_tpetra/64.cc b/tests/trilinos_tpetra/64.cc
new file mode 100644 (file)
index 0000000..37ef14b
--- /dev/null
@@ -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 <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;
+    };
+}
diff --git a/tests/trilinos_tpetra/64.output b/tests/trilinos_tpetra/64.output
new file mode 100644 (file)
index 0000000..8b3b075
--- /dev/null
@@ -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 (file)
index 0000000..4f6d653
--- /dev/null
@@ -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 <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;
+    };
+}
diff --git a/tests/trilinos_tpetra/66.output b/tests/trilinos_tpetra/66.output
new file mode 100644 (file)
index 0000000..4572eca
--- /dev/null
@@ -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 (file)
index 0000000..16fa531
--- /dev/null
@@ -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 <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;
+    };
+}
diff --git a/tests/trilinos_tpetra/67.output b/tests/trilinos_tpetra/67.output
new file mode 100644 (file)
index 0000000..18be879
--- /dev/null
@@ -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 (file)
index 0000000..7d09773
--- /dev/null
@@ -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 <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;
+    };
+}
diff --git a/tests/trilinos_tpetra/68.output b/tests/trilinos_tpetra/68.output
new file mode 100644 (file)
index 0000000..bf77d2a
--- /dev/null
@@ -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 (file)
index 0000000..0d74e9d
--- /dev/null
@@ -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 <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;
+    };
+}
diff --git a/tests/trilinos_tpetra/69.output b/tests/trilinos_tpetra/69.output
new file mode 100644 (file)
index 0000000..615e171
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.