From da05d8ca2dd90a75b7a783b211f9912ce46a624d Mon Sep 17 00:00:00 2001 From: Martin Kronbichler 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 + 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(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(index_ptr), + static_cast(rhs_index_ptr), + sizeof(int)*n_entries) != 0) + { + needs_deep_copy = true; + break; + } + + for (int i=0; iGraph())); + 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 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(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; iNumGlobalEntries(row)); - - std::vector values (max_row_length); - std::vector 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(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(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.")); + // check if the column indices are the same. If yes, can simply + // copy over the data. + expensive_checks = std::memcmp(static_cast(index_ptr), + static_cast(rhs_index_ptr), + sizeof(int)*n_entries) != 0; + if (!expensive_checks) + for (int i=0; iColMap().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 +#include +#include +#include + + +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 +#include +#include +#include + + +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 +#include +#include +#include +#include + + +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 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