]> https://gitweb.dealii.org/ - dealii.git/commitdiff
TrilinosWrappers::SolverDirect for Epetra data structures 17406/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 31 Jul 2024 06:47:18 +0000 (08:47 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 20 Jan 2025 14:53:38 +0000 (15:53 +0100)
doc/news/changes/minor/20240731Munch [new file with mode: 0644]
include/deal.II/lac/trilinos_solver.h
source/lac/trilinos_solver.cc
tests/trilinos/direct_solver_4.cc [new file with mode: 0644]
tests/trilinos/direct_solver_4.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20240731Munch b/doc/news/changes/minor/20240731Munch
new file mode 100644 (file)
index 0000000..8ef3e83
--- /dev/null
@@ -0,0 +1,6 @@
+New: TrilinosWrappers::SolverDirect can now also
+be used for a matrix of type Epetra_Operator and
+vectors of type Epetra_MultiVector. In addition,
+FullMatrix objects can be used as multi vectors.
+<br>
+(Peter Munch, 2024/07/31)
index ce27cac78a178cf1f91ad4420c57fa542efdb6a0..0fb267bed1cdd532bfe69c2efcfe863f18df2135 100644 (file)
@@ -625,6 +625,28 @@ namespace TrilinosWrappers
           dealii::LinearAlgebra::distributed::Vector<double>       &x,
           const dealii::LinearAlgebra::distributed::Vector<double> &b);
 
+    /**
+     * Solve the linear system <tt>Ax=b</tt> where A is an operator,
+     * and the vectors x and b are native Trilinos vector types.
+     * This function can be used when A is a LinearOperators derived
+     * from a TrilinosPayload.
+     */
+    void
+    solve(const Epetra_Operator    &A,
+          Epetra_MultiVector       &x,
+          const Epetra_MultiVector &b);
+
+    /**
+     * Solve the linear system <tt>AX=B</tt> where A is an operator,
+     * and the X and B are FullMatrix types. The matrices are
+     * converted to multi vector and the resulting problem is solved
+     * in a batched way.
+     */
+    void
+    solve(const SparseMatrix       &sparse_matrix,
+          FullMatrix<double>       &solution,
+          const FullMatrix<double> &rhs);
+
     /**
      * Access to object that controls convergence.
      */
index d99cf9287f9568ab325641a5b10cb1ef9a1f2daf..0975c8f04741a13a0ab68da64ab0a2122638842d 100644 (file)
@@ -844,6 +844,69 @@ namespace TrilinosWrappers
   }
 
 
