]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Tpetra: Implement serial constructors, operator(), el(), and diag_element() for Spars... 16616/head
authorDaniel Arndt <arndtd@ornl.gov>
Fri, 9 Feb 2024 04:39:34 +0000 (23:39 -0500)
committerDaniel Arndt <arndtd@ornl.gov>
Fri, 9 Feb 2024 20:46:59 +0000 (15:46 -0500)
15 files changed:
include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
tests/trilinos_tpetra/01.cc [new file with mode: 0644]
tests/trilinos_tpetra/01.output [new file with mode: 0644]
tests/trilinos_tpetra/02.cc [new file with mode: 0644]
tests/trilinos_tpetra/02.output [new file with mode: 0644]
tests/trilinos_tpetra/03a.cc [new file with mode: 0644]
tests/trilinos_tpetra/03a.output [new file with mode: 0644]
tests/trilinos_tpetra/04.cc [new file with mode: 0644]
tests/trilinos_tpetra/04.output [new file with mode: 0644]
tests/trilinos_tpetra/05.cc [new file with mode: 0644]
tests/trilinos_tpetra/05.output [new file with mode: 0644]
tests/trilinos_tpetra/09.cc [new file with mode: 0644]
tests/trilinos_tpetra/09.output [new file with mode: 0644]
tests/trilinos_tpetra/CMakeLists.txt [new file with mode: 0644]

index 772202b2e279e54411ea5e3ba9aeeff27ea4acd6..b145e7a22023140eb5e33da84a335ff8ec10f1cf 100644 (file)
@@ -152,6 +152,28 @@ namespace LinearAlgebra
        */
       SparseMatrix(const SparsityPattern<MemorySpace> &sparsity_pattern);
 
+      /**
+       * Generate a matrix that is completely stored locally, having #m rows and
+       * #n columns.
+       *
+       * The number of columns entries per row is specified as the maximum
+       * number of entries argument.
+       */
+      SparseMatrix(const size_type    m,
+                   const size_type    n,
+                   const unsigned int n_max_entries_per_row);
+
+      /**
+       * Generate a matrix that is completely stored locally, having #m rows and
+       * #n columns.
+       *
+       * The vector <tt>n_entries_per_row</tt> specifies the number of entries
+       * in each row.
+       */
+      SparseMatrix(const size_type                  m,
+                   const size_type                  n,
+                   const std::vector<unsigned int> &n_entries_per_row);
+
       /**
        * Move constructor. Create a new sparse matrix by stealing the internal
        * data of the `other` object.
@@ -587,7 +609,50 @@ namespace LinearAlgebra
           const bool       elide_zero_values = false);
 
       /** @} */
+      /**
+       * @name Entry Access
+       */
+      /** @{ */
 
+      /**
+       * Return the value of the entry (<i>i,j</i>).  This may be an expensive
+       * operation and you should always take care where to call this function.
+       * As in the deal.II sparse matrix class, we throw an exception if the
+       * respective entry doesn't exist in the sparsity pattern of this class,
+       * which is requested from Trilinos. Moreover, an exception will be thrown
+       * when the requested element is not saved on the calling process.
+       */
+      Number
+      operator()(const size_type i, const size_type j) const;
+
+      /**
+       * Return the value of the matrix entry (<i>i,j</i>). If this entry does
+       * not exist in the sparsity pattern, then zero is returned. While this
+       * may be convenient in some cases, note that it is simple to write
+       * algorithms that are slow compared to an optimal solution, since the
+       * sparsity of the matrix is not used.  On the other hand, if you want to
+       * be sure the entry exists, you should use operator() instead.
+       *
+       * The lack of error checking in this function can also yield surprising
+       * results if you have a parallel matrix. In that case, just because you
+       * get a zero result from this function does not mean that either the
+       * entry does not exist in the sparsity pattern or that it does but has a
+       * value of zero. Rather, it could also be that it simply isn't stored on
+       * the current processor; in that case, it may be stored on a different
+       * processor, and possibly so with a nonzero value.
+       */
+      Number
+      el(const size_type i, const size_type j) const;
+
+      /**
+       * Return the main diagonal element in the <i>i</i>th row. This function
+       * throws an error if the matrix is not quadratic and it also throws an
+       * error if <i>(i,i)</i> is not element of the local matrix.
+       */
+      Number
+      diag_element(const size_type i) const;
+
+      /** @} */
       /**
        * @name Multiplications
        */
