From da05d8ca2dd90a75b7a783b211f9912ce46a624d Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 18 Dec 2014 12:03:18 +0100
Subject: [PATCH] Fix bug in Trilinos SparseMatrix::add/copy_from

We assume that the pointers to a Trilinos sparse matrix remain valid when a matrix
is copied from another one or added and the two matrices share the same sparsity pattern.
This is used in step-32. However, the fix introduced in 8d64256 could not ensure this
(while fixing a few other cases). The check if the two sparsity patterns point to the same
memory is not good (we can have the same sparsity pattern memory but different addresses).

Therefore, the add and copy_from methods need to be more careful to not set up a new
matrix structure when this is not necessary. In particular, the two methods now check
whether the column indices are the same in any case but without requiring calling
methods that change the memory layout.
---
 source/lac/trilinos_sparse_matrix.cc  | 223 ++++++++++++++------------
 tests/trilinos/add_matrices_04.cc     | 105 ++++++++++++
 tests/trilinos/add_matrices_04.output |  36 +++++
 tests/trilinos/add_matrices_05.cc     | 110 +++++++++++++
 tests/trilinos/add_matrices_05.output |  56 +++++++
 tests/trilinos/add_matrices_06.cc     | 123 ++++++++++++++
 tests/trilinos/add_matrices_06.output |   5 +
 7 files changed, 551 insertions(+), 107 deletions(-)
 create mode 100644 tests/trilinos/add_matrices_04.cc
 create mode 100644 tests/trilinos/add_matrices_04.output
 create mode 100644 tests/trilinos/add_matrices_05.cc
 create mode 100644 tests/trilinos/add_matrices_05.output
 create mode 100644 tests/trilinos/add_matrices_06.cc
 create mode 100644 tests/trilinos/add_matrices_06.output

diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc
index aaebb8ea13..ab4b29bc92 100644
--- a/source/lac/trilinos_sparse_matrix.cc
+++ b/source/lac/trilinos_sparse_matrix.cc
@@ -387,31 +387,71 @@ namespace TrilinosWrappers
 
 
   void
-  SparseMatrix::copy_from (const SparseMatrix &m)
+  SparseMatrix::copy_from (const SparseMatrix &rhs)
   {
     nonlocal_matrix.reset();
     nonlocal_matrix_exporter.reset();
 
-    // check whether we need to update the partitioner or can just copy the
-    // data: in case we have the same distribution, we can just copy the data.
-    if (&matrix->Graph() == &m.matrix->Graph() && m.matrix->Filled())
+    // check whether we need to update the whole matrix layout (we have
+    // different maps or if we detect a row where the columns of the two
+    // matrices do not match)
+    bool needs_deep_copy =
+      !matrix->RowMap().SameAs(rhs.matrix->RowMap()) ||
+      !matrix->ColMap().SameAs(rhs.matrix->ColMap()) ||
+      !matrix->DomainMap().SameAs(rhs.matrix->DomainMap()) ||
+      n_nonzero_elements() != rhs.n_nonzero_elements();
+    if (!needs_deep_copy)
       {
-        std::memcpy(matrix->ExpertExtractValues(), m.matrix->ExpertExtractValues(),
-                    matrix->NumMyNonzeros()*sizeof(*matrix->ExpertExtractValues()));
+        const std::pair<size_type, size_type>
+        local_range = rhs.local_range();
+
+        int ierr;
+        // Try to copy all the rows of the matrix one by one. In case of error
+        // (i.e., the column indices are different), we need to abort and blow
+        // away the matrix.
+        for (size_type row=local_range.first; row < local_range.second; ++row)
+          {
+            const int row_local =
+              matrix->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row));
+
+            int n_entries, rhs_n_entries;
+            TrilinosScalar *value_ptr, *rhs_value_ptr;
+            int *index_ptr, *rhs_index_ptr;
+            ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,
+                                                 rhs_value_ptr, rhs_index_ptr);
+            Assert (ierr == 0, ExcTrilinosError(ierr));
+
+            ierr = matrix->ExtractMyRowView (row_local, n_entries, value_ptr,
+                                             index_ptr);
+            Assert (ierr == 0, ExcTrilinosError(ierr));
+
+            if (n_entries != rhs_n_entries ||
+                std::memcmp(static_cast<void *>(index_ptr),
+                            static_cast<void *>(rhs_index_ptr),
+                            sizeof(int)*n_entries) != 0)
+              {
+                needs_deep_copy = true;
+                break;
+              }
+
+            for (int i=0; i<n_entries; ++i)
+              value_ptr[i] = rhs_value_ptr[i];
+          }
       }
