MFEM  v3.3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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::MFEM_SPARSEMAT: break;
38 #ifdef MFEM_USE_MPI
39  break;
40 #else
41  MFEM_ABORT("cannot use HYPRE parallel matrix format: "
42  "MFEM is not built with HYPRE support");
43 #endif
46 #ifdef MFEM_USE_PETSC
47  break;
48 #else
49  MFEM_ABORT("cannot use PETSc matrix formats: "
50  "MFEM is not built with PETSc support");
51 #endif
52  default:
53  MFEM_ABORT("invalid Operator::Type, type_id = " << (int)type_id);
54  }
55  return tid;
56 }
57 
58 #ifdef MFEM_USE_MPI
59 void OperatorHandle::MakeSquareBlockDiag(MPI_Comm comm, HYPRE_Int glob_size,
60  HYPRE_Int *row_starts,
61  SparseMatrix *diag)
62 {
63  if (own_oper) { delete oper; }
64 
65  switch (type_id)
66  {
68  oper = new HypreParMatrix(comm, glob_size, row_starts, diag);
69  break;
70 #ifdef MFEM_USE_PETSC
73  // Assuming that PetscInt is the same size as HYPRE_Int, checked above.
74  oper = new PetscParMatrix(comm, glob_size, row_starts, diag, type_id);
75  break;
76 #endif
77  default: MFEM_ABORT(not_supported_msg << type_id);
78  }
79  own_oper = true;
80 }
81 
83 MakeRectangularBlockDiag(MPI_Comm comm, HYPRE_Int glob_num_rows,
84  HYPRE_Int glob_num_cols, HYPRE_Int *row_starts,
85  HYPRE_Int *col_starts, SparseMatrix *diag)
86 {
87  if (own_oper) { delete oper; }
88 
89  switch (type_id)
90  {
92  oper = new HypreParMatrix(comm, glob_num_rows, glob_num_cols,
93  row_starts, col_starts, diag);
94  break;
95 #ifdef MFEM_USE_PETSC
98  // Assuming that PetscInt is the same size as HYPRE_Int, checked above.
99  oper = new PetscParMatrix(comm, glob_num_rows, glob_num_cols,
100  row_starts, col_starts, diag, type_id);
101  break;
102 #endif
103  default: MFEM_ABORT(not_supported_msg << type_id);
104  }
105  own_oper = true;
106 }
107 #endif // MFEM_USE_MPI
108 
110 {
111  MFEM_VERIFY(A.Type() == P.Type(), "type mismatch in A and P");
112  Clear();
113  switch (A.Type())
114  {
116  {
118  SparseMatrix *RA = mfem::Mult(*R, *A.As<SparseMatrix>());
119  delete R;
120  pSet(mfem::Mult(*RA, *P.As<SparseMatrix>()));
121  delete RA;
122  break;
123  }
124 #ifdef MFEM_USE_MPI
127  break;
128 #ifdef MFEM_USE_PETSC
131  {
133  break;
134  }
135 #endif
136 #endif
137  default: MFEM_ABORT(not_supported_msg << type_id);
138  }
139 }
140 
142  OperatorHandle &P)
143 {
144  MFEM_VERIFY(A.Type() == Rt.Type(), "type mismatch in A and Rt");
145  MFEM_VERIFY(A.Type() == P.Type(), "type mismatch in A and P");
146  Clear();
147  switch (A.Type())
148  {
150  {
152  *P.As<SparseMatrix>()));
153  break;
154  }
155 #ifdef MFEM_USE_MPI
158  P.As<HypreParMatrix>()));
159  break;
160 #ifdef MFEM_USE_PETSC
163  {
165  P.As<PetscParMatrix>()));
166  break;
167  }
168 #endif
169 #endif
170  default: MFEM_ABORT(not_supported_msg << type_id);
171  }
172 }
173 
175 {
176  if (own_oper) { delete oper; }
177  if (Type() == A.Type())
178  {
179  oper = A.Ptr();
180  own_oper = false;
181  return;
182  }
183  oper = NULL;
184  switch (Type()) // target type id
185  {
188  switch (A.Type()) // source type id
189  {
191 #ifdef MFEM_USE_PETSC
192  oper = new PetscParMatrix(A.As<HypreParMatrix>(), Type());
193 #endif
194  break;
195  default: break;
196  }
197  break;
198  default: break;
199  }
200  MFEM_VERIFY(oper != NULL, "conversion from type id = " << A.Type()
201  << " to type id = " << Type() << " is not supported");
202  own_oper = true;
203 }
204 
206  const Array<int> &ess_dof_list)
207 {
208  Clear();
209  switch (A.Type())
210  {
212  {
213  const int preserve_diag = 1;
214  SparseMatrix *sA = A.As<SparseMatrix>();
215  SparseMatrix *Ae = new SparseMatrix(sA->Height());
216  for (int i = 0; i < ess_dof_list.Size(); i++)
217  {
218  sA->EliminateRowCol(ess_dof_list[i], *Ae, preserve_diag);
219  }
220  Ae->Finalize();
221  pSet(Ae);
222  break;
223  }
225  {
226 #ifdef MFEM_USE_MPI
227  pSet(A.As<HypreParMatrix>()->EliminateRowsCols(ess_dof_list));
228 #else
229  MFEM_ABORT("type id = Hypre_ParCSR requires MFEM_USE_MPI");
230 #endif
231  break;
232  }
235  {
236 #ifdef MFEM_USE_PETSC
237  pSet(A.As<PetscParMatrix>()->EliminateRowsCols(ess_dof_list));
238 #else
239  MFEM_ABORT("type id = Operator::PETSC_* requires MFEM_USE_PETSC");
240 #endif
241  break;
242  }
243  default: MFEM_ABORT(not_supported_msg << A.Type());
244  }
245 }
246 
248  const Array<int> &ess_dof_list,
249  const Vector &X, Vector &B) const
250 {
251  switch (Type())
252  {
254  {
255  A_e.As<SparseMatrix>()->AddMult(X, B, -1.);
256  As<SparseMatrix>()->PartMult(ess_dof_list, X, B);
257  break;
258  }
260  {
261 #ifdef MFEM_USE_MPI
262  mfem::EliminateBC(*As<HypreParMatrix>(), *A_e.As<HypreParMatrix>(),
263  ess_dof_list, X, B);
264 #else
265  MFEM_ABORT("type id = Hypre_ParCSR requires MFEM_USE_MPI");
266 #endif
267  break;
268  }
271  {
272 #ifdef MFEM_USE_PETSC
273  mfem::EliminateBC(*As<PetscParMatrix>(), *A_e.As<PetscParMatrix>(),
274  ess_dof_list, X, B);
275 #else
276  MFEM_ABORT("type id = Operator::PETSC_* requires MFEM_USE_PETSC");
277 #endif
278  break;
279  }
280  default: MFEM_ABORT(not_supported_msg << Type());
281  }
282 }
283 
284 } // namespace mfem
int Size() const
Logical size of the array.
Definition: array.hpp:110
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:1286
void ConvertFrom(OperatorHandle &A)
Convert the given OperatorHandle A to the currently set type id.
Definition: handle.cpp:174
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:90
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1476
Wrapper for PETSc&#39;s matrix class.
Definition: petsc.hpp:133
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
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:85
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:247
void EliminateRowCol(int rc, const double sol, Vector &rhs, int d=0)
Definition: sparsemat.cpp:1162
Operator * oper
Definition: handle.hpp:38
static const char not_supported_msg[]
Definition: handle.hpp:36
void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Definition: hypre.cpp:1530
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:3016
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
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:205
Operator::Type CheckType(Operator::Type tid)
Definition: handle.cpp:32
void EliminateRowsCols(const Array< int > &rows_cols, const PetscParVector &X, PetscParVector &B)
Eliminate rows and columns from the matrix, and rows from the vector B. Modify B with the BC values i...
Definition: petsc.cpp:1219
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:408
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:82
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:123
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:110
void MakeSquareBlockDiag(MPI_Comm comm, HYPRE_Int glob_size, HYPRE_Int *row_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel square block-diagonal matrix using the currently set type...
Definition: handle.cpp:59
void pSet(OpType *A, bool own_A=true)
Definition: handle.hpp:45
Operator::Type type_id
Definition: handle.hpp:39
void MakeRectangularBlockDiag(MPI_Comm comm, HYPRE_Int glob_num_rows, HYPRE_Int glob_num_cols, HYPRE_Int *row_starts, HYPRE_Int *col_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel rectangular block-diagonal matrix using the currently set...
Definition: handle.cpp:83
Vector data type.
Definition: vector.hpp:41
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
void MakePtAP(OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product P^t A P.
Definition: handle.cpp:109
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:141