From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Mon, 22 Sep 2014 11:56:48 +0000 (+0200)
Subject: Fix TrilinosWrappers::SparseMatrix::add(factor, other_matrix).
X-Git-Tag: v8.2.0-rc1~134^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F168%2Fhead;p=dealii.git

Fix TrilinosWrappers::SparseMatrix::add(factor, other_matrix).
---

diff --git a/doc/news/changes.h b/doc/news/changes.h
index bb77d0f02a..85e7a1c6f2 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -305,6 +305,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: The function TrilinosWrappers::SparseMatrix::add(double factor,
+  SparseMatrix &rhs) produced wrong results and ran into an exception if the
+  rhs matrix included off-processor column entries. This is now fixed.
+  <br>
+  (Martin Kronbichler, 2014/09/21)
+  </li>
+
   <li> New: The function Threads::Task::joinable() can be used to verify whether
   a task object can be joined or not.
   <br>
diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc
index 587483a02a..018cfbca98 100644
--- a/source/lac/trilinos_sparse_matrix.cc
+++ b/source/lac/trilinos_sparse_matrix.cc
@@ -1530,21 +1530,14 @@ namespace TrilinosWrappers
 
     int ierr;
 
-    // If both matrices have been transformed
-    // to local index space (in Trilinos
-    // speak: they are filled), we're having
-    // matrices based on the same indices
-    // with the same number of nonzeros
-    // (actually, we'd need sparsity pattern,
-    // but that is too expensive to check),
-    // we can extract views of the column
-    // data on both matrices and simply
-    // manipulate the values that are
-    // addressed by the pointers.
+    // 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 &&
-        this->local_range() == local_range &&
-        matrix->NumMyNonzeros() == rhs.matrix->NumMyNonzeros())
+        &matrix->Graph() == &rhs.matrix->Graph())
       for (size_type row=local_range.first;
            row < local_range.second; ++row)
         {
@@ -1558,14 +1551,11 @@ namespace TrilinosWrappers
           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.
+          // 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,
@@ -1580,8 +1570,7 @@ namespace TrilinosWrappers
           matrix->ExtractMyRowView (row_local, n_entries, value_ptr);
 #endif
 
-          AssertThrow (n_entries == rhs_n_entries,
-                       ExcDimensionMismatch (n_entries, rhs_n_entries));
+          AssertDimension (n_entries, rhs_n_entries);
 
           for (TrilinosWrappers::types::int_type i=0; i<n_entries; ++i)
             {
@@ -1592,13 +1581,10 @@ namespace TrilinosWrappers
 #endif
             }
         }
-    // If we have different sparsity patterns
-    // (expressed by a different number of
-    // nonzero elements), we have to be more
-    // careful 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.
+    // 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
       {
         int max_row_length = 0;
@@ -1608,58 +1594,36 @@ namespace TrilinosWrappers
             = std::max (max_row_length,rhs.matrix->NumGlobalEntries(row));
 
         std::vector<TrilinosScalar> values (max_row_length);
-
-        if (matrix->Filled() == true && rhs.matrix->Filled() == true &&
-            this->local_range() == local_range)
-          for (size_type row=local_range.first;
-               row < local_range.second; ++row)
-            {
-              std::vector<int> column_indices (max_row_length);
-              const int row_local =
-                matrix->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row));
-              int n_entries;
-
-              ierr = rhs.matrix->ExtractMyRowCopy (row_local, max_row_length,
-                                                   n_entries,
-                                                   &values[0],
-                                                   &column_indices[0]);
-              Assert (ierr == 0, ExcTrilinosError(ierr));
-
-              for (TrilinosWrappers::types::int_type i=0; i<n_entries; ++i)
-                values[i] *= factor;
-
-              TrilinosScalar *value_ptr = &values[0];
-
-              ierr = matrix->SumIntoMyValues (row_local, n_entries, value_ptr,
-                                              &column_indices[0]);
-              Assert (ierr == 0, ExcTrilinosError(ierr));
-            }
-        else
+        std::vector<TrilinosWrappers::types::int_type> column_indices (max_row_length);
+        for (size_type row=local_range.first;
+             row < local_range.second; ++row)
           {
-            //TODO check that is normal that column_indices in the if is an
-            //int while the column_indices in the else is a
-            //TrilinosWrappers::types::int_type
-            std::vector<TrilinosWrappers::types::int_type> column_indices (max_row_length);
-            for (size_type row=local_range.first;
-                 row < local_range.second; ++row)
-              {
-                int n_entries;
-                ierr = rhs.matrix->Epetra_CrsMatrix::ExtractGlobalRowCopy
-                       ((TrilinosWrappers::types::int_type)row, max_row_length,
-                        n_entries, &values[0], &column_indices[0]);
-                Assert (ierr == 0, ExcTrilinosError(ierr));
-
-                for (TrilinosWrappers::types::int_type i=0; i<n_entries; ++i)
-                  values[i] *= factor;
-
-                ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues
-                       ((TrilinosWrappers::types::int_type)row, n_entries,
-                        &values[0], &column_indices[0]);
-                Assert (ierr == 0, ExcTrilinosError(ierr));
-              }
-            compress ();
+            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."));
           }
