MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
handle.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 #include "handle.hpp"
13 #include "sparsemat.hpp"
14 #ifdef MFEM_USE_MPI
15 #include "petsc.hpp"
16 #endif
17 
18 // Make sure that hypre and PETSc use the same size indices.
19 #if defined(MFEM_USE_MPI) && defined(MFEM_USE_PETSC)
20 #if (defined(HYPRE_BIGINT) && !defined(PETSC_USE_64BIT_INDICES)) || \
21  (!defined(HYPRE_BIGINT) && defined(PETSC_USE_64BIT_INDICES))
22 #error HYPRE and PETSC do not use the same size integers!
23 #endif
24 #endif
25 
26 namespace mfem
27 {
28 
30  "Operator::Type is not supported: type_id = ";
31 
33 {
34  switch (tid)
35  {
36  case Operator::ANY_TYPE: break;
37  case Operator::MFEM_SPARSEMAT: break;
39 #ifdef MFEM_USE_MPI
40  break;
41 #else
42  MFEM_ABORT("cannot use HYPRE parallel matrix format: "
43  "MFEM is not built with HYPRE support");
44 #endif
47 #ifdef MFEM_USE_PETSC
48  break;
49 #else
50  MFEM_ABORT("cannot use PETSc matrix formats: "
51  "MFEM is not built with PETSc support");
52 #endif
53  default:
54  MFEM_ABORT("invalid Operator::Type, type_id = " << (int)type_id);
55  }
56  return tid;
57 }
58 
59 #ifdef MFEM_USE_MPI
60 void OperatorHandle::MakeSquareBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_size,
61  HYPRE_BigInt *row_starts,
62  SparseMatrix *diag)
63 {
64  if (own_oper) { delete oper; }
65 
66  switch (type_id)
67  {
68  case Operator::ANY_TYPE: // --> MFEM_SPARSEMAT
70  // As a parallel Operator, the SparseMatrix simply represents the local
71  // Operator, without the need of any communication.
72  pSet(diag, false);
73  return;
75  oper = new HypreParMatrix(comm, glob_size, row_starts, diag);
76  break;
77 #ifdef MFEM_USE_PETSC
80  // Assuming that PetscInt is the same size as HYPRE_Int, checked above.
81  oper = new PetscParMatrix(comm, glob_size, (PetscInt*)row_starts, diag,
82  type_id);
83  break;
84 #endif
85  default: MFEM_ABORT(not_supported_msg << type_id);
86  }
87  own_oper = true;
88 }
89 
91 MakeRectangularBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_num_rows,
92  HYPRE_BigInt glob_num_cols, HYPRE_BigInt *row_starts,
93  HYPRE_BigInt *col_starts, SparseMatrix *diag)
94 {
95  if (own_oper) { delete oper; }
96 
97  switch (type_id)
98  {
99  case Operator::ANY_TYPE: // --> MFEM_SPARSEMAT
101  // As a parallel Operator, the SparseMatrix simply represents the local
102  // Operator, without the need of any communication.
103  pSet(diag, false);
104  return;
106  oper = new HypreParMatrix(comm, glob_num_rows, glob_num_cols,
107  row_starts, col_starts, diag);
108  break;
109 #ifdef MFEM_USE_PETSC
112  // Assuming that PetscInt is the same size as HYPRE_Int, checked above.
113  oper = new PetscParMatrix(comm, glob_num_rows, glob_num_cols,
114  (PetscInt*)row_starts, (PetscInt*)col_starts, diag, type_id);
115  break;
116 #endif
117  default: MFEM_ABORT(not_supported_msg << type_id);
118  }
119  own_oper = true;
120 }
121 #endif // MFEM_USE_MPI
122 
124 {
125  if (A.Type() != Operator::ANY_TYPE)
126  {
127  MFEM_VERIFY(A.Type() == P.Type(), "type mismatch in A and P");
128  }
129  Clear();
130  switch (A.Type())
131  {
132  case Operator::ANY_TYPE:
133  pSet(new RAPOperator(*P.Ptr(), *A.Ptr(), *P.Ptr()));
134  break;
136  {
138  SparseMatrix *RA = mfem::Mult(*R, *A.As<SparseMatrix>());
139  delete R;
140  pSet(mfem::Mult(*RA, *P.As<SparseMatrix>()));
141  delete RA;
142  break;
143  }
144 #ifdef MFEM_USE_MPI
147  break;
148 #ifdef MFEM_USE_PETSC
151  {
153  break;
154  }
155 #endif
156 #endif
157  default: MFEM_ABORT(not_supported_msg << A.Type());
158  }
159 }
160 
162  OperatorHandle &P)
163 {
164  if (A.Type() != Operator::ANY_TYPE)
165  {
166  MFEM_VERIFY(A.Type() == Rt.Type(), "type mismatch in A and Rt");
167  MFEM_VERIFY(A.Type() == P.Type(), "type mismatch in A and P");
168  }
169  Clear();
170  switch (A.Type())
171  {
172  case Operator::ANY_TYPE:
173  pSet(new RAPOperator(*Rt.Ptr(), *A.Ptr(), *P.Ptr()));
174  break;
176  {
178  *P.As<SparseMatrix>()));
179  break;
180  }
181 #ifdef MFEM_USE_MPI
184  P.As<HypreParMatrix>()));
185  break;
186 #ifdef MFEM_USE_PETSC
189  {
191  P.As<PetscParMatrix>()));
192  break;
193  }
194 #endif
195 #endif
196  default: MFEM_ABORT(not_supported_msg << A.Type());
197  }
198 }
199 
201 {
202  if (own_oper) { delete oper; }
203  if (Type() == A.Type() || Type() == Operator::ANY_TYPE)
204  {
205  oper = A.Ptr();
206  own_oper = false;
207  return;
208  }
209  oper = NULL;
210  switch (Type()) // target type id
211  {
213  {
214  oper = A.Is<SparseMatrix>();
215  break;
216  }
218  {
219 #ifdef MFEM_USE_MPI
220  oper = A.Is<HypreParMatrix>();
221 #endif
222  break;
223  }
226  {
227  switch (A.Type()) // source type id
228  {
230 #ifdef MFEM_USE_PETSC
231  oper = new PetscParMatrix(A.As<HypreParMatrix>(), Type());
232 #endif
233  break;
234  default: break;
235  }
236 #ifdef MFEM_USE_PETSC
237  if (!oper)
238  {
239  PetscParMatrix *pA = A.Is<PetscParMatrix>();
240  if (pA->GetType() == Type()) { oper = pA; }
241  }
242 #endif
243  break;
244  }
245  default: break;
246  }
247  MFEM_VERIFY(oper != NULL, "conversion from type id = " << A.Type()
248  << " to type id = " << Type() << " is not supported");
249  own_oper = true;
250 }
251 
253  const Array<int> &ess_dof_list)
254 {
255  Clear();
256  switch (A.Type())
257  {
258  case Operator::ANY_TYPE:
259  {
260  bool own_A = A.OwnsOperator();
261  A.SetOperatorOwner(false);
262  A.Reset(new ConstrainedOperator(A.Ptr(), ess_dof_list, own_A));
263  // Keep this object empty - this will be OK if this object is only
264  // used as the A_e parameter in a call to A.EliminateBC().
265  break;
266  }
268  {
269  const Matrix::DiagonalPolicy preserve_diag = Matrix::DIAG_KEEP;
270  SparseMatrix *sA = A.As<SparseMatrix>();
271  SparseMatrix *Ae = new SparseMatrix(sA->Height());
272  for (int i = 0; i < ess_dof_list.Size(); i++)
273  {
274  sA->EliminateRowCol(ess_dof_list[i], *Ae, preserve_diag);
275  }
276  Ae->Finalize();
277  pSet(Ae);
278  break;
279  }
281  {
282 #ifdef MFEM_USE_MPI
283  pSet(A.As<HypreParMatrix>()->EliminateRowsCols(ess_dof_list));
284 #else
285  MFEM_ABORT("type id = Hypre_ParCSR requires MFEM_USE_MPI");
286 #endif
287  break;
288  }
291  {
292 #ifdef MFEM_USE_PETSC
293  pSet(A.As<PetscParMatrix>()->EliminateRowsCols(ess_dof_list));
294 #else
295  MFEM_ABORT("type id = Operator::PETSC_* requires MFEM_USE_PETSC");
296 #endif
297  break;
298  }
299  default: MFEM_ABORT(not_supported_msg << A.Type());
300  }
301 }
302 
303 void OperatorHandle::EliminateRows(const Array<int> &ess_dof_list)
304 {
305  switch (Type())
306  {
308  {
309 #ifdef MFEM_USE_MPI
310  this->As<HypreParMatrix>()->EliminateRows(ess_dof_list);
311 #else
312  MFEM_ABORT("type id = Hypre_ParCSR requires MFEM_USE_MPI");
313 #endif
314  break;
315  }
316  default:
317  MFEM_ABORT(not_supported_msg << Type());
318  }
319 }
320 
321 void OperatorHandle::EliminateCols(const Array<int> &ess_dof_list)
322 {
323  switch (Type())
324  {
326  {
327 #ifdef MFEM_USE_MPI
328  auto Ae = this->As<HypreParMatrix>()->EliminateCols(ess_dof_list);
329  delete Ae;
330 #else
331  MFEM_ABORT("type id = Hypre_ParCSR requires MFEM_USE_MPI");
332 #endif
333  break;
334  }
335  default:
336  MFEM_ABORT(not_supported_msg << Type());
337  }
338 }
339 
341  const Array<int> &ess_dof_list,
342  const Vector &X, Vector &B) const
343 {
344  switch (Type())
345  {
346  case Operator::ANY_TYPE:
347  {
348  ConstrainedOperator *A = Is<ConstrainedOperator>();
349  MFEM_VERIFY(A != NULL, "EliminateRowsCols() is not called");
350  A->EliminateRHS(X, B);
351  break;
352  }
354  {
355  A_e.As<SparseMatrix>()->AddMult(X, B, -1.);
356  As<SparseMatrix>()->PartMult(ess_dof_list, X, B);
357  break;
358  }
360  {
361 #ifdef MFEM_USE_MPI
362  mfem::EliminateBC(*As<HypreParMatrix>(), *A_e.As<HypreParMatrix>(),
363  ess_dof_list, X, B);
364 #else
365  MFEM_ABORT("type id = Hypre_ParCSR requires MFEM_USE_MPI");
366 #endif
367  break;
368  }
371  {
372 #ifdef MFEM_USE_PETSC
373  mfem::EliminateBC(*As<PetscParMatrix>(), *A_e.As<PetscParMatrix>(),
374  ess_dof_list, X, B);
375 #else
376  MFEM_ABORT("type id = Operator::PETSC_* requires MFEM_USE_PETSC");
377 #endif
378  break;
379  }
380  default: MFEM_ABORT(not_supported_msg << Type());
381  }
382 }
383 
384 } // namespace mfem
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:2052
void ConvertFrom(OperatorHandle &A)
Convert the given OperatorHandle A to the currently set type id.
Definition: handle.cpp:200
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:99
void EliminateBC(const HypreParMatrix &A, const HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Eliminate essential BC specified by ess_dof_list from the solution X to the r.h.s. B.
Definition: hypre.cpp:2857
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:2485
HYPRE_Int PetscInt
Definition: petsc.hpp:47
Wrapper for PETSc&#39;s matrix class.
Definition: petsc.hpp:307
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:94
void EliminateBC(const OperatorHandle &A_e, const Array< int > &ess_dof_list, const Vector &X, Vector &B) const
Eliminate essential dofs from the solution X into the r.h.s. B.
Definition: handle.cpp:340
Type GetType() const
Definition: petsc.cpp:2227
void EliminateRows(const Array< int > &ess_dof_list)
Eliminate the rows corresponding to the essential dofs ess_dof_list.
Definition: handle.cpp:303
void EliminateRowsCols(const Array< int > &rows_cols, const PetscParVector &X, PetscParVector &B, double diag=1.)
Eliminate rows and columns from the matrix, and rows from the vector B. Modify B with the BC values i...
Definition: petsc.cpp:2171
Operator * oper
Definition: handle.hpp:38
void EliminateCols(const Array< int > &ess_dof_list)
Eliminate columns corresponding to the essential dofs ess_dof_list.
Definition: handle.cpp:321
static const char not_supported_msg[]
Definition: handle.hpp:36
ID for class SparseMatrix.
Definition: operator.hpp:255
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:254
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition: handle.hpp:103
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2057
Data type sparse matrix.
Definition: sparsemat.hpp:41
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void EliminateRowsCols(OperatorHandle &A, const Array< int > &ess_dof_list)
Reset the OperatorHandle to be the eliminated part of A after elimination of the essential dofs ess_d...
Definition: handle.cpp:252
Operator::Type CheckType(Operator::Type tid)
Definition: handle.cpp:32
The operator x -&gt; R*A*P*x constructed through the actions of R^T, A and P.
Definition: operator.hpp:756
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:82
void MakeRectangularBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_num_rows, HYPRE_BigInt glob_num_cols, HYPRE_BigInt *row_starts, HYPRE_BigInt *col_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel rectangular block-diagonal matrix using the currently set...
Definition: handle.cpp:91
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:252
void MakeSquareBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_size, HYPRE_BigInt *row_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel square block-diagonal matrix using the currently set type...
Definition: handle.cpp:60
bool OwnsOperator() const
Return true if the OperatorHandle owns the held Operator.
Definition: handle.hpp:112
void EliminateRowCol(int rc, const double sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1708
HYPRE_Int HYPRE_BigInt
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:414
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:119
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:115
void pSet(OpType *A, bool own_A=true)
Definition: handle.hpp:45
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate &quot;essential boundary condition&quot; values specified in x from the given right-hand side b...
Definition: operator.cpp:452
Operator::Type type_id
Definition: handle.hpp:39
Keep the diagonal value.
Definition: operator.hpp:50
ID for class PetscParMatrix, MATAIJ format.
Definition: operator.hpp:257
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:46
Vector data type.
Definition: vector.hpp:60
ID for class HypreParMatrix.
Definition: operator.hpp:256
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:834
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
ID for class PetscParMatrix, MATIS format.
Definition: operator.hpp:258
void MakePtAP(OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product P^t A P.
Definition: handle.cpp:123
void MakeRAP(OperatorHandle &Rt, OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product R A P, where R = Rt^t.
Definition: handle.cpp:161
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:140