From: Martin Kronbichler <kronbichler@lnm.mw.tum.de> Date: Wed, 5 Aug 2015 06:52:51 +0000 (+0200) Subject: Improve documentation of TrilinosWrappers::SparseMatrix::set X-Git-Tag: v8.4.0-rc2~674^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a5acacf18a66f1666ff39bccfb7426aab0f2c57c;p=dealii.git Improve documentation of TrilinosWrappers::SparseMatrix::set --- diff --git a/include/deal.II/lac/trilinos_sparse_matrix.h b/include/deal.II/lac/trilinos_sparse_matrix.h index 6abe7d3cfb..bfc51149a1 100644 --- a/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_sparse_matrix.h @@ -1112,15 +1112,21 @@ namespace TrilinosWrappers * * This function is able to insert new elements into the matrix as long as * compress() has not been called, so the sparsity pattern will be - * extended. When compress() is called for the first time, then this is no - * longer possible and an insertion of elements at positions which have - * not been initialized will throw an exception. Note that in case - * elements need to be inserted, it is mandatory that elements are - * inserted only once. Otherwise, the elements will actually be added in - * the end (since it is not possible to efficiently find values to the - * same entry before compress() has been called). In the case that an - * element is set more than once, initialize the matrix with a sparsity - * pattern first. + * extended. When compress() is called for the first time (or in case the + * matrix is initialized from a sparsity pattern), no new elements can be + * added and an insertion of elements at positions which have not been + * initialized will throw an exception. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Epetra_FECrsMatrix data + * structure: If the same matrix entry is inserted more than once, the + * matrix entries will be added upon calling compress() (since Epetra does + * not track values to the same entry before the final compress() is + * called), even if VectorOperation::insert is specified as argument to + * compress(). In the case you cannot make sure that matrix entries are + * only set once, initialize the matrix with a sparsity pattern to fix the + * matrix structure before inserting elements. */ void set (const size_type i, const size_type j, @@ -1137,14 +1143,26 @@ namespace TrilinosWrappers * * This function is able to insert new elements into the matrix as long as * compress() has not been called, so the sparsity pattern will be - * extended. When compress() is called for the first time, then this is no - * longer possible and an insertion of elements at positions which have - * not been initialized will throw an exception. + * extended. After compress() has been called for the first time or the + * matrix has been initialized from a sparsity pattern, extending the + * sparsity pattern is no longer possible and an insertion of elements at + * positions which have not been initialized will throw an exception. * * The optional parameter <tt>elide_zero_values</tt> can be used to * specify whether zero values should be inserted anyway or they should be * filtered away. The default value is <tt>false</tt>, i.e., even zero * values are inserted/replaced. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Epetra_FECrsMatrix data + * structure: If the same matrix entry is inserted more than once, the + * matrix entries will be added upon calling compress() (since Epetra does + * not track values to the same entry before the final compress() is + * called), even if VectorOperation::insert is specified as argument to + * compress(). In the case you cannot make sure that matrix entries are + * only set once, initialize the matrix with a sparsity pattern to fix the + * matrix structure before inserting elements. */ void set (const std::vector<size_type> &indices, const FullMatrix<TrilinosScalar> &full_matrix, @@ -1166,14 +1184,26 @@ namespace TrilinosWrappers * * This function is able to insert new elements into the matrix as long as * compress() has not been called, so the sparsity pattern will be - * extended. When compress() is called for the first time, then this is no - * longer possible and an insertion of elements at positions which have - * not been initialized will throw an exception. + * extended. After compress() has been called for the first time or the + * matrix has been initialized from a sparsity pattern, extending the + * sparsity pattern is no longer possible and an insertion of elements at + * positions which have not been initialized will throw an exception. * * The optional parameter <tt>elide_zero_values</tt> can be used to * specify whether zero values should be inserted anyway or they should be * filtered away. The default value is <tt>false</tt>, i.e., even zero * values are inserted/replaced. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Epetra_FECrsMatrix data + * structure: If the same matrix entry is inserted more than once, the + * matrix entries will be added upon calling compress() (since Epetra does + * not track values to the same entry before the final compress() is + * called), even if VectorOperation::insert is specified as argument to + * compress(). In the case you cannot make sure that matrix entries are + * only set once, initialize the matrix with a sparsity pattern to fix the + * matrix structure before inserting elements. */ void set (const size_type row, const std::vector<size_type> &col_indices, @@ -1186,14 +1216,26 @@ namespace TrilinosWrappers * * This function is able to insert new elements into the matrix as long as * compress() has not been called, so the sparsity pattern will be - * extended. When compress() is called for the first time, then this is no - * longer possible and an insertion of elements at positions which have - * not been initialized will throw an exception. + * extended. After compress() has been called for the first time or the + * matrix has been initialized from a sparsity pattern, extending the + * sparsity pattern is no longer possible and an insertion of elements at + * positions which have not been initialized will throw an exception. * * The optional parameter <tt>elide_zero_values</tt> can be used to * specify whether zero values should be inserted anyway or they should be * filtered away. The default value is <tt>false</tt>, i.e., even zero * values are inserted/replaced. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Epetra_FECrsMatrix data + * structure: If the same matrix entry is inserted more than once, the + * matrix entries will be added upon calling compress() (since Epetra does + * not track values to the same entry before the final compress() is + * called), even if VectorOperation::insert is specified as argument to + * compress(). In the case you cannot make sure that matrix entries are + * only set once, initialize the matrix with a sparsity pattern to fix the + * matrix structure before inserting elements. */ void set (const size_type row, const size_type n_cols, @@ -2509,10 +2551,11 @@ namespace TrilinosWrappers SparseMatrix::size_type SparseMatrix::n () const { + Assert(column_space_map.get() != 0, ExcInternalError()); #ifndef DEAL_II_WITH_64BIT_INDICES - return matrix->NumGlobalCols(); + return column_space_map->NumGlobalElements(); #else - return matrix->NumGlobalCols64(); + return column_space_map->NumGlobalElements64(); #endif } diff --git a/tests/trilinos/sparse_matrix_set_03.cc b/tests/trilinos/sparse_matrix_set_03.cc new file mode 100644 index 0000000000..0df5d18eed --- /dev/null +++ b/tests/trilinos/sparse_matrix_set_03.cc @@ -0,0 +1,154 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 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 setting off-processor entries of Epetra_CrsMatrix. It turns out that +// the underlying Epetra data structures actually add the entries even though +// we want to insert them only. + +#include "../tests.h" +#include <deal.II/base/utilities.h> +#include <deal.II/lac/trilinos_sparse_matrix.h> +#include <fstream> +#include <iostream> + + +void test () +{ + const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + const unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + IndexSet rows(5*n_procs); + rows.add_range(5*myid, 5*(myid+1)); + rows.compress(); + IndexSet columns(5); + columns.add_range(5*myid/n_procs, 5*(myid+1)/n_procs); + columns.compress(); + + { + TrilinosWrappers::SparseMatrix m(rows, columns, MPI_COMM_WORLD); + for (unsigned int i=5*myid; i<5*(myid+1); ++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); + + m.compress (VectorOperation::insert); + deallog << "Matrix norm set local once: " << m.frobenius_norm() << std::endl; + } + + { + TrilinosWrappers::SparseMatrix m(rows, columns, MPI_COMM_WORLD); + for (unsigned int i=5*myid; i<5*(myid+1); ++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); + + for (unsigned int i=5*myid; i<5*(myid+1); ++i) + for (unsigned int j=0; j<m.n(); ++j) + if ((i+2*j+1) % 3 == 0) + m.set (i,j, i*j+1.); + + m.compress (VectorOperation::insert); + deallog << "Matrix norm set local twice: " << m.frobenius_norm() << std::endl; + } + + { + TrilinosWrappers::SparseMatrix m(rows, columns, MPI_COMM_WORLD); + for (unsigned int i=5*myid; i<5*(myid+1); ++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); + + if (myid == 1) + m.set(1, 3, 10.); + + m.compress (VectorOperation::insert); + deallog << "Matrix norm set non-local once: " << m.frobenius_norm() << std::endl; + } + + { + TrilinosWrappers::SparseMatrix m(rows, columns, MPI_COMM_WORLD); + for (unsigned int i=5*myid; i<5*(myid+1); ++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); + + m.set(1, 3, 10.); + + m.compress (VectorOperation::insert); + deallog << "Matrix norm set non-local twice: " << m.frobenius_norm() << std::endl; + } + + { + TrilinosWrappers::SparseMatrix m(rows, columns, MPI_COMM_WORLD); + for (unsigned int i=5*myid; i<5*(myid+1); ++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); + + m.set(1, 3, 10.); + m.set(2, 3, 2*3*0.5+0.5); + + m.compress (VectorOperation::insert); + deallog << "Matrix norm set non-local twice: " << m.frobenius_norm() << std::endl; + + for (unsigned int i=5*myid; i<5*(myid+1); ++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); + + m.set(1, 3, 10.); + m.set(2, 3, 2*3*0.5+0.5); + + m.compress (VectorOperation::insert); + deallog << "Matrix norm set twice, 2nd attempt: " << m.frobenius_norm() << std::endl; + + for (unsigned int i=5*myid; i<5*(myid+1); ++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); + + m.set(1, 3, 10.); + m.set(2, 3, 2*3*0.5+0.5); + m.set(2, 3, 2*3*0.5+0.5); + + m.compress (VectorOperation::insert); + deallog << "Matrix norm set twice-twice: " << m.frobenius_norm() << std::endl; + } + + deallog << "OK" << std::endl; +} + + + +int main (int argc,char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::ofstream logfile ("output"); + deallog.attach (logfile); + deallog.depth_console (0); + deallog.threshold_double (1.e-10); + + test(); + } + else + { + test(); + } +} diff --git a/tests/trilinos/sparse_matrix_set_03.mpirun=1.output b/tests/trilinos/sparse_matrix_set_03.mpirun=1.output new file mode 100644 index 0000000000..b28c15d3ae --- /dev/null +++ b/tests/trilinos/sparse_matrix_set_03.mpirun=1.output @@ -0,0 +1,9 @@ + +DEAL::Matrix norm set local once: 9.04157 +DEAL::Matrix norm set local twice: 27.1247 +DEAL::Matrix norm set non-local once: 9.04157 +DEAL::Matrix norm set non-local twice: 13.4815 +DEAL::Matrix norm set non-local twice: 14.7817 +DEAL::Matrix norm set twice, 2nd attempt: 13.4815 +DEAL::Matrix norm set twice-twice: 13.4815 +DEAL::OK diff --git a/tests/trilinos/sparse_matrix_set_03.mpirun=2.output b/tests/trilinos/sparse_matrix_set_03.mpirun=2.output new file mode 100644 index 0000000000..f4a36ade35 --- /dev/null +++ b/tests/trilinos/sparse_matrix_set_03.mpirun=2.output @@ -0,0 +1,9 @@ + +DEAL::Matrix norm set local once: 29.9082 +DEAL::Matrix norm set local twice: 89.7246 +DEAL::Matrix norm set non-local once: 31.5357 +DEAL::Matrix norm set non-local twice: 35.9792 +DEAL::Matrix norm set non-local twice: 37.3162 +DEAL::Matrix norm set twice, 2nd attempt: 31.5357 +DEAL::Matrix norm set twice-twice: 31.5357 +DEAL::OK diff --git a/tests/trilinos/trilinos_sparsity_pattern_04.cc b/tests/trilinos/trilinos_sparsity_pattern_04.cc index e184f49af7..89d2fc6ff2 100644 --- a/tests/trilinos/trilinos_sparsity_pattern_04.cc +++ b/tests/trilinos/trilinos_sparsity_pattern_04.cc @@ -196,11 +196,6 @@ void test () 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, testing_max_num_threads()); if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)