-    else
+
+    if (needs_deep_copy)
       {
-        column_space_map.reset (new Epetra_Map (m.domain_partitioner()));
+        column_space_map.reset (new Epetra_Map (rhs.domain_partitioner()));
 
         // release memory before reallocation
         matrix.reset ();
-        matrix.reset (new Epetra_FECrsMatrix(*m.matrix));
-      }
+        matrix.reset (new Epetra_FECrsMatrix(*rhs.matrix));
 
-    if (m.nonlocal_matrix.get() != 0)
-      nonlocal_matrix.reset(new Epetra_CrsMatrix(Copy, m.nonlocal_matrix->Graph()));
+        matrix->FillComplete(*column_space_map, matrix->RowMap());
+      }
 
-    matrix->FillComplete(*column_space_map, matrix->RowMap());
+    if (rhs.nonlocal_matrix.get() != 0)
+      nonlocal_matrix.reset(new Epetra_CrsMatrix(Copy, rhs.nonlocal_matrix->Graph()));
   }
 
 
@@ -1522,108 +1562,77 @@ namespace TrilinosWrappers
   SparseMatrix::add (const TrilinosScalar  factor,
                      const SparseMatrix   &rhs)
   {
-    Assert (rhs.m() == m(), ExcDimensionMismatch (rhs.m(), m()));
-    Assert (rhs.n() == n(), ExcDimensionMismatch (rhs.n(), n()));
+    AssertDimension (rhs.m(), m());
+    AssertDimension (rhs.n(), n());
+    AssertDimension (rhs.local_range().first, local_range().first);
+    AssertDimension (rhs.local_range().second, local_range().second);
+    Assert(matrix->RowMap().SameAs(rhs.matrix->RowMap()),
+           ExcMessage("Can only add matrices with same distribution of rows"));
+    Assert(matrix->Filled() && rhs.matrix->Filled(),
+           ExcMessage("Addition of matrices only allowed if matrices are "
+                      "filled, i.e., compress() has been called"));
 
     const std::pair<size_type, size_type>
     local_range = rhs.local_range();
+    const bool same_col_map = matrix->ColMap().SameAs(rhs.matrix->ColMap());
 
     int ierr;
-
-    // If both matrices have been transformed to local index space (in
-    // Trilinos speak: they are filled) and the matrices are based on the same
-    // sparsity pattern, we can extract views of the column data on both
-    // matrices and simply manipulate the values that are addressed by the
-    // pointers.
-    if (matrix->Filled() == true &&
-        rhs.matrix->Filled() == true &&
-        &matrix->Graph() == &rhs.matrix->Graph())
-      for (size_type row=local_range.first;
-           row < local_range.second; ++row)
-        {
-          Assert (matrix->NumGlobalEntries(row) ==
-                  rhs.matrix->NumGlobalEntries(row),
-                  ExcDimensionMismatch(matrix->NumGlobalEntries(row),
-                                       rhs.matrix->NumGlobalEntries(row)));
-
-          const TrilinosWrappers::types::int_type row_local =
-            matrix->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row));
-          int n_entries, rhs_n_entries;
-          TrilinosScalar *value_ptr, *rhs_value_ptr;
-
-          // In debug mode, we want to check whether the indices really are
-          // the same in the calling matrix and the input matrix. The reason
-          // for doing this only in debug mode is that both extracting indices
-          // and comparing indices is relatively slow compared to just working
-          // with the values.
-#ifdef DEBUG
-          int *index_ptr, *rhs_index_ptr;
-          ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,
-                                               rhs_value_ptr, rhs_index_ptr);
-          Assert (ierr == 0, ExcTrilinosError(ierr));
-
-          ierr = matrix->ExtractMyRowView (row_local, n_entries, value_ptr,
-                                           index_ptr);
-          Assert (ierr == 0, ExcTrilinosError(ierr));
-#else
-          rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,rhs_value_ptr);
-          matrix->ExtractMyRowView (row_local, n_entries, value_ptr);
-#endif
-
-          AssertDimension (n_entries, rhs_n_entries);
-
-          for (TrilinosWrappers::types::int_type i=0; i<n_entries; ++i)
-            {
-              *value_ptr++ += *rhs_value_ptr++ * factor;
-#ifdef DEBUG
-              Assert (*index_ptr++ == *rhs_index_ptr++,
-                      ExcInternalError());
-#endif
-            }
-        }
-    // If we have different sparsity patterns, we have to be more careful (in
-    // particular when we use multiple processors) and extract a copy of the
-    // row data, multiply it by the factor and then add it to the matrix using
-    // the respective add() function.
-    else
+    for (size_type row=local_range.first; row < local_range.second; ++row)
       {
-        int max_row_length = 0;
-        for (size_type row=local_range.first;
-             row < local_range.second; ++row)
-          max_row_length
-            = std::max (max_row_length,rhs.matrix->NumGlobalEntries(row));
-
-        std::vector<TrilinosScalar> values (max_row_length);
-        std::vector<TrilinosWrappers::types::int_type> column_indices (max_row_length);
-        for (size_type row=local_range.first;
-             row < local_range.second; ++row)
+        const int row_local =
+          matrix->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row));
+
+        // First get a view to the matrix columns of both matrices. Note that
+        // the data is in local index spaces so we need to be careful not only
+        // to compare column indices in case they are derived from the same
+        // map.
+        int n_entries, rhs_n_entries;
+        TrilinosScalar *value_ptr, *rhs_value_ptr;
+        int *index_ptr, *rhs_index_ptr;
+        ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,
+                                             rhs_value_ptr, rhs_index_ptr);
+        Assert (ierr == 0, ExcTrilinosError(ierr));
+
+        ierr = matrix->ExtractMyRowView (row_local, n_entries, value_ptr,
+                                         index_ptr);
+        Assert (ierr == 0, ExcTrilinosError(ierr));
+        bool expensive_checks = (n_entries != rhs_n_entries || !same_col_map);
+        if (!expensive_checks)
           {
-            int n_entries;
-            ierr = rhs.matrix->Epetra_CrsMatrix::ExtractGlobalRowCopy
-                   (static_cast<TrilinosWrappers::types::int_type>(row),
-                    max_row_length, n_entries, &values[0], &column_indices[0]);
-            Assert (ierr == 0, ExcTrilinosError(ierr));
-
-            // Filter away zero elements
-            unsigned int n_actual_entries = 0.;
-            for (int i=0; i<n_entries; ++i)
-              if (std::abs(values[i]) != 0.)
-                {
-                  column_indices[n_actual_entries] = column_indices[i];
-                  values[n_actual_entries++] = values[i] * factor;
-                }
-
-            ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues
-                   (static_cast<TrilinosWrappers::types::int_type>(row),
-                    n_actual_entries, &values[0], &column_indices[0]);
-            Assert (ierr == 0,
-                    ExcMessage("Adding the entries from the other matrix "
-                               "failed, possibly because the sparsity pattern "
-                               "of that matrix includes more elements than the "
-                               "calling matrix, which is not allowed."));
+            // check if the column indices are the same. If yes, can simply
+            // copy over the data.
+            expensive_checks = std::memcmp(static_cast<void *>(index_ptr),
+                                           static_cast<void *>(rhs_index_ptr),
+                                           sizeof(int)*n_entries) != 0;
+            if (!expensive_checks)
+              for (int i=0; i<n_entries; ++i)
+                value_ptr[i] += rhs_value_ptr[i] * factor;
+          }
+        // Now to the expensive case where we need to check all column indices
+        // against each other (transformed into global index space) and where
+        // we need to make sure that all entries we are about to add into the
+        // lhs matrix actually exist
+        if (expensive_checks)
+          {
+            for (int i=0; i<rhs_n_entries; ++i)
+              {
+                if (rhs_value_ptr[i] == 0.)
+                  continue;
+                const TrilinosWrappers::types::int_type rhs_global_col =
+                  global_column_index(*rhs.matrix, rhs_index_ptr[i]);
+                int local_col = matrix->ColMap().LID(rhs_global_col);
+                int *local_index = Utilities::lower_bound(index_ptr,
+                                                          index_ptr+n_entries,
+                                                          local_col);
+                Assert(local_index != index_ptr + n_entries &&
+                       *local_index == local_col,
+                       ExcMessage("Adding the entries from the other matrix "
+                                  "failed, because the sparsity pattern "
+                                  "of that matrix includes more elements than the "
+                                  "calling matrix, which is not allowed."));
+                value_ptr[local_index-index_ptr] += factor * rhs_value_ptr[i];
+              }
           }