@@ -802,9 +867,39 @@ namespace LinearAlgebra
                        "Are you trying to put the product of the "
                        "matrix with a vector into a vector that has "
                        "ghost elements?");
+
+      /**
+       * Exception
+       */
+      DeclException2(ExcInvalidIndex,
+                     size_type,
+                     size_type,
+                     << "The entry with index <" << arg1 << ',' << arg2
+                     << "> does not exist.");
+
+      /**
+       * Exception
+       */
+      DeclException4(ExcAccessToNonLocalElement,
+                     size_type,
+                     size_type,
+                     size_type,
+                     size_type,
+                     << "You tried to access element (" << arg1 << '/' << arg2
+                     << ')'
+                     << " of a distributed matrix, but only rows in range ["
+                     << arg3 << ',' << arg4
+                     << "] are stored locally and can be accessed.");
+
       /** @} */
 
     private:
+      /**
+       * Helper function for operator() and el().
+       */
+      Number
+      element(const size_type i, const size_type j, const bool no_error) const;
+
       /**
        * Pointer to the user-supplied Tpetra Trilinos mapping of the matrix
        * columns that assigns parts of the matrix to the individual processes.
index 6e9ae879b54af81f21312a245fe5511328364a42..8273703b263f6120584d8b78f0f3605e9e8f6394 100644 (file)
@@ -230,6 +230,49 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    SparseMatrix<Number, MemorySpace>::SparseMatrix(
+      const size_type    m,
+      const size_type    n,
+      const unsigned int n_max_entries_per_row)
+      : column_space_map(Utilities::Trilinos::internal::make_rcp<MapType>(
+          n,
+          0,
+          Utilities::Trilinos::tpetra_comm_self()))
+      , matrix(Utilities::Trilinos::internal::make_rcp<MatrixType>(
+          Utilities::Trilinos::internal::make_rcp<MapType>(
+            m,
+            0,
+            Utilities::Trilinos::tpetra_comm_self()),
+          column_space_map,
+          n_max_entries_per_row))
+      , compressed(false)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    SparseMatrix<Number, MemorySpace>::SparseMatrix(
+      const size_type                  m,
+      const size_type                  n,
+      const std::vector<unsigned int> &n_entries_per_row)
+      : column_space_map(Utilities::Trilinos::internal::make_rcp<MapType>(
+          n,
+          0,
+          Utilities::Trilinos::tpetra_comm_self()))
+      , compressed(false)
+    {
+      std::vector<size_t> entries_per_row_size_type(n_entries_per_row.begin(),
+                                                    n_entries_per_row.end());
+      matrix = Utilities::Trilinos::internal::make_rcp<MatrixType>(
+        Utilities::Trilinos::internal::make_rcp<MapType>(
+          m, 0, Utilities::Trilinos::tpetra_comm_self()),
+        column_space_map,
+        Teuchos::ArrayView<size_t>{entries_per_row_size_type});
+    }
+
+
+
     template <typename Number, typename MemorySpace>
     SparseMatrix<Number, MemorySpace>::SparseMatrix(
       SparseMatrix<Number, MemorySpace> &&other) noexcept
@@ -826,6 +869,8 @@ namespace LinearAlgebra
         }
     }
 
+
+
     template <typename Number, typename MemorySpace>
     void
     SparseMatrix<Number, MemorySpace>::resume_fill()
@@ -837,6 +882,101 @@ namespace LinearAlgebra
         }
     }
 
+
+    template <typename Number, typename MemorySpace>
+    Number
+    SparseMatrix<Number, MemorySpace>::element(const size_type i,
+                                               const size_type j,
+                                               const bool      no_error) const
+    {
+      // Extract local indices in the matrix.
+      const int      trilinos_i = matrix->getRowMap()->getLocalElement(i);
+      const int      trilinos_j = matrix->getColMap()->getLocalElement(j);
+      TrilinosScalar value      = 0.;
+
+      if (trilinos_i == Teuchos::OrdinalTraits<int>::invalid() ||
+          trilinos_j == Teuchos::OrdinalTraits<int>::invalid())
+        {
+          if (no_error)
+            return {};
+          Assert(false,
+                 ExcAccessToNonLocalElement(
+                   i, j, local_range().first, local_range().second - 1));
+        }
+      else
+        {
+          Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
+
+          // Prepare pointers for extraction of a view of the row.
+          size_t nnz_present = matrix->getNumEntriesInLocalRow(trilinos_i);
+          typename MatrixType::nonconst_local_inds_host_view_type col_indices(
+            "indices", nnz_present);
+          typename MatrixType::nonconst_values_host_view_type values(
+            "values", nnz_present);
+
+          matrix->getLocalRowCopy(trilinos_i, col_indices, values, nnz_present);
+
+          // Search the index where we look for the value, and then finally get
+          // it.
+          int local_col_index = 0;
+          for (; local_col_index < static_cast<int>(nnz_present);
+               ++local_col_index)
+            {
+              if (col_indices(local_col_index) == trilinos_j)
+                break;
+            }
+
+          if (local_col_index == static_cast<int>(nnz_present))
+            {
+              if (no_error)
+                return {};
+              Assert(false, ExcInvalidIndex(i, j));
+            }
+          else
+            value = values(local_col_index);
+        }
+
+      return value;
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    Number
+    SparseMatrix<Number, MemorySpace>::operator()(const size_type i,
+                                                  const size_type j) const
+    {
+      return element(i, j, /* no_error */ false);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    Number
+    SparseMatrix<Number, MemorySpace>::el(const size_type i,
+                                          const size_type j) const
+    {
+      return element(i, j, /* no_error */ true);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    Number
+    SparseMatrix<Number, MemorySpace>::diag_element(const size_type i) const
+    {
+      Assert(m() == n(), ExcNotQuadratic());
+
+#  ifdef DEBUG
+      // use operator() in debug mode because it checks if this is a valid
+      // element (in parallel)
+      return operator()(i, i);
+#  else
+      // Trilinos doesn't seem to have a more efficient way to access the
+      // diagonal than by just using the standard el(i,j) function.
+      return el(i, i);
+#  endif
+    }
   } // namespace TpetraWrappers
 
 } // namespace LinearAlgebra
