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()));
}
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);
-
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
projects / dealii.git / commitdiff