From 8d64256d8f4dccaea27b94bc2eb8858ce23ef103 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 22 Sep 2014 13:56:48 +0200 Subject: [PATCH] Fix TrilinosWrappers::SparseMatrix::add(factor, other_matrix). --- doc/news/changes.h | 7 + source/lac/trilinos_sparse_matrix.cc | 122 ++++++---------- tests/trilinos/add_matrices_01.cc | 133 ++++++++++++++++++ ..._matrices_01.with_mpi=true.mpirun=3.output | 4 + tests/trilinos/add_matrices_02.cc | 125 ++++++++++++++++ ..._matrices_02.with_mpi=true.mpirun=3.output | 4 + tests/trilinos/add_matrices_03.cc | 128 +++++++++++++++++ ..._matrices_03.with_mpi=true.mpirun=3.output | 4 + 8 files changed, 448 insertions(+), 79 deletions(-) create mode 100644 tests/trilinos/add_matrices_01.cc create mode 100644 tests/trilinos/add_matrices_01.with_mpi=true.mpirun=3.output create mode 100644 tests/trilinos/add_matrices_02.cc create mode 100644 tests/trilinos/add_matrices_02.with_mpi=true.mpirun=3.output create mode 100644 tests/trilinos/add_matrices_03.cc create mode 100644 tests/trilinos/add_matrices_03.with_mpi=true.mpirun=3.output 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.

Specific improvements

    +
  1. 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. +
    + (Martin Kronbichler, 2014/09/21) +
    2. +
  2. New: The function Threads::Task::joinable() can be used to verify whether a task object can be joined or not.
    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; iNumGlobalEntries(row)); std::vector 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 column_indices (max_row_length); - const int row_local = - matrix->RowMap().LID(static_cast(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; iSumIntoMyValues (row_local, n_entries, value_ptr, - &column_indices[0]); - Assert (ierr == 0, ExcTrilinosError(ierr)); - } - else + std::vector 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 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; iEpetra_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(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; iEpetra_CrsMatrix::SumIntoGlobalValues + (static_cast(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 +#include +#include +#include +#include + + +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 +#include +#include +#include +#include + + +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 +#include +#include +#include +#include + + +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