diff --git a/tests/trilinos_tpetra/01.cc b/tests/trilinos_tpetra/01.cc
new file mode 100644 (file)
index 0000000..e39faab
--- /dev/null
@@ -0,0 +1,100 @@
+// ---------------------------------------------------------------------
+//
+// 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 setting elements in a trilinos matrix using
+// LinearAlgebra::TpetraWrappers::SparseMatrix<double>::set()
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &m)
+{
+  // first set a few entries
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        m.set(i, j, i * j * .5 + .5);
+
+  m.compress(VectorOperation::insert);
+
+  // then make sure we retrieve the same ones
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        {
+          AssertThrow(m(i, j) == i * j * .5 + .5, ExcInternalError());
+          AssertThrow(m.el(i, j) == i * j * .5 + .5, ExcInternalError());
+        }
+      else
+        {
+          AssertThrow(m.el(i, j) == 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> m(5U, 5U, 3U);
+        test(m);
+      }
+    }
+  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/01.output b/tests/trilinos_tpetra/01.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos_tpetra/02.cc b/tests/trilinos_tpetra/02.cc
new file mode 100644 (file)
index 0000000..e620df8
--- /dev/null
@@ -0,0 +1,106 @@
+// ---------------------------------------------------------------------
+//
+// 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 setting elements in a petsc matrix using
+// LinearAlgebra::TpetraWrappers::SparseMatrix<double>::add()
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &m)
+{
+  // first set a few entries
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        m.set(i, j, 0.);
+
+  // then add values to these entries
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        m.add(i, j, i * j * .5 + .5);
+
+  m.compress(VectorOperation::add);
+
+  // then make sure we retrieve the same ones
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        {
+          AssertThrow(m(i, j) == i * j * .5 + .5, ExcInternalError());
+          AssertThrow(m.el(i, j) == i * j * .5 + .5, ExcInternalError());
+        }
+      else
+        {
+          AssertThrow(m.el(i, j) == 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> m(5U, 5U, 3U);
+        test(m);
+      }
+    }
+  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/02.output b/tests/trilinos_tpetra/02.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos_tpetra/03a.cc b/tests/trilinos_tpetra/03a.cc
new file mode 100644 (file)
index 0000000..3c63425
--- /dev/null
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+//
+// 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 setting elements in a petsc matrix using set() and add()
+// intermixed. this poses PETSc some problems, since one has to flush some
+// buffer in between these two types of operations
+//
+// in contrast to trilinos_03, we set and add the same elements here twice, to
+// get double the original value
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &m)
+{
+  // first set a few entries
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        m.set(i, j, i * j * .5 + .5);
+
+  // then add the same elements again
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        m.add(i, j, i * j * .5 + .5);
+
+  m.compress(VectorOperation::add);
+
+  // then make sure we retrieve the same ones
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        {
+          AssertThrow(m(i, j) == 2 * (i * j * .5 + .5), ExcInternalError());
+          AssertThrow(m.el(i, j) == 2 * (i * j * .5 + .5), ExcInternalError());
+        }
+      else
+        {
+          AssertThrow(m.el(i, j) == 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> m(5U, 5U, 3U);
+        test(m);
+      }
+    }
+  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/03a.output b/tests/trilinos_tpetra/03a.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos_tpetra/04.cc b/tests/trilinos_tpetra/04.cc
new file mode 100644 (file)
index 0000000..86a1942
--- /dev/null
@@ -0,0 +1,82 @@
+// ---------------------------------------------------------------------
+//
+// 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 querying matrix sizes
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &m)
+{
+  AssertThrow(m.m() == 5, ExcInternalError());
+  AssertThrow(m.n() == 5, 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> m(5U, 5U, 3U);
+        test(m);
+      }
+    }
+  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/04.output b/tests/trilinos_tpetra/04.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos_tpetra/05.cc b/tests/trilinos_tpetra/05.cc
new file mode 100644 (file)
index 0000000..7561408
--- /dev/null
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// 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 querying the number of nonzero elements in
+// LinearAlgebra::TpetraWrappers::SparseMatrix<double>
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &m)
+{
+  // first set a few entries. count how many
+  // entries we have
+  unsigned int counter = 0;
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        {
+          m.set(i, j, i * j * .5 + .5);
+          ++counter;
+        }
+
+  m.compress(VectorOperation::insert);
+
+  deallog << m.n_nonzero_elements() << std::endl;
+  Assert(m.n_nonzero_elements() == counter, 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> m(5U, 5U, 3U);
+        test(m);
+      }
+    }
+  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/05.output b/tests/trilinos_tpetra/05.output
new file mode 100644 (file)
index 0000000..169db65
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::8
+DEAL::OK
diff --git a/tests/trilinos_tpetra/09.cc b/tests/trilinos_tpetra/09.cc
new file mode 100644 (file)
index 0000000..4016f89
--- /dev/null
@@ -0,0 +1,105 @@
+// ---------------------------------------------------------------------
+//
+// 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>::operator *=
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &m)
+{
+  // first set a few entries
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        m.set(i, j, i * j * .5 + .5);
+
+  m.compress(VectorOperation::insert);
+
+  // then multiply everything by 1.25 and
+  // make sure we retrieve the values we
+  // expect
+  m *= 1.25;
+
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.m(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0)
+        {
+          AssertThrow(m(i, j) == (i * j * .5 + .5) * 1.25, ExcInternalError());
+          AssertThrow(m.el(i, j) == (i * j * .5 + .5) * 1.25,
+                      ExcInternalError());
+        }
+      else
+        {
+          AssertThrow(m.el(i, j) == 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> m(5U, 5U, 3U);
+        test(m);
+      }
+    }
+  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/09.output b/tests/trilinos_tpetra/09.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos_tpetra/CMakeLists.txt b/tests/trilinos_tpetra/CMakeLists.txt
new file mode 100644 (file)
index 0000000..1a131cb
--- /dev/null
@@ -0,0 +1,6 @@
+cmake_minimum_required(VERSION 3.13.4)
+include(../scripts/setup_testsubproject.cmake)
+project(testsuite CXX)
+if(DEAL_II_WITH_TRILINOS)
+  deal_ii_pickup_tests()
+endif()

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.