+  void
+  SolverDirect::solve(const Epetra_Operator    &A,
+                      Epetra_MultiVector       &x,
+                      const Epetra_MultiVector &b)
+  {
+    // We need an Epetra_LinearProblem object to let the Amesos solver know
+    // about the matrix and vectors.
+    linear_problem =
+      std::make_unique<Epetra_LinearProblem>(const_cast<Epetra_Operator *>(&A),
+                                             &x,
+                                             const_cast<Epetra_MultiVector *>(
+                                               &b));
+
+    do_solve();
+  }
+
+
+  void
+  SolverDirect::solve(const SparseMatrix       &sparse_matrix,
+                      FullMatrix<double>       &solution,
+                      const FullMatrix<double> &rhs)
+  {
+    Assert(sparse_matrix.m() == sparse_matrix.n(), ExcInternalError());
+    Assert(rhs.m() == sparse_matrix.m(), ExcInternalError());
+    Assert(rhs.m() == solution.m(), ExcInternalError());
+    Assert(rhs.n() == solution.n(), ExcInternalError());
+
+    const unsigned int m = rhs.m();
+    const unsigned int n = rhs.n();
+
+    FullMatrix<double> rhs_t(n, m);
+    FullMatrix<double> solution_t(n, m);
+
+    rhs_t.copy_transposed(rhs);
+    solution_t.copy_transposed(solution);
+
+    std::vector<double *> rhs_ptrs(n);
+    std::vector<double *> sultion_ptrs(n);
+
+    for (unsigned int i = 0; i < n; ++i)
+      {
+        rhs_ptrs[i]     = &rhs_t[i][0];
+        sultion_ptrs[i] = &solution_t[i][0];
+      }
+
+    const Epetra_CrsMatrix &mat = sparse_matrix.trilinos_matrix();
+
+    Epetra_MultiVector trilinos_dst(View,
+                                    mat.OperatorRangeMap(),
+                                    sultion_ptrs.data(),
+                                    sultion_ptrs.size());
+    Epetra_MultiVector trilinos_src(View,
+                                    mat.OperatorDomainMap(),
+                                    rhs_ptrs.data(),
+                                    rhs_ptrs.size());
+
+    this->initialize(sparse_matrix);
+    this->solve(mat, trilinos_dst, trilinos_src);
+
+    solution.copy_transposed(solution_t);
+  }
+
+
   void
   SolverDirect::solve(const SparseMatrix &A,
                       MPI::Vector        &x,
diff --git a/tests/trilinos/direct_solver_4.cc b/tests/trilinos/direct_solver_4.cc
new file mode 100644 (file)
index 0000000..1fe009f
--- /dev/null
@@ -0,0 +1,112 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// tests Trilinos direct solvers for multi vectors
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+
+  const unsigned int n = 5;
+  const unsigned int m = 3;
+
+  // initialize the matrix
+  FullMatrix<double> A(n, n);
+
+  for (unsigned int i = 0; i < n; ++i)
+    A[i][i] = 2.0;
+
+  for (unsigned int i = 0; i < n - 1; ++i)
+    {
+      A[i][i + 1] = -1.0;
+      A[i + 1][i] = -1.0;
+    }
+
+  // initialize the right-hand-side vector
+  FullMatrix<double> B(n, m);
+
+  for (unsigned int i = 0; i < m; ++i)
+    B[i][i] = 1;
+
+  // solve with direct solver
+  {
+    auto A_copy = A;
+    A_copy.gauss_jordan();
+
+    FullMatrix<double> X(n, m);
+    A_copy.mmult(X, B);
+
+    X.print(deallog.get_file_stream());
+  }
+
+  deallog << std::endl;
+
+  // solve with sparse direct solver
+  {
+    std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
+      entries;
+
+    for (unsigned int i = 0; i < A.m(); ++i)
+      for (unsigned int j = 0; j < A.n(); ++j)
+        if (A[i][j] != 0)
+          entries.emplace_back(i, j);
+
+    SparsityPattern sp(n, n);
+    sp.add_entries(entries);
+    sp.compress();
+
+    TrilinosWrappers::SparseMatrix A_sparse;
+    A_sparse.reinit(sp);
+
+    for (unsigned int i = 0; i < A.m(); ++i)
+      for (unsigned int j = 0; j < A.n(); ++j)
+        if (A[i][j] != 0)
+          A_sparse.set(i, j, A[i][j]);
+
+    TrilinosWrappers::SolverDirect solver;
+
+    FullMatrix<double> X(n, m);
+    solver.solve(A_sparse, X, B);
+    X.print(deallog.get_file_stream());
+  }
+}
diff --git a/tests/trilinos/direct_solver_4.output b/tests/trilinos/direct_solver_4.output
new file mode 100644 (file)
index 0000000..004f556
--- /dev/null
@@ -0,0 +1,14 @@
+
+0.83 0.67 0.50 
+0.67 1.3  1.0  
+0.50 1.0  1.5  
+0.33 0.67 1.0  
+0.17 0.33 0.50 
+DEAL::
+DEAL::Starting value 0.00000
+DEAL::Convergence step 0 value 0.00000
+0.83 0.67 0.50 
+0.67 1.3  1.0  
+0.50 1.0  1.5  
+0.33 0.67 1.0  
+0.17 0.33 0.50 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.