From cee65e3fca6f68eae43805021fc53a7c3cad3a61 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 31 Jul 2024 08:47:18 +0200 Subject: [PATCH] TrilinosWrappers::SolverDirect for Epetra data structures --- doc/news/changes/minor/20240731Munch | 6 ++ include/deal.II/lac/trilinos_solver.h | 22 +++++ source/lac/trilinos_solver.cc | 63 +++++++++++++++ tests/trilinos/direct_solver_4.cc | 112 ++++++++++++++++++++++++++ tests/trilinos/direct_solver_4.output | 14 ++++ 5 files changed, 217 insertions(+) create mode 100644 doc/news/changes/minor/20240731Munch create mode 100644 tests/trilinos/direct_solver_4.cc create mode 100644 tests/trilinos/direct_solver_4.output diff --git a/doc/news/changes/minor/20240731Munch b/doc/news/changes/minor/20240731Munch new file mode 100644 index 0000000000..8ef3e8389e --- /dev/null +++ b/doc/news/changes/minor/20240731Munch @@ -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. +
+(Peter Munch, 2024/07/31) diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h index ce27cac78a..0fb267bed1 100644 --- a/include/deal.II/lac/trilinos_solver.h +++ b/include/deal.II/lac/trilinos_solver.h @@ -625,6 +625,28 @@ namespace TrilinosWrappers dealii::LinearAlgebra::distributed::Vector &x, const dealii::LinearAlgebra::distributed::Vector &b); + /** + * Solve the linear system Ax=b 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 AX=B 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 &solution, + const FullMatrix &rhs); + /** * Access to object that controls convergence. */ diff --git a/source/lac/trilinos_solver.cc b/source/lac/trilinos_solver.cc index d99cf9287f..0975c8f047 100644 --- a/source/lac/trilinos_solver.cc +++ b/source/lac/trilinos_solver.cc @@ -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(const_cast(&A), + &x, + const_cast( + &b)); + + do_solve(); + } + + + void + SolverDirect::solve(const SparseMatrix &sparse_matrix, + FullMatrix &solution, + const FullMatrix &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 rhs_t(n, m); + FullMatrix solution_t(n, m); + + rhs_t.copy_transposed(rhs); + solution_t.copy_transposed(solution); + + std::vector rhs_ptrs(n); + std::vector 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 index 0000000000..1fe009f226 --- /dev/null +++ b/tests/trilinos/direct_solver_4.cc @@ -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 + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#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 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 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 X(n, m); + A_copy.mmult(X, B); + + X.print(deallog.get_file_stream()); + } + + deallog << std::endl; + + // solve with sparse direct solver + { + std::vector> + 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 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 index 0000000000..004f5563a8 --- /dev/null +++ b/tests/trilinos/direct_solver_4.output @@ -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 -- 2.39.5