From 076b455639533dd34fa600713641f1b8a9237d57 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 11 May 2015 08:24:00 -0500 Subject: [PATCH] Make MatrixOut also work for Trilinos matrices. --- doc/news/changes.h | 5 + include/deal.II/lac/matrix_out.h | 142 ++++++----- tests/lac/matrix_out_02.cc | 62 +++++ .../matrix_out_02.with_trilinos=true.output | 232 ++++++++++++++++++ 4 files changed, 380 insertions(+), 61 deletions(-) create mode 100644 tests/lac/matrix_out_02.cc create mode 100644 tests/lac/matrix_out_02.with_trilinos=true.output diff --git a/doc/news/changes.h b/doc/news/changes.h index aee4e94d26..ec481a65f8 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -493,6 +493,11 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: MatrixOut now also works with Trilinos matrices. +
    + (Wolfgang Bangerth, 2015/05/11) +
    2. +
  2. Changed: TrilinosWrappers::Vector, TrilinosWrappers::BlockVector, PETScWrappers::Vector, and PETScWrappers::BlockVector are deprecated. Either use the MPI or the deal.II version of the Vector/BlockVector. diff --git a/include/deal.II/lac/matrix_out.h b/include/deal.II/lac/matrix_out.h index 00f21e256b..7be4cbbbab 100644 --- a/include/deal.II/lac/matrix_out.h +++ b/include/deal.II/lac/matrix_out.h @@ -21,6 +21,16 @@ #include #include +#ifdef DEAL_II_WITH_TRILINOS +# include +# include +#endif + +#ifdef DEAL_II_WITH_PETSC +# include +# include +#endif + DEAL_II_NAMESPACE_OPEN /** @@ -120,12 +130,12 @@ public: * the fields of this structure for more information. * * Note that this function requires that we can extract elements of the - * matrix, which is done using the get_element() function declared below. By - * adding specializations, you can extend this class to other matrix classes - * which are not presently supported. Furthermore, we need to be able to - * extract the size of the matrix, for which we assume that the matrix type - * offers member functions m() and n(), which return the - * number of rows and columns, respectively. + * matrix, which is done using the get_element() function declared in an + * internal namespace. By adding specializations, you can extend this class + * to other matrix classes which are not presently supported. Furthermore, + * we need to be able to extract the size of the matrix, for which we assume + * that the matrix type offers member functions m() and + * n(), which return the number of rows and columns, respectively. */ template void build_patches (const Matrix &matrix, @@ -165,32 +175,6 @@ private: */ virtual std::vector get_dataset_names () const; - /** - * Return the element with given indices of a sparse matrix. - */ - template - static double get_element (const SparseMatrix &matrix, - const size_type i, - const size_type j); - - /** - * Return the element with given indices of a block sparse matrix. - */ - template - static double get_element (const BlockSparseMatrix &matrix, - const size_type i, - const size_type j); - - /** - * Return the element with given indices from any matrix type for which no - * specialization of this function was declared above. This will call - * operator() on the matrix. - */ - template - static double get_element (const Matrix &matrix, - const size_type i, - const size_type j); - /** * Get the value of the matrix at gridpoint (i,j). Depending on the * given flags, this can mean different things, for example if only absolute @@ -209,40 +193,76 @@ private: /* ---------------------- Template and inline functions ------------- */ -template -inline -double -MatrixOut::get_element (const SparseMatrix &matrix, - const size_type i, - const size_type j) +namespace internal { - return matrix.el(i,j); -} + namespace MatrixOut + { + namespace + { + /** + * Return the element with given indices of a sparse matrix. + */ + template + double get_element (const dealii::SparseMatrix &matrix, + const types::global_dof_index i, + const types::global_dof_index j) + { + return matrix.el(i,j); + } + /** + * Return the element with given indices of a block sparse matrix. + */ + template + double get_element (const dealii::BlockSparseMatrix &matrix, + const types::global_dof_index i, + const types::global_dof_index j) + { + return matrix.el(i,j); + } + +#ifdef DEAL_II_WITH_TRILINOS + /** + * Return the element with given indices of a sparse matrix. + */ + double get_element (const TrilinosWrappers::SparseMatrix &matrix, + const types::global_dof_index i, + const types::global_dof_index j) + { + return matrix.el(i,j); + } -template -inline -double -MatrixOut::get_element (const BlockSparseMatrix &matrix, - const size_type i, - const size_type j) -{ - return matrix.el(i,j); -} + /** + * Return the element with given indices of a Trilinos block sparse + * matrix. + */ + double get_element (const TrilinosWrappers::BlockSparseMatrix &matrix, + const types::global_dof_index i, + const types::global_dof_index j) + { + return matrix.el(i,j); + } +#endif -template -inline -double -MatrixOut::get_element (const Matrix &matrix, - const size_type i, - const size_type j) -{ - return matrix(i,j); + /** + * Return the element with given indices from any matrix type for which no + * specialization of this function was declared above. This will call + * operator() on the matrix. + */ + template + double get_element (const Matrix &matrix, + const types::global_dof_index i, + const types::global_dof_index j) + { + return matrix(i,j); + } + } + } } @@ -261,9 +281,9 @@ MatrixOut::get_gridpoint_value (const Matrix &matrix, if (options.block_size == 1) { if (options.show_absolute_values == true) - return std::fabs(get_element (matrix, i, j)); + return std::fabs(internal::MatrixOut::get_element (matrix, i, j)); else - return get_element (matrix, i, j); + return internal::MatrixOut::get_element (matrix, i, j); } // if blocksize greater than one, @@ -277,9 +297,9 @@ MatrixOut::get_gridpoint_value (const Matrix &matrix, col < std::min(size_type(matrix.m()), size_type((j+1)*options.block_size)); ++col, ++n_elements) if (options.show_absolute_values == true) - average += std::fabs(get_element (matrix, row, col)); + average += std::fabs(internal::MatrixOut::get_element (matrix, row, col)); else - average += get_element (matrix, row, col); + average += internal::MatrixOut::get_element (matrix, row, col); average /= n_elements; return average; } diff --git a/tests/lac/matrix_out_02.cc b/tests/lac/matrix_out_02.cc new file mode 100644 index 0000000000..c3a8f9899b --- /dev/null +++ b/tests/lac/matrix_out_02.cc @@ -0,0 +1,62 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 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. +// +// --------------------------------------------------------------------- + + +// like matrix_out.cc, but test for Trilinos matrices +// +// also test some of the other options of the MatrixOut::Options class + + +#include "../tests.h" +#include +#include +#include +#include +#include + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, numbers::invalid_unsigned_int); + + std::ofstream logfile("output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + // test for a rectangular sparse + // matrix + if (true) + { + TrilinosWrappers::SparsityPattern sparsity (4,8,7); + for (unsigned int i=0; i<4; ++i) + for (unsigned int j=0; j<8; ++j) + if (i==j+1) + sparsity.add (i,j); + sparsity.compress (); + + TrilinosWrappers::SparseMatrix sparse_matrix(sparsity); + for (unsigned int i=0; i<4; ++i) + for (unsigned int j=0; j<8; ++j) + if (i==j+1) + sparse_matrix.set(i,j, i+3*j); + + MatrixOut matrix_out; + matrix_out.build_patches (sparse_matrix, "sparse_matrix", + MatrixOut::Options (true, 1, true)); + matrix_out.write_gnuplot (logfile); + } +} diff --git a/tests/lac/matrix_out_02.with_trilinos=true.output b/tests/lac/matrix_out_02.with_trilinos=true.output new file mode 100644 index 0000000000..c0fa8d8c6e --- /dev/null +++ b/tests/lac/matrix_out_02.with_trilinos=true.output @@ -0,0 +1,232 @@ + +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 +0.00000 -1.00000 0.00000 + +1.00000 0.00000 0.00000 +1.00000 -1.00000 0.00000 + + +1.00000 0.00000 0.00000 +1.00000 -1.00000 0.00000 + +2.00000 0.00000 0.00000 +2.00000 -1.00000 0.00000 + + +2.00000 0.00000 0.00000 +2.00000 -1.00000 0.00000 + +3.00000 0.00000 0.00000 +3.00000 -1.00000 0.00000 + + +3.00000 0.00000 0.00000 +3.00000 -1.00000 0.00000 + +4.00000 0.00000 0.00000 +4.00000 -1.00000 0.00000 + + +4.00000 0.00000 0.00000 +4.00000 -1.00000 0.00000 + +5.00000 0.00000 0.00000 +5.00000 -1.00000 0.00000 + + +5.00000 0.00000 0.00000 +5.00000 -1.00000 0.00000 + +6.00000 0.00000 0.00000 +6.00000 -1.00000 0.00000 + + +6.00000 0.00000 0.00000 +6.00000 -1.00000 0.00000 + +7.00000 0.00000 0.00000 +7.00000 -1.00000 0.00000 + + +7.00000 0.00000 0.00000 +7.00000 -1.00000 0.00000 + +8.00000 0.00000 0.00000 +8.00000 -1.00000 0.00000 + + +0.00000 -1.00000 1.00000 +0.00000 -2.00000 1.00000 + +1.00000 -1.00000 1.00000 +1.00000 -2.00000 1.00000 + + +1.00000 -1.00000 0.00000 +1.00000 -2.00000 0.00000 + +2.00000 -1.00000 0.00000 +2.00000 -2.00000 0.00000 + + +2.00000 -1.00000 0.00000 +2.00000 -2.00000 0.00000 + +3.00000 -1.00000 0.00000 +3.00000 -2.00000 0.00000 + + +3.00000 -1.00000 0.00000 +3.00000 -2.00000 0.00000 + +4.00000 -1.00000 0.00000 +4.00000 -2.00000 0.00000 + + +4.00000 -1.00000 0.00000 +4.00000 -2.00000 0.00000 + +5.00000 -1.00000 0.00000 +5.00000 -2.00000 0.00000 + + +5.00000 -1.00000 0.00000 +5.00000 -2.00000 0.00000 + +6.00000 -1.00000 0.00000 +6.00000 -2.00000 0.00000 + + +6.00000 -1.00000 0.00000 +6.00000 -2.00000 0.00000 + +7.00000 -1.00000 0.00000 +7.00000 -2.00000 0.00000 + + +7.00000 -1.00000 0.00000 +7.00000 -2.00000 0.00000 + +8.00000 -1.00000 0.00000 +8.00000 -2.00000 0.00000 + + +0.00000 -2.00000 0.00000 +0.00000 -3.00000 0.00000 + +1.00000 -2.00000 0.00000 +1.00000 -3.00000 0.00000 + + +1.00000 -2.00000 5.00000 +1.00000 -3.00000 5.00000 + +2.00000 -2.00000 5.00000 +2.00000 -3.00000 5.00000 + + +2.00000 -2.00000 0.00000 +2.00000 -3.00000 0.00000 + +3.00000 -2.00000 0.00000 +3.00000 -3.00000 0.00000 + + +3.00000 -2.00000 0.00000 +3.00000 -3.00000 0.00000 + +4.00000 -2.00000 0.00000 +4.00000 -3.00000 0.00000 + + +4.00000 -2.00000 0.00000 +4.00000 -3.00000 0.00000 + +5.00000 -2.00000 0.00000 +5.00000 -3.00000 0.00000 + + +5.00000 -2.00000 0.00000 +5.00000 -3.00000 0.00000 + +6.00000 -2.00000 0.00000 +6.00000 -3.00000 0.00000 + + +6.00000 -2.00000 0.00000 +6.00000 -3.00000 0.00000 + +7.00000 -2.00000 0.00000 +7.00000 -3.00000 0.00000 + + +7.00000 -2.00000 0.00000 +7.00000 -3.00000 0.00000 + +8.00000 -2.00000 0.00000 +8.00000 -3.00000 0.00000 + + +0.00000 -3.00000 0.00000 +0.00000 -4.00000 0.00000 + +1.00000 -3.00000 0.00000 +1.00000 -4.00000 0.00000 + + +1.00000 -3.00000 0.00000 +1.00000 -4.00000 0.00000 + +2.00000 -3.00000 0.00000 +2.00000 -4.00000 0.00000 + + +2.00000 -3.00000 9.00000 +2.00000 -4.00000 9.00000 + +3.00000 -3.00000 9.00000 +3.00000 -4.00000 9.00000 + + +3.00000 -3.00000 0.00000 +3.00000 -4.00000 0.00000 + +4.00000 -3.00000 0.00000 +4.00000 -4.00000 0.00000 + + +4.00000 -3.00000 0.00000 +4.00000 -4.00000 0.00000 + +5.00000 -3.00000 0.00000 +5.00000 -4.00000 0.00000 + + +5.00000 -3.00000 0.00000 +5.00000 -4.00000 0.00000 + +6.00000 -3.00000 0.00000 +6.00000 -4.00000 0.00000 + + +6.00000 -3.00000 0.00000 +6.00000 -4.00000 0.00000 + +7.00000 -3.00000 0.00000 +7.00000 -4.00000 0.00000 + + +7.00000 -3.00000 0.00000 +7.00000 -4.00000 0.00000 + +8.00000 -3.00000 0.00000 +8.00000 -4.00000 0.00000 + + -- 2.39.5