+        compress (VectorOperation::add);
+
       }
   }
 
diff --git a/tests/trilinos/add_matrices_01.cc b/tests/trilinos/add_matrices_01.cc
new file mode 100644
index 0000000000..9b29911c1d
--- /dev/null
+++ b/tests/trilinos/add_matrices_01.cc
@@ -0,0 +1,133 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test SparseMatrix::add(factor, SparseMatrix) based on different sparsity
+// pattern (and where the sparsity pattern of the calling matrix is more
+// inclusive which is necessary for the operation to succeed)
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <fstream>
+#include <iostream>
+
+
+void test ()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+  if (myid==0) deallog << "numproc=" << numproc << std::endl;
+
+
+  // each processor owns 3 indices
+  IndexSet local_owned(numproc*3);
+  local_owned.add_range(myid*3,myid*3+3);
+
+  // Create sparsity patterns
+  TrilinosWrappers::SparsityPattern sp1(local_owned, MPI_COMM_WORLD),
+                   sp2(local_owned, MPI_COMM_WORLD);
+
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    for (unsigned int j=0; j<local_owned.size(); ++j)
+      if ((i+j) % 2 == 1)
+        {
+          sp1.add (i,j);
+          if (j%2 == 0)
+            sp2.add(i,j);
+        }
+
+  sp1.compress ();
+  sp2.compress();
+
+  // create matrices by adding some elements into the respective positions
+  TrilinosWrappers::SparseMatrix m1(sp1), m2(sp2);
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    for (unsigned int j=0; j<local_owned.size(); ++j)
+      if ((i+j) % 2 == 1)
+        {
+          m1.add (i,j, i+j);
+          if (j%2 == 0)
+            m2.add(i,j, i+2*j+1);
+        }
+  m1.compress(VectorOperation::add);
+  m2.compress(VectorOperation::add);
+
+  m1.add(2, m2);
+
+  // Check for correctness of entries (all floating point comparisons should
+  // be exact)
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    for (unsigned int j=0; j<local_owned.size(); ++j)
+      if ((i+j) % 2 == 1 && j%2 == 0)
+        {
+          Assert(m1.el(i,j) == (double)i+j+2*i+4*j+2, ExcInternalError());
+        }
+      else if ((i+j) % 2 == 1)
+        {
+          Assert(m1.el(i,j) == (double)i+j, ExcInternalError());
+        }
+      else
+        {
+          Assert(m1.el(i,j) == 0., ExcInternalError());
+        }
+
+  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
+    {
+      test ();
+    }
+  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_01.with_mpi=true.mpirun=3.output b/tests/trilinos/add_matrices_01.with_mpi=true.mpirun=3.output
new file mode 100644
index 0000000000..763b6fa8be
--- /dev/null
+++ b/tests/trilinos/add_matrices_01.with_mpi=true.mpirun=3.output
@@ -0,0 +1,4 @@
+
+DEAL::OK
+proc=3
+DEAL::OK
diff --git a/tests/trilinos/add_matrices_02.cc b/tests/trilinos/add_matrices_02.cc
new file mode 100644
index 0000000000..7fc4101f6f
--- /dev/null
+++ b/tests/trilinos/add_matrices_02.cc
@@ -0,0 +1,125 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test SparseMatrix::add(factor, SparseMatrix) based on different sparsity
+// patterns (neither sparsity pattern is more inclusive than the other, but
+// the entry in the rhs matrix in the 'offending' position is zero).
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <fstream>
+#include <iostream>
+
+
+void test ()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+  if (myid==0) deallog << "numproc=" << numproc << std::endl;
+
+
+  // each processor owns 3 indices
+  IndexSet local_owned(numproc*3);
+  local_owned.add_range(myid*3,myid*3+3);
+
+  // Create sparsity patterns
+  TrilinosWrappers::SparsityPattern sp1(local_owned, MPI_COMM_WORLD),
+                   sp2(local_owned, MPI_COMM_WORLD);
+
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    {
+      sp1.add(i,i);
+      sp2.add(i,i);
+    }
+  if (myid == 0)
+    {
+      sp1.add(0,1);
+      sp2.add(1,0);
+    }
+
+  sp1.compress ();
+  sp2.compress();
+
+  // create matrices by adding some elements into the respective positions
+  TrilinosWrappers::SparseMatrix m1(sp1), m2(sp2);
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    {
+      m1.add (i,i, i+2);
+      m2.add (i,i, 4);
+    }
+  if (myid == 0)
+    m1.add(0,1,3);
+
+  m1.compress(VectorOperation::add);
+  m2.compress(VectorOperation::add);
+
+  m1.add(2, m2);
+
+  // Check for correctness of entries (all floating point comparisons should
+  // be exact)
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    Assert(m1.el(i,i) == (double)i+2+8, ExcInternalError());
+  if (myid == 0)
+    Assert(m1.el(0,1) == 3., ExcInternalError());
+
+  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
+    {
+      test ();
+    }
+  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_02.with_mpi=true.mpirun=3.output b/tests/trilinos/add_matrices_02.with_mpi=true.mpirun=3.output
new file mode 100644
index 0000000000..763b6fa8be
--- /dev/null
+++ b/tests/trilinos/add_matrices_02.with_mpi=true.mpirun=3.output
@@ -0,0 +1,4 @@
+
+DEAL::OK
+proc=3
+DEAL::OK
diff --git a/tests/trilinos/add_matrices_03.cc b/tests/trilinos/add_matrices_03.cc
new file mode 100644
index 0000000000..7cd4e15581
--- /dev/null
+++ b/tests/trilinos/add_matrices_03.cc
@@ -0,0 +1,128 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test SparseMatrix::add(factor, SparseMatrix) based on matrices of the same
+// sparsity pattern
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <fstream>
+#include <iostream>
+
+
+void test ()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+  if (myid==0) deallog << "numproc=" << numproc << std::endl;
+
+
+  // each processor owns 3 indices
+  IndexSet local_owned(numproc*3);
+  local_owned.add_range(myid*3,myid*3+3);
+
+  // Create sparsity patterns
+  TrilinosWrappers::SparsityPattern sp(local_owned, MPI_COMM_WORLD);
+
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    for (unsigned int j=0; j<local_owned.size(); ++j)
+      if ((i+j) % 2 == 1)
+        {
+          sp.add (i,j);
+        }
+
+  sp.compress ();
+
+  // create matrices by adding some elements into the respective positions
+  TrilinosWrappers::SparseMatrix m1(sp), m2(sp);
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    for (unsigned int j=0; j<local_owned.size(); ++j)
+      if ((i+j) % 2 == 1)
+        {
+          m1.add (i,j, i+j);
+          if (j%2 == 0)
+            m2.add(i,j, i+2*j+1);
+        }
+  m1.compress(VectorOperation::add);
+  m2.compress(VectorOperation::add);
+
+  m1.add(2, m2);
+
+  // Check for correctness of entries (all floating point comparisons should
+  // be exact)
+  for (unsigned int i=myid*3; i<myid*3+3; ++i)
+    for (unsigned int j=0; j<local_owned.size(); ++j)
+      if ((i+j) % 2 == 1 && j%2 == 0)
+        {
+          Assert(m1.el(i,j) == (double)i+j+2*i+4*j+2, ExcInternalError());
+        }
+      else if ((i+j) % 2 == 1)
+        {
+          Assert(m1.el(i,j) == (double)i+j, ExcInternalError());
+        }
+      else
+        {
+          Assert(m1.el(i,j) == 0., ExcInternalError());
+        }
+
+  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
+    {
+      test ();
+    }
+  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_03.with_mpi=true.mpirun=3.output b/tests/trilinos/add_matrices_03.with_mpi=true.mpirun=3.output
new file mode 100644
index 0000000000..763b6fa8be
--- /dev/null
+++ b/tests/trilinos/add_matrices_03.with_mpi=true.mpirun=3.output
@@ -0,0 +1,4 @@
+
+DEAL::OK
+proc=3
+DEAL::OK