--- /dev/null
+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)
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.
*/
}
+ 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,
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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());
+ }
+}
--- /dev/null
+
+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