}
+/**
+ * @relates LinearOperator
+ *
+ * Return a LinearOperator that is the identity of the vector space @p Range.
+ *
+ * The function takes a LinearOperator @p op and uses its range initializer
+ * to create an identiy operator. In contrast to the function above, this
+ * function also ensures that the underlying Payload matches that of the
+ * input @p op.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+identity_operator(const LinearOperator<Range, Domain, Payload> &op)
+{
+ const LinearOperator<Range, Domain, Payload> id_op
+ = identity_operator<Range,Payload>(op.reinit_range_vector);
+ LinearOperator<Range, Domain, Payload> return_op (
+ op.identity_payload()
+ );
+
+ return_op.reinit_range_vector = id_op.reinit_range_vector;
+ return_op.reinit_domain_vector = id_op.reinit_domain_vector;
+ return_op.vmult = id_op.vmult;
+ return_op.vmult_add = id_op.vmult_add;
+ return_op.Tvmult = id_op.Tvmult;
+ return_op.Tvmult_add = id_op.Tvmult_add;
+
+ return return_op;
+}
+
+
/**
* @relates LinearOperator
*
LinearOperator<Range, Domain, Payload>
null_operator(const LinearOperator<Range, Domain, Payload> &op)
{
- auto return_op = op;
+ LinearOperator<Range, Domain, Payload> return_op (
+ op.null_payload()
+ );
return_op.is_null_operator = true;
{ }
+ /**
+ * Returns a payload configured for identity operations
+ */
+ EmptyPayload
+ identity_payload () const
+ {
+ return *this;
+ }
+
+
+ /**
+ * Returns a payload configured for null operations
+ */
+ EmptyPayload
+ null_payload () const
+ {
+ return *this;
+ }
+
+
/**
* Returns a payload configured for transpose operations
*/
//@{
- /**
- * @relates LinearOperator
- *
- * Return a LinearOperator that is the identity of the vector space @p Range.
- *
- * This function is the equivalent of the dealii::identity_operator, but
- * ensures full compatibility with Trilinos operations by preselecting the
- * appropriate template parameters.
- *
- * @author Jean-Paul Pelteret, 2016
- *
- * @ingroup TrilinosWrappers
- */
- template <typename Range>
- inline LinearOperator<Range, Range, TrilinosWrappers::internal::LinearOperator::TrilinosPayload>
- identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
- {
- typedef TrilinosWrappers::internal::LinearOperator::TrilinosPayload Payload;
- return dealii::identity_operator<Range, Payload>(reinit_vector);
- }
-
-
/**
* @relates LinearOperator
*
*/
virtual ~TrilinosPayload() {}
+ /**
+ * Returns a payload configured for identity operations
+ */
+ TrilinosPayload identity_payload () const;
+
+ /**
+ * Returns a payload configured for null operations
+ */
+ TrilinosPayload null_payload () const;
+
/**
* Returns a payload configured for transpose operations
*/
TrilinosPayload transpose_payload () const;
-
/**
* Returns a payload configured for inverse operations
*
+ TrilinosPayload
+ TrilinosPayload::identity_payload () const
+ {
+ TrilinosPayload return_op (*this);
+
+ return_op.vmult = [](Range &tril_dst,
+ const Range &tril_src)
+ {
+ tril_dst = tril_src;
+ };
+
+ return_op.Tvmult = [](Range &tril_dst,
+ const Range &tril_src)
+ {
+ tril_dst = tril_src;
+ };
+
+ return_op.inv_vmult = [](Range &tril_dst,
+ const Range &tril_src)
+ {
+ tril_dst = tril_src;
+ };
+
+ return_op.inv_Tvmult = [](Range &tril_dst,
+ const Range &tril_src)
+ {
+ tril_dst = tril_src;
+ };
+
+ return return_op;
+ }
+
+
+
+ TrilinosPayload
+ TrilinosPayload::null_payload () const
+ {
+ TrilinosPayload return_op (*this);
+
+ return_op.vmult = [](Range &tril_dst,
+ const Domain &)
+ {
+ tril_dst.Scale(0);
+ };
+
+ return_op.Tvmult = [](Domain &tril_dst,
+ const Range &)
+ {
+ tril_dst.Scale(0);
+ };
+
+ return_op.inv_vmult = [](Domain &tril_dst,
+ const Range &)
+ {
+ AssertThrow(false, ExcMessage("Cannot compute inverse of null operator"));
+ tril_dst.Scale(0);
+ };
+
+ return_op.inv_Tvmult = [](Range &tril_dst,
+ const Domain &)
+ {
+ AssertThrow(false, ExcMessage("Cannot compute inverse of null operator"));
+ tril_dst.Scale(0);
+ };
+
+ return return_op;
+ }
+
+
+
TrilinosPayload
TrilinosPayload::transpose_payload () const
{
--- /dev/null
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2017 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 the identity_operator LinearOperator and its payload initialization
+
+
+#include "../tests.h"
+#include "../testmatrix.h"
+
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/linear_operator_tools.h>
+
+using namespace dealii;
+
+int main (int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ initlog();
+
+ // Need to get the linear operator to do some nonsense work to verify
+ // that it can indeed be used as expected
+ const unsigned int size = 32;
+ unsigned int dim = (size-1)*(size-1);
+
+ // Make matrix
+ FDMatrix testproblem(size, size);
+
+ {
+ SparsityPattern sp (dim, dim);
+ testproblem.five_point_structure(sp);
+ sp.compress();
+
+ SparseMatrix<double> A;
+ A.reinit(sp);
+ testproblem.five_point(A);
+
+ Vector<double> f(dim);
+ Vector<double> u(dim);
+
+ const auto lo_A = linear_operator<Vector<double> >(A);
+ const auto lo_id = identity_operator<Vector<double> >(lo_A.reinit_range_vector);
+ const auto lo_A_plus_id = lo_A + lo_id;
+
+ u = lo_A_plus_id * f;
+ }
+
+ {
+ DynamicSparsityPattern csp (dim, dim);
+ testproblem.five_point_structure(csp);
+
+ TrilinosWrappers::SparseMatrix A;
+ A.reinit(csp);
+ testproblem.five_point(A);
+
+ TrilinosWrappers::Vector f(dim);
+ TrilinosWrappers::Vector u(dim);
+
+ A.compress (VectorOperation::insert);
+ f.compress (VectorOperation::insert);
+ u.compress (VectorOperation::insert);
+
+ const auto lo_A = linear_operator<TrilinosWrappers::Vector>(A);
+ const auto lo_id_1 = identity_operator<TrilinosWrappers::Vector>(lo_A.reinit_range_vector);
+ const auto lo_id_2 = identity_operator<TrilinosWrappers::Vector>(lo_A);
+ const auto lo_A_plus_id_1 = lo_A + lo_id_1; // Not a good idea. See below.
+ const auto lo_A_plus_id_2 = lo_A + lo_id_2;
+
+ // u = lo_A_plus_id_1 * f; // This is not allowed -- different payloads
+ u = lo_A_plus_id_2 * f;
+ }
+
+ deallog << "OK" << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2017 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 the null_operator LinearOperator and its payload initialization
+
+
+#include "../tests.h"
+#include "../testmatrix.h"
+
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/linear_operator_tools.h>
+
+using namespace dealii;
+
+int main (int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ initlog();
+
+ // Need to get the linear operator to do some nonsense work to verify
+ // that it can indeed be used as expected
+ const unsigned int size = 32;
+ unsigned int dim = (size-1)*(size-1);
+
+ // Make matrix
+ FDMatrix testproblem(size, size);
+
+ {
+ SparsityPattern sp (dim, dim);
+ testproblem.five_point_structure(sp);
+ sp.compress();
+
+ SparseMatrix<double> A;
+ A.reinit(sp);
+ testproblem.five_point(A);
+
+ Vector<double> f(dim);
+ Vector<double> u(dim);
+
+ const auto lo_A = linear_operator<Vector<double> >(A);
+ const auto lo_null = null_operator(lo_A);
+ const auto lo_A_plus_null = lo_A + lo_null;
+
+ u = lo_A_plus_null * f;
+ }
+
+ {
+ DynamicSparsityPattern csp (dim, dim);
+ testproblem.five_point_structure(csp);
+
+ TrilinosWrappers::SparseMatrix A;
+ A.reinit(csp);
+ testproblem.five_point(A);
+
+ TrilinosWrappers::Vector f(dim);
+ TrilinosWrappers::Vector u(dim);
+
+ A.compress (VectorOperation::insert);
+ f.compress (VectorOperation::insert);
+ u.compress (VectorOperation::insert);
+
+ const auto lo_A = linear_operator<TrilinosWrappers::Vector>(A);
+ const auto lo_null = null_operator(lo_A);
+ const auto lo_A_plus_null = lo_A + lo_null;
+
+ u = lo_A_plus_null * f;
+ }
+
+ deallog << "OK" << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::OK
projects / dealii.git / commitdiff