<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>
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)
{
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,
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)
{
#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;
= 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);
+
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+proc=3
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+// $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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+proc=3
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+// $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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+proc=3
+DEAL::OK
projects / dealii.git / commitdiff