-        compress (VectorOperation::add);
-
       }
   }
 
diff --git a/tests/trilinos/add_matrices_04.cc b/tests/trilinos/add_matrices_04.cc
new file mode 100644
index 0000000000..0cc21083c6
--- /dev/null
+++ b/tests/trilinos/add_matrices_04.cc
@@ -0,0 +1,105 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check SparseMatrix::add(SparseMatrix) in a few variants
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+
+
+void test (TrilinosWrappers::SparseMatrix &m)
+{
+  TrilinosWrappers::SparseMatrix m2(m.m(), m.n(), 0);
+
+  // first set a few entries one-by-one
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+      if ((i+2*j+1) % 3 == 0)
+        {
+          m.set (i,j, i*j*.5+.5);
+          m2.set (i,j, 1.);
+        }
+
+  m.compress (VectorOperation::insert);
+  m2.compress(VectorOperation::insert);
+
+  m.print(deallog.get_file_stream());
+  deallog << std::endl;
+
+  m.add(2, m2);
+  m.print(deallog.get_file_stream());
+
+  deallog << std::endl;
+
+  m.add(-2, m2);
+  m.print(deallog.get_file_stream());
+
+  m.copy_from(m2);
+  m.add(-1., m2);
+  
+  deallog << std::endl << "Frobenius norm: " << m.frobenius_norm() << std::endl;
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+  try
+    {
+      {
+        TrilinosWrappers::SparseMatrix m (5U,6U,3U);
+
+        test (m);
+      }
+    }
+  catch (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/add_matrices_04.output b/tests/trilinos/add_matrices_04.output
new file mode 100644
index 0000000000..20a27af6d3
--- /dev/null
+++ b/tests/trilinos/add_matrices_04.output
@@ -0,0 +1,36 @@
+
+(0,1) 0.500000
+(0,4) 0.500000
+(1,2) 1.50000
+(1,5) 3.00000
+(2,0) 0.500000
+(2,3) 3.50000
+(3,1) 2.00000
+(3,4) 6.50000
+(4,2) 4.50000
+(4,5) 10.5000
+DEAL::
+(0,1) 2.50000
+(0,4) 2.50000
+(1,2) 3.50000
+(1,5) 5.00000
+(2,0) 2.50000
+(2,3) 5.50000
+(3,1) 4.00000
+(3,4) 8.50000
+(4,2) 6.50000
+(4,5) 12.5000
+DEAL::
+(0,1) 0.500000
+(0,4) 0.500000
+(1,2) 1.50000
+(1,5) 3.00000
+(2,0) 0.500000
+(2,3) 3.50000
+(3,1) 2.00000
+(3,4) 6.50000
+(4,2) 4.50000
+(4,5) 10.5000
+DEAL::
+DEAL::Frobenius norm: 0
+DEAL::OK
diff --git a/tests/trilinos/add_matrices_05.cc b/tests/trilinos/add_matrices_05.cc
new file mode 100644
index 0000000000..0f1fc5bf57
--- /dev/null
+++ b/tests/trilinos/add_matrices_05.cc
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check SparseMatrix::add(SparseMatrix) where the other matrix has either
+// more or less entries (but with reasonable entries in any case)
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+
+
+void test (TrilinosWrappers::SparseMatrix &m)
+{
+  TrilinosWrappers::SparseMatrix m2(m.m(), m.n(), 0);
+
+  // first set a few entries one-by-one
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+      if ((i+2*j+1) % 3 == 0)
+        {
+          m.set (i,j, i*j*.5+.5);
+          m2.set (i,j, 1.);
+        }
+      else if (j % 2 == 0)
+        m2.set(i,j, 0);
+
+  m.compress (VectorOperation::insert);
+  m2.compress(VectorOperation::insert);
+
+  deallog << "Matrix nonzeros: " << m.n_nonzero_elements() <<  " "
+          << m2.n_nonzero_elements() << std::endl;
+
+  m.print(deallog.get_file_stream());
+  deallog << std::endl;
+
+  m.add(1, m2);
+  m.print(deallog.get_file_stream());
+  deallog << std::endl;
+
+  m.add(-1, m2);
+  m2.add(-1, m);
+  m2.print(deallog.get_file_stream());
+  deallog << std::endl;
+
+  m.add(-1, m2);
+  m.print(deallog.get_file_stream());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+  try
+    {
+      {
+        TrilinosWrappers::SparseMatrix m (5U,6U,3U);
+
+        test (m);
+      }
+    }
+  catch (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/add_matrices_05.output b/tests/trilinos/add_matrices_05.output
new file mode 100644
index 0000000000..9ef10410d8
--- /dev/null
+++ b/tests/trilinos/add_matrices_05.output
@@ -0,0 +1,56 @@
+
+DEAL::Matrix nonzeros: 10 20
+(0,1) 0.500000
+(0,4) 0.500000
+(1,2) 1.50000
+(1,5) 3.00000
+(2,0) 0.500000
+(2,3) 3.50000
+(3,1) 2.00000
+(3,4) 6.50000
+(4,2) 4.50000
+(4,5) 10.5000
+DEAL::
+(0,1) 1.50000
+(0,4) 1.50000
+(1,2) 2.50000
+(1,5) 4.00000
+(2,0) 1.50000
+(2,3) 4.50000
+(3,1) 3.00000
+(3,4) 7.50000
+(4,2) 5.50000
+(4,5) 11.5000
+DEAL::
+(0,0) 0.00000
+(0,1) 0.500000
+(0,2) 0.00000
+(0,4) 0.500000
+(1,0) 0.00000
+(1,2) -0.500000
+(1,4) 0.00000
+(1,5) -2.00000
+(2,0) 0.500000
+(2,2) 0.00000
+(2,3) -2.50000
+(2,4) 0.00000
+(3,0) 0.00000
+(3,1) -1.00000
+(3,2) 0.00000
+(3,4) -5.50000
+(4,0) 0.00000
+(4,2) -3.50000
+(4,4) 0.00000
+(4,5) -9.50000
+DEAL::
+(0,1) 0.00000
+(0,4) 0.00000
+(1,2) 2.00000
+(1,5) 5.00000
+(2,0) 0.00000
+(2,3) 6.00000
+(3,1) 3.00000
+(3,4) 12.0000
+(4,2) 8.00000
+(4,5) 20.0000
+DEAL::OK
diff --git a/tests/trilinos/add_matrices_06.cc b/tests/trilinos/add_matrices_06.cc
new file mode 100644
index 0000000000..117a1db946
--- /dev/null
+++ b/tests/trilinos/add_matrices_06.cc
@@ -0,0 +1,123 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that SparseMatrix::add(SparseMatrix) and
+// SparseMatrix::copy_from(SparseMatrix) do not destroy memory when called
+// from matrices with the same sparsity pattern and thus, structures depending
+// on these matrices such as PreconditionJacobi do not get destroyed
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <fstream>
+#include <iostream>
+
+
+void test (TrilinosWrappers::SparseMatrix &m)
+{
+  TrilinosWrappers::SparseMatrix m2(m.m(), m.n(), 0), m3(m.m(), m.n(), 0);
+
+  // first set a few entries one-by-one
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+      if ((i+2*j+1) % 3 == 0 || i == j)
+        {
+          m.set (i,j, i*j*.5+.5);
+          m2.set (i,j, 1.);
+          m3.set (i,j, -0.1);
+        }
+
+  m.compress ();
+  m2.compress();
+  m3.compress();
+
+  deallog << "Matrix nonzeros: " << m.n_nonzero_elements() <<  " "
+          << m2.n_nonzero_elements() <<  " "
+          << m3.n_nonzero_elements() << std::endl;
+
+  m.copy_from(m2);
+  m.add(0.12, m3);
+
+  Vector<double> src(m.m()), dst(m.m());
+  src = 1.;
+
+  TrilinosWrappers::PreconditionJacobi prec;
+  prec.initialize(m);
+
+  m.vmult(dst, src);
+  deallog << "Vector norm: " << dst.l2_norm() << " ";
+
+  prec.vmult(dst, src);
+  deallog << dst.l2_norm() << std::endl;
+
+  m.copy_from(m2);
+  m.add(0.06, m3);
+
+  m.vmult(dst, src);
+  deallog << "Vector norm: " << dst.l2_norm() << " ";
+
+  prec.vmult(dst, src);
+  deallog << dst.l2_norm() << std::endl;
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+  try
+    {
+      {
+        TrilinosWrappers::SparseMatrix m (5U,5U,3U);
+
+        test (m);
+      }
+    }
+  catch (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/add_matrices_06.output b/tests/trilinos/add_matrices_06.output
new file mode 100644
index 0000000000..f1e4081029
--- /dev/null
+++ b/tests/trilinos/add_matrices_06.output
@@ -0,0 +1,5 @@
+
+DEAL::Matrix nonzeros: 13 13 13
+DEAL::Vector norm: 5.84509 2.26323
+DEAL::Vector norm: 5.88058 2.26323
+DEAL::OK
-- 
2.39.5