MFEM  v4.6.0
Finite element discretization library
mumps.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_MUMPS
13 #define MFEM_MUMPS
14 
15 #include "../config/config.hpp"
16 
17 #ifdef MFEM_USE_MUMPS
18 #ifdef MFEM_USE_MPI
19 
20 #include "operator.hpp"
21 #include "hypre.hpp"
22 #include <mpi.h>
23 
24 #include "dmumps_c.h"
25 
26 namespace mfem
27 {
28 
29 /**
30  * @brief MUMPS: A Parallel Sparse Direct Solver
31  *
32  * Interface for the distributed MUMPS solver
33  */
34 class MUMPSSolver : public Solver
35 {
36 public:
37  enum MatType
38  {
42  };
43 
45  {
46  AUTOMATIC = 0,
47  AMD,
48  AMF,
54  };
55 
56  /**
57  * @brief Constructor with MPI_Comm parameter.
58  */
59  MUMPSSolver(MPI_Comm comm_);
60 
61  /**
62  * @brief Constructor with a HypreParMatrix Operator.
63  */
64  MUMPSSolver(const Operator &op);
65 
66  /**
67  * @brief Set the Operator and perform factorization
68  *
69  * @a op needs to be of type HypreParMatrix.
70  *
71  * @param op Operator used in factorization and solve
72  */
73  void SetOperator(const Operator &op);
74 
75  /**
76  * @brief Solve y = Op^{-1} x.
77  *
78  * @param x RHS vector
79  * @param y Solution vector
80  */
81  void Mult(const Vector &x, Vector &y) const;
82  void ArrayMult(const Array<const Vector *> &X, Array<Vector *> &Y) const;
83 
84  /**
85  * @brief Transpose Solve y = Op^{-T} x.
86  *
87  * @param x RHS vector
88  * @param y Solution vector
89  */
90  void MultTranspose(const Vector &x, Vector &y) const;
92  Array<Vector *> &Y) const;
93 
94  /**
95  * @brief Set the error print level for MUMPS
96  *
97  * @param print_lvl Print level
98  *
99  * @note This method has to be called before SetOperator
100  */
101  void SetPrintLevel(int print_lvl);
102 
103  /**
104  * @brief Set the matrix type
105  *
106  * Supported matrix types: General, symmetric indefinite and
107  * symmetric positive definite
108  *
109  * @param mtype Matrix type
110  *
111  * @note This method has to be called before SetOperator
112  */
113  void SetMatrixSymType(MatType mtype);
114 
115  /**
116  * @brief Set the reordering strategy
117  *
118  * Supported reorderings are: AUTOMATIC, AMD, AMF, PORD, METIS, PARMETIS,
119  * SCOTCH, and PTSCOTCH
120  *
121  * @param method Reordering method
122  *
123  * @note This method has to be called before SetOperator
124  */
126 
127  /**
128  * @brief Set the flag controlling reuse of the symbolic factorization
129  * for multiple operators
130  *
131  * @param reuse Flag to reuse symbolic factorization
132  *
133  * @note This method has to be called before repeated calls to SetOperator
134  */
135  void SetReorderingReuse(bool reuse);
136 
137  /**
138  * @brief Set the tolerance for activating block low-rank (BLR) approximate
139  * factorization
140  *
141  * @param tol Tolerance
142  *
143  * @note This method has to be called before SetOperator
144  */
145 #if MFEM_MUMPS_VERSION >= 510
146  void SetBLRTol(double tol);
147 #endif
148 
149  // Destructor
150  ~MUMPSSolver();
151 
152 private:
153  // MPI communicator
154  MPI_Comm comm;
155 
156  // Number of procs
157  int numProcs;
158 
159  // MPI rank
160  int myid;
161 
162  // Parameter controlling the matrix type
163  MatType mat_type;
164 
165  // Parameter controlling the printing level
166  int print_level;
167 
168  // Parameter controlling the reordering strategy
169  ReorderingStrategy reorder_method;
170 
171  // Parameter controlling whether or not to reuse the symbolic factorization
172  // for multiple calls to SetOperator
173  bool reorder_reuse;
174 
175 #if MFEM_MUMPS_VERSION >= 510
176  // Parameter controlling the Block Low-Rank (BLR) feature in MUMPS
177  double blr_tol;
178 #endif
179 
180  // Local row offsets
181  int row_start;
182 
183  // MUMPS object
184  DMUMPS_STRUC_C *id;
185 
186  // Method for initialization
187  void Init(MPI_Comm comm_);
188 
189  // Method for setting MUMPS internal parameters
190  void SetParameters();
191 
192  // Method for configuring storage for distributed/centralized RHS and
193  // solution
194  void InitRhsSol(int nrhs) const;
195 
196 #if MFEM_MUMPS_VERSION >= 530
197  // Row offests array on all procs
198  Array<int> row_starts;
199 
200  // Row maps and storage for distributed RHS and solution
201  int *irhs_loc, *isol_loc;
202  mutable double *rhs_loc, *sol_loc;
203 
204  // These two methods are needed to distribute the local solution
205  // vectors returned by MUMPS to the original MFEM parallel partition
206  int GetRowRank(int i, const Array<int> &row_starts_) const;
207  void RedistributeSol(const int *rmap, const double *x, const int lx_loc,
208  Array<Vector *> &Y) const;
209 #else
210  // Arrays needed for MPI_Gatherv and MPI_Scatterv
211  int *recv_counts, *displs;
212  mutable double *rhs_glob;
213 #endif
214 }; // mfem::MUMPSSolver class
215 
216 } // namespace mfem
217 
218 #endif // MFEM_USE_MPI
219 #endif // MFEM_USE_MUMPS
220 #endif // MFEM_MUMPS
void ArrayMultTranspose(const Array< const Vector *> &X, Array< Vector *> &Y) const
Action of the transpose operator on a matrix: Y=A^t(X).
Definition: mumps.cpp:399
void SetOperator(const Operator &op)
Set the Operator and perform factorization.
Definition: mumps.cpp:98
void SetPrintLevel(int print_lvl)
Set the error print level for MUMPS.
Definition: mumps.cpp:410
MUMPSSolver(MPI_Comm comm_)
Constructor with MPI_Comm parameter.
Definition: mumps.cpp:40
void SetMatrixSymType(MatType mtype)
Set the matrix type.
Definition: mumps.cpp:415
void Mult(const Vector &x, Vector &y) const
Solve y = Op^{-1} x.
Definition: mumps.cpp:328
void ArrayMult(const Array< const Vector *> &X, Array< Vector *> &Y) const
Operator application on a matrix: Y=A(X).
Definition: mumps.cpp:337
void SetReorderingStrategy(ReorderingStrategy method)
Set the reordering strategy.
Definition: mumps.cpp:420
void SetBLRTol(double tol)
Set the tolerance for activating block low-rank (BLR) approximate factorization.
Definition: mumps.cpp:431
void SetReorderingReuse(bool reuse)
Set the flag controlling reuse of the symbolic factorization for multiple operators.
Definition: mumps.cpp:425
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:54
void MultTranspose(const Vector &x, Vector &y) const
Transpose Solve y = Op^{-T} x.
Definition: mumps.cpp:389
Vector data type.
Definition: vector.hpp:58
Base class for solvers.
Definition: operator.hpp:682
Abstract operator.
Definition: operator.hpp:24
MUMPS: A Parallel Sparse Direct Solver.
Definition: mumps.hpp:34