MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
petsc.cpp
Go to the documentation of this file.
1 // Copyright (c) 2016, 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 // Author: Stefano Zampini <stefano.zampini@gmail.com>
13 
14 #include "../config/config.hpp"
15 
16 #ifdef MFEM_USE_MPI
17 #ifdef MFEM_USE_PETSC
18 
19 #include "linalg.hpp"
20 #include "../fem/fem.hpp"
21 #if defined(PETSC_HAVE_HYPRE)
22 #include "petscmathypre.h"
23 #endif
24 
25 #include <fstream>
26 #include <iomanip>
27 #include <cmath>
28 #include <cstdlib>
29 // Note: there are additional #include statements below.
30 
31 // Error handling
32 // Prints PETSc's stacktrace and then calls MFEM_ABORT
33 // We cannot use PETSc's CHKERRQ since it returns a PetscErrorCode
34 #define PCHKERRQ(obj,err) do { \
35  if ((err)) \
36  { \
37  PetscError(PetscObjectComm((PetscObject)(obj)),__LINE__,_MFEM_FUNC_NAME, \
38  __FILE__,(err),PETSC_ERROR_REPEAT,NULL); \
39  MFEM_ABORT("Error in PETSc. See stacktrace above."); \
40  } \
41  } while(0);
42 #define CCHKERRQ(comm,err) do { \
43  if ((err)) \
44  { \
45  PetscError(comm,__LINE__,_MFEM_FUNC_NAME, \
46  __FILE__,(err),PETSC_ERROR_REPEAT,NULL); \
47  MFEM_ABORT("Error in PETSc. See stacktrace above."); \
48  } \
49  } while(0);
50 
51 // Callback functions: these functions will be called by PETSc
52 static PetscErrorCode __mfem_ts_monitor(TS,PetscInt,PetscReal,Vec,void*);
53 static PetscErrorCode __mfem_ts_rhsfunction(TS,PetscReal,Vec,Vec,void*);
54 static PetscErrorCode __mfem_ts_rhsjacobian(TS,PetscReal,Vec,Mat,Mat,
55  void*);
56 static PetscErrorCode __mfem_ts_ifunction(TS,PetscReal,Vec,Vec,Vec,void*);
57 static PetscErrorCode __mfem_ts_ijacobian(TS,PetscReal,Vec,Vec,
58  PetscReal,Mat,
59  Mat,void*);
60 static PetscErrorCode __mfem_snes_monitor(SNES,PetscInt,PetscReal,void*);
61 static PetscErrorCode __mfem_snes_jacobian(SNES,Vec,Mat,Mat,void*);
62 static PetscErrorCode __mfem_snes_function(SNES,Vec,Vec,void*);
63 static PetscErrorCode __mfem_ksp_monitor(KSP,PetscInt,PetscReal,void*);
64 static PetscErrorCode __mfem_pc_shell_apply(PC,Vec,Vec);
65 static PetscErrorCode __mfem_pc_shell_apply_transpose(PC,Vec,Vec);
66 static PetscErrorCode __mfem_pc_shell_setup(PC);
67 static PetscErrorCode __mfem_pc_shell_destroy(PC);
68 static PetscErrorCode __mfem_pc_shell_view(PC,PetscViewer);
69 static PetscErrorCode __mfem_mat_shell_apply(Mat,Vec,Vec);
70 static PetscErrorCode __mfem_mat_shell_apply_transpose(Mat,Vec,Vec);
71 static PetscErrorCode __mfem_mat_shell_destroy(Mat);
72 static PetscErrorCode __mfem_array_container_destroy(void*);
73 static PetscErrorCode __mfem_matarray_container_destroy(void*);
74 
75 // auxiliary functions
76 static PetscErrorCode Convert_Array_IS(MPI_Comm,bool,const mfem::Array<int>*,
77  PetscInt,IS*);
78 static PetscErrorCode Convert_Vmarks_IS(MPI_Comm,mfem::Array<Mat>&,
79  const mfem::Array<int>*,PetscInt,IS*);
80 static PetscErrorCode MakeShellPC(PC,mfem::Solver&,bool);
81 static PetscErrorCode MakeShellPCWithFactory(PC,
82  mfem::PetscPreconditionerFactory*);
83 
84 // Equivalent functions are present in PETSc source code
85 // if PETSc has been compiled with hypre support
86 // We provide them here in case PETSC_HAVE_HYPRE is not defined
87 #if !defined(PETSC_HAVE_HYPRE)
88 static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix*,Mat*);
89 static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix*,Mat*);
90 #endif
91 
92 // structs used by PETSc code
93 typedef struct
94 {
95  mfem::Operator *op;
96 } __mfem_mat_shell_ctx;
97 
98 typedef struct
99 {
100  mfem::Solver *op;
101  mfem::PetscPreconditionerFactory *factory;
102  bool ownsop;
103  unsigned long int numprec;
104 } __mfem_pc_shell_ctx;
105 
106 typedef struct
107 {
108  mfem::Operator *op; // The nonlinear operator
109  mfem::PetscBCHandler *bchandler; // Handling of essential bc
110  mfem::Vector *work; // Work vector
111  mfem::Operator::Type jacType; // OperatorType for the Jacobian
112 } __mfem_snes_ctx;
113 
114 typedef struct
115 {
116  mfem::TimeDependentOperator *op; // The time-dependent operator
117  mfem::PetscBCHandler *bchandler; // Handling of essential bc
118  mfem::Vector *work; // Work vector
119  mfem::Operator::Type jacType; // OperatorType for the Jacobian
120  enum mfem::PetscODESolver::Type type;
121  PetscReal cached_shift;
122  PetscObjectState cached_ijacstate;
123  PetscObjectState cached_rhsjacstate;
124 } __mfem_ts_ctx;
125 
126 // use global scope ierr to check PETSc errors inside mfem calls
127 static PetscErrorCode ierr;
128 
129 using namespace std;
130 
131 namespace mfem
132 {
133 
134 // PetscParVector methods
135 
136 void PetscParVector::_SetDataAndSize_()
137 {
138  const PetscScalar *array;
139  PetscInt n;
140 
141  ierr = VecGetArrayRead(x,&array); PCHKERRQ(x,ierr);
142  ierr = VecGetLocalSize(x,&n); PCHKERRQ(x,ierr);
143  SetDataAndSize((PetscScalar*)array,n);
144  ierr = VecRestoreArrayRead(x,&array); PCHKERRQ(x,ierr);
145 }
146 
147 PetscInt PetscParVector::GlobalSize() const
148 {
149  PetscInt N;
150  ierr = VecGetSize(x,&N); PCHKERRQ(x,ierr);
151  return N;
152 }
153 
154 PetscParVector::PetscParVector(MPI_Comm comm, const Vector &_x) : Vector()
155 {
156  ierr = VecCreate(comm,&x); CCHKERRQ(comm,ierr);
157  ierr = VecSetSizes(x,_x.Size(),PETSC_DECIDE); PCHKERRQ(x,ierr);
158  ierr = VecSetType(x,VECSTANDARD); PCHKERRQ(x,ierr);
159  _SetDataAndSize_();
160 }
161 
162 PetscParVector::PetscParVector(MPI_Comm comm, PetscInt glob_size,
163  PetscInt *col) : Vector()
164 {
165  ierr = VecCreate(comm,&x); CCHKERRQ(comm,ierr);
166  if (col)
167  {
168  PetscMPIInt myid;
169  MPI_Comm_rank(comm, &myid);
170  ierr = VecSetSizes(x,col[myid+1]-col[myid],PETSC_DECIDE); PCHKERRQ(x,ierr);
171  }
172  else
173  {
174  ierr = VecSetSizes(x,PETSC_DECIDE,glob_size); PCHKERRQ(x,ierr);
175  }
176  ierr = VecSetType(x,VECSTANDARD); PCHKERRQ(x,ierr);
177  _SetDataAndSize_();
178 }
179 
180 PetscParVector::~PetscParVector()
181 {
182  MPI_Comm comm = PetscObjectComm((PetscObject)x);
183  ierr = VecDestroy(&x); CCHKERRQ(comm,ierr);
184 }
185 
186 PetscParVector::PetscParVector(MPI_Comm comm, PetscInt glob_size,
187  PetscScalar *_data, PetscInt *col) : Vector()
188 {
189  MFEM_VERIFY(col,"Missing distribution");
190  PetscMPIInt myid;
191  MPI_Comm_rank(comm, &myid);
192  ierr = VecCreateMPIWithArray(comm,1,col[myid+1]-col[myid],glob_size,_data,
193  &x); CCHKERRQ(comm,ierr)
194  _SetDataAndSize_();
195 }
196 
197 PetscParVector::PetscParVector(const PetscParVector &y) : Vector()
198 {
199  ierr = VecDuplicate(y.x,&x); PCHKERRQ(x,ierr);
200  _SetDataAndSize_();
201 }
202 
203 PetscParVector::PetscParVector(MPI_Comm comm, const Operator &op,
204  bool transpose, bool allocate) : Vector()
205 {
206  PetscInt loc = transpose ? op.Height() : op.Width();
207  if (allocate)
208  {
209  ierr = VecCreate(comm,&x);
210  CCHKERRQ(comm,ierr);
211  ierr = VecSetSizes(x,loc,PETSC_DECIDE);
212  PCHKERRQ(x,ierr);
213  ierr = VecSetType(x,VECSTANDARD);
214  PCHKERRQ(x,ierr);
215  ierr = VecSetUp(x);
216  PCHKERRQ(x,ierr);
217  }
218  else
219  {
220  ierr = VecCreateMPIWithArray(comm,1,loc,PETSC_DECIDE,NULL,
221  &x); CCHKERRQ(comm,ierr);
222  }
223  _SetDataAndSize_();
224 }
225 
226 PetscParVector::PetscParVector(const PetscParMatrix &A,
227  bool transpose, bool allocate) : Vector()
228 {
229  Mat pA = const_cast<PetscParMatrix&>(A);
230  if (!transpose)
231  {
232  ierr = MatCreateVecs(pA,&x,NULL); PCHKERRQ(pA,ierr);
233  }
234  else
235  {
236  ierr = MatCreateVecs(pA,NULL,&x); PCHKERRQ(pA,ierr);
237  }
238  if (!allocate)
239  {
240  ierr = VecReplaceArray(x,NULL); PCHKERRQ(x,ierr);
241  }
242  _SetDataAndSize_();
243 }
244 
245 PetscParVector::PetscParVector(Vec y, bool ref) : Vector()
246 {
247  if (ref)
248  {
249  ierr = PetscObjectReference((PetscObject)y); PCHKERRQ(y,ierr);
250  }
251  x = y;
252  _SetDataAndSize_();
253 }
254 
255 PetscParVector::PetscParVector(ParFiniteElementSpace *pfes) : Vector()
256 {
257 
258  HYPRE_Int* offsets = pfes->GetTrueDofOffsets();
259  MPI_Comm comm = pfes->GetComm();
260  ierr = VecCreate(comm,&x); CCHKERRQ(comm,ierr);
261 
262  PetscMPIInt myid = 0;
263  if (!HYPRE_AssumedPartitionCheck())
264  {
265  MPI_Comm_rank(comm,&myid);
266  }
267  ierr = VecSetSizes(x,offsets[myid+1]-offsets[myid],PETSC_DECIDE);
268  PCHKERRQ(x,ierr);
269  ierr = VecSetType(x,VECSTANDARD); PCHKERRQ(x,ierr);
270  _SetDataAndSize_();
271 }
272 
273 Vector * PetscParVector::GlobalVector() const
274 {
275  VecScatter scctx;
276  Vec vout;
277  PetscScalar *array;
278  PetscInt size;
279 
280  ierr = VecScatterCreateToAll(x,&scctx,&vout); PCHKERRQ(x,ierr);
281  ierr = VecScatterBegin(scctx,x,vout,INSERT_VALUES,SCATTER_FORWARD);
282  PCHKERRQ(x,ierr);
283  ierr = VecScatterEnd(scctx,x,vout,INSERT_VALUES,SCATTER_FORWARD);
284  PCHKERRQ(x,ierr);
285  ierr = VecScatterDestroy(&scctx); PCHKERRQ(x,ierr);
286  ierr = VecGetArray(vout,&array); PCHKERRQ(x,ierr);
287  ierr = VecGetLocalSize(vout,&size); PCHKERRQ(x,ierr);
288  Array<PetscScalar> data(size);
289  data.Assign(array);
290  ierr = VecRestoreArray(vout,&array); PCHKERRQ(x,ierr);
291  ierr = VecDestroy(&vout); PCHKERRQ(x,ierr);
292  Vector *v = new Vector(data, internal::to_int(size));
293  v->MakeDataOwner();
294  data.LoseData();
295  return v;
296 }
297 
298 PetscParVector& PetscParVector::operator=(PetscScalar d)
299 {
300  ierr = VecSet(x,d); PCHKERRQ(x,ierr);
301  return *this;
302 }
303 
304 PetscParVector& PetscParVector::operator=(const PetscParVector &y)
305 {
306  ierr = VecCopy(y.x,x); PCHKERRQ(x,ierr);
307  return *this;
308 }
309 
310 void PetscParVector::PlaceArray(PetscScalar *temp_data)
311 {
312  ierr = VecPlaceArray(x,temp_data); PCHKERRQ(x,ierr);
313 }
314 
315 void PetscParVector::ResetArray()
316 {
317  ierr = VecResetArray(x); PCHKERRQ(x,ierr);
318 }
319 
320 void PetscParVector::Randomize(PetscInt seed)
321 {
322  PetscRandom rctx;
323 
324  ierr = PetscRandomCreate(PetscObjectComm((PetscObject)x),&rctx);
325  PCHKERRQ(x,ierr);
326  ierr = PetscRandomSetSeed(rctx,(unsigned long)seed); PCHKERRQ(x,ierr);
327  ierr = PetscRandomSeed(rctx); PCHKERRQ(x,ierr);
328  ierr = VecSetRandom(x,rctx); PCHKERRQ(x,ierr);
329  ierr = PetscRandomDestroy(&rctx); PCHKERRQ(x,ierr);
330 }
331 
332 void PetscParVector::Print(const char *fname, bool binary) const
333 {
334  if (fname)
335  {
336  PetscViewer view;
337 
338  if (binary)
339  {
340  ierr = PetscViewerBinaryOpen(PetscObjectComm((PetscObject)x),fname,
341  FILE_MODE_WRITE,&view);
342  }
343  else
344  {
345  ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)x),fname,&view);
346  }
347  PCHKERRQ(x,ierr);
348  ierr = VecView(x,view); PCHKERRQ(x,ierr);
349  ierr = PetscViewerDestroy(&view); PCHKERRQ(x,ierr);
350  }
351  else
352  {
353  ierr = VecView(x,NULL); PCHKERRQ(x,ierr);
354  }
355 }
356 
357 // PetscParMatrix methods
358 
359 PetscInt PetscParMatrix::GetNumRows() const
360 {
361  PetscInt N;
362  ierr = MatGetLocalSize(A,&N,NULL); PCHKERRQ(A,ierr);
363  return N;
364 }
365 
366 PetscInt PetscParMatrix::GetNumCols() const
367 {
368  PetscInt N;
369  ierr = MatGetLocalSize(A,NULL,&N); PCHKERRQ(A,ierr);
370  return N;
371 }
372 
373 PetscInt PetscParMatrix::M() const
374 {
375  PetscInt N;
376  ierr = MatGetSize(A,&N,NULL); PCHKERRQ(A,ierr);
377  return N;
378 }
379 
380 PetscInt PetscParMatrix::N() const
381 {
382  PetscInt N;
383  ierr = MatGetSize(A,NULL,&N); PCHKERRQ(A,ierr);
384  return N;
385 }
386 
387 PetscInt PetscParMatrix::NNZ() const
388 {
389  MatInfo info;
390  ierr = MatGetInfo(A,MAT_GLOBAL_SUM,&info); PCHKERRQ(A,ierr);
391  return (PetscInt)info.nz_used;
392 }
393 
394 void PetscParMatrix::Init()
395 {
396  A = NULL;
397  X = Y = NULL;
398  height = width = 0;
399 }
400 
401 PetscParMatrix::PetscParMatrix()
402 {
403  Init();
404 }
405 
406 PetscParMatrix::PetscParMatrix(const HypreParMatrix *ha, Operator::Type tid)
407 {
408  Init();
409  height = ha->Height();
410  width = ha->Width();
411  ConvertOperator(ha->GetComm(),*ha,&A,tid);
412 }
413 
414 PetscParMatrix::PetscParMatrix(MPI_Comm comm, const Operator *op,
415  Operator::Type tid)
416 {
417  Init();
418  height = op->Height();
419  width = op->Width();
420  ConvertOperator(comm,*op,&A,tid);
421 }
422 
423 PetscParMatrix::PetscParMatrix(MPI_Comm comm, PetscInt glob_size,
424  PetscInt *row_starts, SparseMatrix *diag,
425  Operator::Type tid)
426 {
427  Init();
428  BlockDiagonalConstructor(comm,row_starts,row_starts,diag,
429  tid==PETSC_MATAIJ,&A);
430  // update base class
431  height = GetNumRows();
432  width = GetNumCols();
433 }
434 
435 PetscParMatrix::PetscParMatrix(MPI_Comm comm, PetscInt global_num_rows,
436  PetscInt global_num_cols, PetscInt *row_starts,
437  PetscInt *col_starts, SparseMatrix *diag,
438  Operator::Type tid)
439 {
440  Init();
441  BlockDiagonalConstructor(comm,row_starts,col_starts,diag,
442  tid==PETSC_MATAIJ,&A);
443  // update base class
444  height = GetNumRows();
445  width = GetNumCols();
446 }
447 
448 PetscParMatrix& PetscParMatrix::operator=(const HypreParMatrix& B)
449 {
450  if (A)
451  {
452  MPI_Comm comm = PetscObjectComm((PetscObject)A);
453  ierr = MatDestroy(&A); CCHKERRQ(comm,ierr);
454  if (X) { delete X; }
455  if (Y) { delete Y; }
456  X = Y = NULL;
457  }
458  height = B.Height();
459  width = B.Width();
460 #if defined(PETSC_HAVE_HYPRE)
461  ierr = MatCreateFromParCSR(B,MATAIJ,PETSC_USE_POINTER,&A);
462 #else
463  ierr = MatConvert_hypreParCSR_AIJ(B,&A); CCHKERRQ(B.GetComm(),ierr);
464 #endif
465  return *this;
466 }
467 
468 PetscParMatrix& PetscParMatrix::operator=(const PetscParMatrix& B)
469 {
470  if (A)
471  {
472  MPI_Comm comm = PetscObjectComm((PetscObject)A);
473  ierr = MatDestroy(&A); CCHKERRQ(comm,ierr);
474  if (X) { delete X; }
475  if (Y) { delete Y; }
476  X = Y = NULL;
477  }
478  height = B.Height();
479  width = B.Width();
480  ierr = MatDuplicate(B,MAT_COPY_VALUES,&A); CCHKERRQ(B.GetComm(),ierr);
481  return *this;
482 }
483 
484 PetscParMatrix& PetscParMatrix::operator+=(const PetscParMatrix& B)
485 {
486  if (!A)
487  {
488  ierr = MatDuplicate(B,MAT_COPY_VALUES,&A); CCHKERRQ(B.GetComm(),ierr);
489  }
490  else
491  {
492  MFEM_VERIFY(height == B.Height(),"Invalid number of local rows");
493  MFEM_VERIFY(width == B.Width(), "Invalid number of local columns");
494  ierr = MatAXPY(A,1.0,B,DIFFERENT_NONZERO_PATTERN); CCHKERRQ(B.GetComm(),ierr);
495  }
496  return *this;
497 }
498 
499 void PetscParMatrix::
500 BlockDiagonalConstructor(MPI_Comm comm,
501  PetscInt *row_starts, PetscInt *col_starts,
502  SparseMatrix *diag, bool assembled, Mat* Ad)
503 {
504  Mat A;
505  PetscInt lrsize,lcsize,rstart,cstart;
506  PetscMPIInt myid = 0,commsize;
507 
508  ierr = MPI_Comm_size(comm,&commsize); CCHKERRQ(comm,ierr);
509  if (!HYPRE_AssumedPartitionCheck())
510  {
511  ierr = MPI_Comm_rank(comm,&myid); CCHKERRQ(comm,ierr);
512  }
513  lrsize = row_starts[myid+1]-row_starts[myid];
514  rstart = row_starts[myid];
515  lcsize = col_starts[myid+1]-col_starts[myid];
516  cstart = col_starts[myid];
517 
518  if (!assembled)
519  {
520  IS is;
521  ierr = ISCreateStride(comm,diag->Height(),rstart,1,&is); CCHKERRQ(comm,ierr);
522  ISLocalToGlobalMapping rl2g,cl2g;
523  ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g); PCHKERRQ(is,ierr);
524  ierr = ISDestroy(&is); CCHKERRQ(comm,ierr);
525  if (row_starts != col_starts)
526  {
527  ierr = ISCreateStride(comm,diag->Width(),cstart,1,&is);
528  CCHKERRQ(comm,ierr);
529  ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g); PCHKERRQ(is,ierr);
530  ierr = ISDestroy(&is); CCHKERRQ(comm,ierr);
531  }
532  else
533  {
534  ierr = PetscObjectReference((PetscObject)rl2g); PCHKERRQ(rl2g,ierr);
535  cl2g = rl2g;
536  }
537 
538  // Create the PETSc object (MATIS format)
539  ierr = MatCreate(comm,&A); CCHKERRQ(comm,ierr);
540  ierr = MatSetSizes(A,lrsize,lcsize,PETSC_DECIDE,PETSC_DECIDE);
541  PCHKERRQ(A,ierr);
542  ierr = MatSetType(A,MATIS); PCHKERRQ(A,ierr);
543  ierr = MatSetLocalToGlobalMapping(A,rl2g,cl2g); PCHKERRQ(A,ierr);
544  ierr = ISLocalToGlobalMappingDestroy(&rl2g); PCHKERRQ(A,ierr)
545  ierr = ISLocalToGlobalMappingDestroy(&cl2g); PCHKERRQ(A,ierr)
546 
547  // Copy SparseMatrix into PETSc SeqAIJ format
548  Mat lA;
549  ierr = MatISGetLocalMat(A,&lA); PCHKERRQ(A,ierr);
550  if (sizeof(PetscInt) == sizeof(int))
551  {
552  ierr = MatSeqAIJSetPreallocationCSR(lA,diag->GetI(),diag->GetJ(),
553  diag->GetData()); PCHKERRQ(lA,ierr);
554  }
555  else
556  {
557  MFEM_ABORT("64bit indices not yet supported");
558  }
559  }
560  else
561  {
562  PetscScalar *da;
563  PetscInt *dii,*djj,*oii,
564  m = diag->Height()+1, nnz = diag->NumNonZeroElems();
565 
566  diag->SortColumnIndices();
567  // if we can take ownership of the SparseMatrix arrays, we can avoid this
568  // step
569  ierr = PetscMalloc1(m,&dii); CCHKERRQ(PETSC_COMM_SELF,ierr);
570  ierr = PetscMalloc1(nnz,&djj); CCHKERRQ(PETSC_COMM_SELF,ierr);
571  ierr = PetscMalloc1(nnz,&da); CCHKERRQ(PETSC_COMM_SELF,ierr);
572  if (sizeof(PetscInt) == sizeof(int))
573  {
574  ierr = PetscMemcpy(dii,diag->GetI(),m*sizeof(PetscInt));
575  CCHKERRQ(PETSC_COMM_SELF,ierr);
576  ierr = PetscMemcpy(djj,diag->GetJ(),nnz*sizeof(PetscInt));
577  CCHKERRQ(PETSC_COMM_SELF,ierr);
578  ierr = PetscMemcpy(da,diag->GetData(),nnz*sizeof(PetscScalar));
579  CCHKERRQ(PETSC_COMM_SELF,ierr);
580  }
581  else
582  {
583  MFEM_ABORT("64bit indices not yet supported");
584  }
585  ierr = PetscCalloc1(m,&oii);
586  CCHKERRQ(PETSC_COMM_SELF,ierr);
587  if (commsize > 1)
588  {
589  ierr = MatCreateMPIAIJWithSplitArrays(comm,lrsize,lcsize,PETSC_DECIDE,
590  PETSC_DECIDE,
591  dii,djj,da,oii,NULL,NULL,&A);
592  CCHKERRQ(comm,ierr);
593  }
594  else
595  {
596  ierr = MatCreateSeqAIJWithArrays(comm,lrsize,lcsize,dii,djj,da,&A);
597  CCHKERRQ(comm,ierr);
598  }
599 
600  void *ptrs[4] = {dii,djj,da,oii};
601  const char *names[4] = {"_mfem_csr_dii",
602  "_mfem_csr_djj",
603  "_mfem_csr_da",
604  "_mfem_csr_oii",
605  };
606  for (PetscInt i=0; i<4; i++)
607  {
608  PetscContainer c;
609 
610  ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
611  ierr = PetscContainerSetPointer(c,ptrs[i]); CCHKERRQ(comm,ierr);
612  ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
613  CCHKERRQ(comm,ierr);
614  ierr = PetscObjectCompose((PetscObject)A,names[i],(PetscObject)c);
615  CCHKERRQ(comm,ierr);
616  ierr = PetscContainerDestroy(&c); CCHKERRQ(comm,ierr);
617  }
618  }
619 
620  // Tell PETSc the matrix is ready to be used
621  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); PCHKERRQ(A,ierr);
622  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); PCHKERRQ(A,ierr);
623 
624  *Ad = A;
625 }
626 
627 // TODO ADD THIS CONSTRUCTOR
628 //PetscParMatrix::PetscParMatrix(MPI_Comm comm, int nrows, PetscInt glob_nrows,
629 // PetscInt glob_ncols, int *I, PetscInt *J,
630 // double *data, PetscInt *rows, PetscInt *cols)
631 //{
632 //}
633 
634 // TODO This should take a reference on op but how?
635 void PetscParMatrix::MakeWrapper(MPI_Comm comm, const Operator* op, Mat *A)
636 {
637  __mfem_mat_shell_ctx *ctx = new __mfem_mat_shell_ctx;
638  ierr = MatCreate(comm,A); CCHKERRQ(comm,ierr);
639  ierr = MatSetSizes(*A,op->Height(),op->Width(),
640  PETSC_DECIDE,PETSC_DECIDE); PCHKERRQ(A,ierr);
641  ierr = MatSetType(*A,MATSHELL); PCHKERRQ(A,ierr);
642  ierr = MatShellSetContext(*A,(void *)ctx); PCHKERRQ(A,ierr);
643  ierr = MatShellSetOperation(*A,MATOP_MULT,
644  (void (*)())__mfem_mat_shell_apply);
645  PCHKERRQ(A,ierr);
646  ierr = MatShellSetOperation(*A,MATOP_MULT_TRANSPOSE,
647  (void (*)())__mfem_mat_shell_apply_transpose);
648  PCHKERRQ(A,ierr);
649  ierr = MatShellSetOperation(*A,MATOP_DESTROY,
650  (void (*)())__mfem_mat_shell_destroy);
651  PCHKERRQ(A,ierr);
652  ierr = MatSetUp(*A); PCHKERRQ(*A,ierr);
653  ctx->op = const_cast<Operator *>(op);
654 }
655 
656 void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
657  Operator::Type tid)
658 {
659  PetscParMatrix *pA = const_cast<PetscParMatrix *>
660  (dynamic_cast<const PetscParMatrix *>(&op));
661  HypreParMatrix *pH = const_cast<HypreParMatrix *>
662  (dynamic_cast<const HypreParMatrix *>(&op));
663  BlockOperator *pB = const_cast<BlockOperator *>
664  (dynamic_cast<const BlockOperator *>(&op));
665  IdentityOperator *pI = const_cast<IdentityOperator *>
666  (dynamic_cast<const IdentityOperator *>(&op));
667 
668  if (pA)
669  {
670  Mat At = NULL;
671  PetscBool ismatis,istrans;
672 
673  ierr = PetscObjectTypeCompare((PetscObject)(pA->A),MATTRANSPOSEMAT,&istrans);
674  CCHKERRQ(pA->GetComm(),ierr);
675  if (!istrans)
676  {
677  if (tid == pA->GetType()) // use same object and return
678  {
679  ierr = PetscObjectReference((PetscObject)(pA->A));
680  CCHKERRQ(pA->GetComm(),ierr);
681  *A = pA->A;
682  return;
683  }
684  ierr = PetscObjectTypeCompare((PetscObject)(pA->A),MATIS,&ismatis);
685  CCHKERRQ(pA->GetComm(),ierr);
686  }
687  else
688  {
689  ierr = MatTransposeGetMat(pA->A,&At); CCHKERRQ(pA->GetComm(),ierr);
690  ierr = PetscObjectTypeCompare((PetscObject)(At),MATIS,&ismatis);
691  CCHKERRQ(pA->GetComm(),ierr);
692  }
693 
694  // Try to convert
695  if (tid == PETSC_MATAIJ)
696  {
697  if (ismatis)
698  {
699  if (istrans)
700  {
701  Mat B;
702 
703  ierr = MatISGetMPIXAIJ(At,MAT_INITIAL_MATRIX,&B); PCHKERRQ(pA->A,ierr);
704  ierr = MatCreateTranspose(B,A); PCHKERRQ(pA->A,ierr);
705  ierr = MatDestroy(&B); PCHKERRQ(pA->A,ierr);
706  }
707  else
708  {
709  ierr = MatISGetMPIXAIJ(pA->A,MAT_INITIAL_MATRIX,A);
710  PCHKERRQ(pA->A,ierr);
711  }
712  }
713  else
714  {
715  PetscMPIInt size;
716  ierr = MPI_Comm_size(comm,&size); CCHKERRQ(comm,ierr);
717 
718  // call MatConvert and see if a converter is available
719  if (istrans)
720  {
721  Mat B;
722  ierr = MatConvert(At,size > 1 ? MATMPIAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
723  PCHKERRQ(pA->A,ierr);
724  ierr = MatCreateTranspose(B,A); PCHKERRQ(pA->A,ierr);
725  ierr = MatDestroy(&B); PCHKERRQ(pA->A,ierr);
726  }
727  else
728  {
729  ierr = MatConvert(pA->A, size > 1 ? MATMPIAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,A);
730  PCHKERRQ(pA->A,ierr);
731  }
732  }
733  }
734  else if (tid == PETSC_MATIS)
735  {
736  if (istrans)
737  {
738  Mat B;
739  ierr = MatConvert(At,MATIS,MAT_INITIAL_MATRIX,&B); PCHKERRQ(pA->A,ierr);
740  ierr = MatCreateTranspose(B,A); PCHKERRQ(pA->A,ierr);
741  ierr = MatDestroy(&B); PCHKERRQ(pA->A,ierr);
742  }
743  else
744  {
745  ierr = MatConvert(pA->A,MATIS,MAT_INITIAL_MATRIX,A); PCHKERRQ(pA->A,ierr);
746  }
747  }
748 #if defined(PETSC_HAVE_HYPRE)
749  else if (tid == PETSC_MATHYPRE)
750  {
751  if (istrans)
752  {
753  Mat B;
754  ierr = MatConvert(At,MATHYPRE,MAT_INITIAL_MATRIX,&B); PCHKERRQ(pA->A,ierr);
755  ierr = MatCreateTranspose(B,A); PCHKERRQ(pA->A,ierr);
756  ierr = MatDestroy(&B); PCHKERRQ(pA->A,ierr);
757  }
758  else
759  {
760  ierr = MatConvert(pA->A,MATHYPRE,MAT_INITIAL_MATRIX,A); PCHKERRQ(pA->A,ierr);
761  }
762  }
763 #endif
764  else if (tid == PETSC_MATSHELL)
765  {
766  MakeWrapper(comm,&op,A);
767  }
768  else
769  {
770  MFEM_ABORT("Unsupported operator type conversion " << tid)
771  }
772  }
773  else if (pH)
774  {
775  if (tid == PETSC_MATAIJ)
776  {
777 #if defined(PETSC_HAVE_HYPRE)
778  ierr = MatCreateFromParCSR(const_cast<HypreParMatrix&>(*pH),MATAIJ,
779  PETSC_USE_POINTER,A);
780 #else
781  ierr = MatConvert_hypreParCSR_AIJ(const_cast<HypreParMatrix&>(*pH),A);
782 #endif
783  CCHKERRQ(pH->GetComm(),ierr);
784  }
785  else if (tid == PETSC_MATIS)
786  {
787 #if defined(PETSC_HAVE_HYPRE)
788  ierr = MatCreateFromParCSR(const_cast<HypreParMatrix&>(*pH),MATIS,
789  PETSC_USE_POINTER,A);
790 #else
791  ierr = MatConvert_hypreParCSR_IS(const_cast<HypreParMatrix&>(*pH),A);
792 #endif
793  CCHKERRQ(pH->GetComm(),ierr);
794  }
795 #if defined(PETSC_HAVE_HYPRE)
796  else if (tid == PETSC_MATHYPRE)
797  {
798  ierr = MatCreateFromParCSR(const_cast<HypreParMatrix&>(*pH),MATHYPRE,
799  PETSC_USE_POINTER,A);
800  CCHKERRQ(pH->GetComm(),ierr);
801  }
802 #endif
803  else if (tid == PETSC_MATSHELL)
804  {
805  MakeWrapper(comm,&op,A);
806  }
807  else
808  {
809  MFEM_ABORT("Conversion from HypreParCSR to operator type = " << tid <<
810  " is not implemented");
811  }
812  }
813  else if (pB)
814  {
815  Mat *mats,*matsl2l = NULL;
816  PetscInt i,j,nr,nc;
817 
818  nr = pB->NumRowBlocks();
819  nc = pB->NumColBlocks();
820  ierr = PetscCalloc1(nr*nc,&mats); CCHKERRQ(PETSC_COMM_SELF,ierr);
821  if (tid == PETSC_MATIS)
822  {
823  ierr = PetscCalloc1(nr,&matsl2l); CCHKERRQ(PETSC_COMM_SELF,ierr);
824  }
825  for (i=0; i<nr; i++)
826  {
827  PetscBool needl2l = PETSC_TRUE;
828 
829  for (j=0; j<nc; j++)
830  {
831  if (!pB->IsZeroBlock(i,j))
832  {
833  ConvertOperator(comm,pB->GetBlock(i,j),&mats[i*nc+j],tid);
834  if (tid == PETSC_MATIS && needl2l)
835  {
836  PetscContainer c;
837  ierr = PetscObjectQuery((PetscObject)mats[i*nc+j],"__mfem_l2l",
838  (PetscObject*)&c);
839  PCHKERRQ(mats[i*nc+j],ierr);
840  // special case for block operators: the local Vdofs should be
841  // ordered as:
842  // [f1_1,...f1_N1,f2_1,...,f2_N2,...,fm_1,...,fm_Nm]
843  // with m fields, Ni the number of Vdofs for the i-th field
844  if (c)
845  {
846  Array<Mat> *l2l = NULL;
847  ierr = PetscContainerGetPointer(c,(void**)&l2l);
848  PCHKERRQ(c,ierr);
849  MFEM_VERIFY(l2l->Size() == 1,"Unexpected size "
850  << l2l->Size() << " for block row " << i );
851  ierr = PetscObjectReference((PetscObject)(*l2l)[0]);
852  PCHKERRQ(c,ierr);
853  matsl2l[i] = (*l2l)[0];
854  needl2l = PETSC_FALSE;
855  }
856  }
857  }
858  }
859  }
860  ierr = MatCreateNest(comm,nr,NULL,nc,NULL,mats,A); CCHKERRQ(comm,ierr);
861  if (tid == PETSC_MATIS)
862  {
863  ierr = MatConvert(*A,MATIS,MAT_INPLACE_MATRIX,A); CCHKERRQ(comm,ierr);
864 
865  mfem::Array<Mat> *vmatsl2l = new mfem::Array<Mat>(nr);
866  for (PetscInt i=0; i<nr; i++) { (*vmatsl2l)[i] = matsl2l[i]; }
867  ierr = PetscFree(matsl2l); CCHKERRQ(PETSC_COMM_SELF,ierr);
868 
869  PetscContainer c;
870  ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
871  ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
872  ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
873  PCHKERRQ(c,ierr);
874  ierr = PetscObjectCompose((PetscObject)(*A),"__mfem_l2l",(PetscObject)c);
875  PCHKERRQ((*A),ierr);
876  ierr = PetscContainerDestroy(&c); CCHKERRQ(comm,ierr);
877  }
878  for (i=0; i<nr*nc; i++) { ierr = MatDestroy(&mats[i]); CCHKERRQ(comm,ierr); }
879  ierr = PetscFree(mats); CCHKERRQ(PETSC_COMM_SELF,ierr);
880  }
881  else if (pI)
882  {
883  MFEM_VERIFY(tid == PETSC_MATAIJ,"Unsupported operation");
884  PetscInt rst;
885 
886  ierr = MatCreate(comm,A); CCHKERRQ(comm,ierr);
887  ierr = MatSetSizes(*A,pI->Height(),pI->Width(),PETSC_DECIDE,PETSC_DECIDE);
888  PCHKERRQ(A,ierr);
889  ierr = MatSetType(*A,MATAIJ); PCHKERRQ(*A,ierr);
890  ierr = MatMPIAIJSetPreallocation(*A,1,NULL,0,NULL); PCHKERRQ(*A,ierr);
891  ierr = MatSeqAIJSetPreallocation(*A,1,NULL); PCHKERRQ(*A,ierr);
892  ierr = MatSetOption(*A,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE); PCHKERRQ(*A,ierr);
893  ierr = MatGetOwnershipRange(*A,&rst,NULL); PCHKERRQ(*A,ierr);
894  for (PetscInt i = rst; i < rst+pI->Height(); i++)
895  {
896  ierr = MatSetValue(*A,i,i,1.,INSERT_VALUES); PCHKERRQ(*A,ierr);
897  }
898  ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY); PCHKERRQ(*A,ierr);
899  ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY); PCHKERRQ(*A,ierr);
900  }
901  else // fallback to general operator
902  {
903  MakeWrapper(comm,&op,A);
904  }
905 }
906 
907 void PetscParMatrix::Destroy()
908 {
909  if (A != NULL)
910  {
911  MPI_Comm comm = MPI_COMM_NULL;
912  ierr = PetscObjectGetComm((PetscObject)A,&comm); PCHKERRQ(A,ierr);
913  ierr = MatDestroy(&A); CCHKERRQ(comm,ierr);
914  }
915  delete X;
916  delete Y;
917  X = Y = NULL;
918 }
919 
920 PetscParMatrix::PetscParMatrix(Mat a, bool ref)
921 {
922  if (ref)
923  {
924  ierr = PetscObjectReference((PetscObject)a); PCHKERRQ(a,ierr);
925  }
926  Init();
927  A = a;
928  height = GetNumRows();
929  width = GetNumCols();
930 }
931 
932 // Computes y = alpha * A * x + beta * y
933 // or y = alpha * A^T* x + beta * y
934 static void MatMultKernel(Mat A,PetscScalar a,Vec X,PetscScalar b,Vec Y,
935  bool transpose)
936 {
937  PetscErrorCode (*f)(Mat,Vec,Vec);
938  PetscErrorCode (*fadd)(Mat,Vec,Vec,Vec);
939  if (transpose)
940  {
941  f = MatMultTranspose;
942  fadd = MatMultTransposeAdd;
943  }
944  else
945  {
946  f = MatMult;
947  fadd = MatMultAdd;
948  }
949  if (a != 0.)
950  {
951  if (b == 1.)
952  {
953  ierr = VecScale(X,a); PCHKERRQ(A,ierr);
954  ierr = (*fadd)(A,X,Y,Y); PCHKERRQ(A,ierr);
955  ierr = VecScale(X,1./a); PCHKERRQ(A,ierr);
956  }
957  else if (b != 0.)
958  {
959  ierr = VecScale(X,a); PCHKERRQ(A,ierr);
960  ierr = VecScale(Y,b); PCHKERRQ(A,ierr);
961  ierr = (*fadd)(A,X,Y,Y); PCHKERRQ(A,ierr);
962  ierr = VecScale(X,1./a); PCHKERRQ(A,ierr);
963  }
964  else
965  {
966  ierr = (*f)(A,X,Y); PCHKERRQ(A,ierr);
967  if (a != 1.)
968  {
969  ierr = VecScale(Y,a); PCHKERRQ(A,ierr);
970  }
971  }
972  }
973  else
974  {
975  if (b == 1.)
976  {
977  // do nothing
978  }
979  else if (b != 0.)
980  {
981  ierr = VecScale(Y,b); PCHKERRQ(A,ierr);
982  }
983  else
984  {
985  ierr = VecSet(Y,0.); PCHKERRQ(A,ierr);
986  }
987  }
988 }
989 
990 void PetscParMatrix::MakeRef(const PetscParMatrix &master)
991 {
992  ierr = PetscObjectReference((PetscObject)master.A); PCHKERRQ(master.A,ierr);
993  Destroy();
994  Init();
995  A = master.A;
996  height = master.height;
997  width = master.width;
998 }
999 
1000 PetscParVector * PetscParMatrix::GetX() const
1001 {
1002  if (!X)
1003  {
1004  MFEM_VERIFY(A,"Mat not present");
1005  X = new PetscParVector(*this,false); PCHKERRQ(A,ierr);
1006  }
1007  return X;
1008 }
1009 
1010 PetscParVector * PetscParMatrix::GetY() const
1011 {
1012  if (!Y)
1013  {
1014  MFEM_VERIFY(A,"Mat not present");
1015  Y = new PetscParVector(*this,true); PCHKERRQ(A,ierr);
1016  }
1017  return Y;
1018 }
1019 
1020 PetscParMatrix * PetscParMatrix::Transpose(bool action)
1021 {
1022  Mat B;
1023  if (action)
1024  {
1025  ierr = MatCreateTranspose(A,&B); PCHKERRQ(A,ierr);
1026  }
1027  else
1028  {
1029  ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&B); PCHKERRQ(A,ierr);
1030  }
1031  return new PetscParMatrix(B,false);
1032 }
1033 
1034 void PetscParMatrix::operator*=(double s)
1035 {
1036  ierr = MatScale(A,s); PCHKERRQ(A,ierr);
1037 }
1038 
1039 void PetscParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
1040 {
1041  MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
1042  << ", expected size = " << Width());
1043  MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
1044  << ", expected size = " << Height());
1045 
1046  PetscParVector *XX = GetX();
1047  PetscParVector *YY = GetY();
1048  XX->PlaceArray(x.GetData());
1049  YY->PlaceArray(y.GetData());
1050  MatMultKernel(A,a,XX->x,b,YY->x,false);
1051  XX->ResetArray();
1052  YY->ResetArray();
1053 }
1054 
1055 void PetscParMatrix::MultTranspose(double a, const Vector &x, double b,
1056  Vector &y) const
1057 {
1058  MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
1059  << ", expected size = " << Height());
1060  MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
1061  << ", expected size = " << Width());
1062 
1063  PetscParVector *XX = GetX();
1064  PetscParVector *YY = GetY();
1065  YY->PlaceArray(x.GetData());
1066  XX->PlaceArray(y.GetData());
1067  MatMultKernel(A,a,YY->x,b,XX->x,true);
1068  XX->ResetArray();
1069  YY->ResetArray();
1070 }
1071 
1072 void PetscParMatrix::Print(const char *fname, bool binary) const
1073 {
1074  if (fname)
1075  {
1076  PetscViewer view;
1077 
1078  if (binary)
1079  {
1080  ierr = PetscViewerBinaryOpen(PetscObjectComm((PetscObject)A),fname,
1081  FILE_MODE_WRITE,&view);
1082  }
1083  else
1084  {
1085  ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)A),fname,&view);
1086  }
1087  PCHKERRQ(A,ierr);
1088  ierr = MatView(A,view); PCHKERRQ(A,ierr);
1089  ierr = PetscViewerDestroy(&view); PCHKERRQ(A,ierr);
1090  }
1091  else
1092  {
1093  ierr = MatView(A,NULL); PCHKERRQ(A,ierr);
1094  }
1095 }
1096 
1097 
1098 PetscParMatrix * RAP(PetscParMatrix *Rt, PetscParMatrix *A, PetscParMatrix *P)
1099 {
1100  Mat pA = *A,pP = *P,pRt = *Rt;
1101  Mat B;
1102  PetscBool Aismatis,Pismatis,Rtismatis;
1103 
1104  MFEM_VERIFY(A->Width() == P->Height(),
1105  "Petsc RAP: Number of local cols of A " << A->Width() <<
1106  " differs from number of local rows of P " << P->Height());
1107  MFEM_VERIFY(A->Height() == Rt->Height(),
1108  "Petsc RAP: Number of local rows of A " << A->Height() <<
1109  " differs from number of local rows of Rt " << Rt->Height());
1110  ierr = PetscObjectTypeCompare((PetscObject)pA,MATIS,&Aismatis);
1111  PCHKERRQ(pA,ierr);
1112  ierr = PetscObjectTypeCompare((PetscObject)pP,MATIS,&Pismatis);
1113  PCHKERRQ(pA,ierr);
1114  ierr = PetscObjectTypeCompare((PetscObject)pRt,MATIS,&Rtismatis);
1115  PCHKERRQ(pA,ierr);
1116  if (Aismatis &&
1117  Pismatis &&
1118  Rtismatis) // handle special case (this code will eventually go into PETSc)
1119  {
1120  Mat lA,lP,lB,lRt;
1121  ISLocalToGlobalMapping cl2gP,cl2gRt;
1122  PetscInt rlsize,clsize,rsize,csize;
1123 
1124  ierr = MatGetLocalToGlobalMapping(pP,NULL,&cl2gP); PCHKERRQ(pA,ierr);
1125  ierr = MatGetLocalToGlobalMapping(pRt,NULL,&cl2gRt); PCHKERRQ(pA,ierr);
1126  ierr = MatGetLocalSize(pP,NULL,&clsize); PCHKERRQ(pP,ierr);
1127  ierr = MatGetLocalSize(pRt,NULL,&rlsize); PCHKERRQ(pRt,ierr);
1128  ierr = MatGetSize(pP,NULL,&csize); PCHKERRQ(pP,ierr);
1129  ierr = MatGetSize(pRt,NULL,&rsize); PCHKERRQ(pRt,ierr);
1130  ierr = MatCreate(A->GetComm(),&B); PCHKERRQ(pA,ierr);
1131  ierr = MatSetSizes(B,rlsize,clsize,rsize,csize); PCHKERRQ(B,ierr);
1132  ierr = MatSetType(B,MATIS); PCHKERRQ(B,ierr);
1133  ierr = MatSetLocalToGlobalMapping(B,cl2gRt,cl2gP); PCHKERRQ(B,ierr);
1134  ierr = MatISGetLocalMat(pA,&lA); PCHKERRQ(pA,ierr);
1135  ierr = MatISGetLocalMat(pP,&lP); PCHKERRQ(pA,ierr);
1136  ierr = MatISGetLocalMat(pRt,&lRt); PCHKERRQ(pA,ierr);
1137  if (lRt == lP)
1138  {
1139  ierr = MatPtAP(lA,lP,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&lB);
1140  PCHKERRQ(lA,ierr);
1141  }
1142  else
1143  {
1144  Mat lR;
1145  ierr = MatTranspose(lRt,MAT_INITIAL_MATRIX,&lR); PCHKERRQ(lRt,ierr);
1146  ierr = MatMatMatMult(lR,lA,lP,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&lB);
1147  PCHKERRQ(lRt,ierr);
1148  ierr = MatDestroy(&lR); PCHKERRQ(lRt,ierr);
1149  }
1150 
1151  // attach lRt matrix to the subdomain local matrix
1152  // it may be used if markers on vdofs have to be mapped on
1153  // subdomain true dofs
1154  {
1155  mfem::Array<Mat> *vmatsl2l = new mfem::Array<Mat>(1);
1156  ierr = PetscObjectReference((PetscObject)lRt); PCHKERRQ(lRt,ierr);
1157  (*vmatsl2l)[0] = lRt;
1158 
1159  PetscContainer c;
1160  ierr = PetscContainerCreate(PetscObjectComm((PetscObject)B),&c);
1161  PCHKERRQ(B,ierr);
1162  ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
1163  ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
1164  PCHKERRQ(c,ierr);
1165  ierr = PetscObjectCompose((PetscObject)B,"__mfem_l2l",(PetscObject)c);
1166  PCHKERRQ(B,ierr);
1167  ierr = PetscContainerDestroy(&c); PCHKERRQ(B,ierr);
1168  }
1169 
1170  // Set local problem
1171  ierr = MatISSetLocalMat(B,lB); PCHKERRQ(lB,ierr);
1172  ierr = MatDestroy(&lB); PCHKERRQ(lA,ierr);
1173  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); PCHKERRQ(B,ierr);
1174  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); PCHKERRQ(B,ierr);
1175  }
1176  else // it raises an error if the PtAP is not supported in PETSc
1177  {
1178  if (pP == pRt)
1179  {
1180  ierr = MatPtAP(pA,pP,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B);
1181  PCHKERRQ(pA,ierr);
1182  }
1183  else
1184  {
1185  Mat pR;
1186  ierr = MatTranspose(pRt,MAT_INITIAL_MATRIX,&pR); PCHKERRQ(Rt,ierr);
1187  ierr = MatMatMatMult(pR,pA,pP,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B);
1188  PCHKERRQ(pRt,ierr);
1189  ierr = MatDestroy(&pR); PCHKERRQ(pRt,ierr);
1190  }
1191  }
1192  return new PetscParMatrix(B);
1193 }
1194 
1195 PetscParMatrix * RAP(PetscParMatrix *A, PetscParMatrix *P)
1196 {
1197  PetscParMatrix *out = RAP(P,A,P);
1198  return out;
1199 }
1200 
1201 PetscParMatrix* PetscParMatrix::EliminateRowsCols(const Array<int> &rows_cols)
1202 {
1203  Mat Ae;
1204 
1205  PetscParVector dummy(GetComm(),0);
1206  ierr = MatDuplicate(A,MAT_COPY_VALUES,&Ae); PCHKERRQ(A,ierr);
1207  EliminateRowsCols(rows_cols,dummy,dummy);
1208  ierr = MatAXPY(Ae,-1.,A,SAME_NONZERO_PATTERN); PCHKERRQ(A,ierr);
1209  return new PetscParMatrix(Ae);
1210 }
1211 
1212 void PetscParMatrix::EliminateRowsCols(const Array<int> &rows_cols,
1213  const HypreParVector &X,
1214  HypreParVector &B)
1215 {
1216  MFEM_ABORT("PetscParMatrix::EliminateRowsCols(). To be implemented");
1217 }
1218 
1219 void PetscParMatrix::EliminateRowsCols(const Array<int> &rows_cols,
1220  const PetscParVector &X,
1221  PetscParVector &B)
1222 {
1223  PetscInt M,N;
1224  ierr = MatGetSize(A,&M,&N); PCHKERRQ(A,ierr);
1225  MFEM_VERIFY(M == N,"Rectangular case unsupported");
1226 
1227  // TODO: what if a diagonal term is not present?
1228  ierr = MatSetOption(A,MAT_NO_OFF_PROC_ZERO_ROWS,PETSC_TRUE); PCHKERRQ(A,ierr);
1229 
1230  // rows need to be in global numbering
1231  PetscInt rst;
1232  ierr = MatGetOwnershipRange(A,&rst,NULL); PCHKERRQ(A,ierr);
1233 
1234  IS dir;
1235  ierr = Convert_Array_IS(GetComm(),true,&rows_cols,rst,&dir); PCHKERRQ(A,ierr);
1236  if (!X.GlobalSize() && !B.GlobalSize())
1237  {
1238  ierr = MatZeroRowsColumnsIS(A,dir,1.,NULL,NULL); PCHKERRQ(A,ierr);
1239  }
1240  else
1241  {
1242  ierr = MatZeroRowsColumnsIS(A,dir,1.,X,B); PCHKERRQ(A,ierr);
1243  }
1244  ierr = ISDestroy(&dir); PCHKERRQ(A,ierr);
1245 }
1246 
1247 Mat PetscParMatrix::ReleaseMat(bool dereference)
1248 {
1249 
1250  Mat B = A;
1251  if (dereference)
1252  {
1253  MPI_Comm comm = GetComm();
1254  ierr = PetscObjectDereference((PetscObject)A); CCHKERRQ(comm,ierr);
1255  }
1256  A = NULL;
1257  return B;
1258 }
1259 
1260 Operator::Type PetscParMatrix::GetType() const
1261 {
1262  PetscBool ok;
1263  MFEM_VERIFY(A, "no associated PETSc Mat object");
1264  PetscObject oA = (PetscObject)(this->A);
1265  // map all of MATAIJ, MATSEQAIJ, and MATMPIAIJ to -> PETSC_MATAIJ
1266  ierr = PetscObjectTypeCompare(oA, MATAIJ, &ok); PCHKERRQ(A,ierr);
1267  if (ok == PETSC_TRUE) { return PETSC_MATAIJ; }
1268  ierr = PetscObjectTypeCompare(oA, MATSEQAIJ, &ok); PCHKERRQ(A,ierr);
1269  if (ok == PETSC_TRUE) { return PETSC_MATAIJ; }
1270  ierr = PetscObjectTypeCompare(oA, MATMPIAIJ, &ok); PCHKERRQ(A,ierr);
1271  if (ok == PETSC_TRUE) { return PETSC_MATAIJ; }
1272  ierr = PetscObjectTypeCompare(oA, MATIS, &ok); PCHKERRQ(A,ierr);
1273  if (ok == PETSC_TRUE) { return PETSC_MATIS; }
1274  ierr = PetscObjectTypeCompare(oA, MATSHELL, &ok); PCHKERRQ(A,ierr);
1275  if (ok == PETSC_TRUE) { return PETSC_MATSHELL; }
1276  ierr = PetscObjectTypeCompare(oA, MATNEST, &ok); PCHKERRQ(A,ierr);
1277  if (ok == PETSC_TRUE) { return PETSC_MATNEST; }
1278 #if defined(PETSC_HAVE_HYPRE)
1279  ierr = PetscObjectTypeCompare(oA, MATHYPRE, &ok); PCHKERRQ(A,ierr);
1280  if (ok == PETSC_TRUE) { return PETSC_MATHYPRE; }
1281 #endif
1282  return PETSC_MATGENERIC;
1283 }
1284 
1285 void EliminateBC(PetscParMatrix &A, PetscParMatrix &Ae,
1286  const Array<int> &ess_dof_list,
1287  const Vector &X, Vector &B)
1288 {
1289  const PetscScalar *array;
1290  Mat pA = const_cast<PetscParMatrix&>(A);
1291 
1292  // B -= Ae*X
1293  Ae.Mult(-1.0, X, 1.0, B);
1294 
1295  Vec diag = const_cast<PetscParVector&>((*A.GetX()));
1296  ierr = MatGetDiagonal(pA,diag); PCHKERRQ(pA,ierr);
1297  ierr = VecGetArrayRead(diag,&array); PCHKERRQ(diag,ierr);
1298  for (int i = 0; i < ess_dof_list.Size(); i++)
1299  {
1300  int r = ess_dof_list[i];
1301  B(r) = array[r] * X(r);
1302  }
1303  ierr = VecRestoreArrayRead(diag,&array); PCHKERRQ(diag,ierr);
1304 }
1305 
1306 // PetscSolver methods
1307 
1308 PetscSolver::PetscSolver() : clcustom(false)
1309 {
1310  obj = NULL;
1311  B = X = NULL;
1312  cid = -1;
1313  operatorset = false;
1314  bchandler = NULL;
1315  private_ctx = NULL;
1316 }
1317 
1318 PetscSolver::~PetscSolver()
1319 {
1320  delete B;
1321  delete X;
1322  FreePrivateContext();
1323 }
1324 
1325 void PetscSolver::SetTol(double tol)
1326 {
1327  SetRelTol(tol);
1328 }
1329 
1330 void PetscSolver::SetRelTol(double tol)
1331 {
1332  if (cid == KSP_CLASSID)
1333  {
1334  KSP ksp = (KSP)obj;
1335  ierr = KSPSetTolerances(ksp,tol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1336  }
1337  else if (cid == SNES_CLASSID)
1338  {
1339  SNES snes = (SNES)obj;
1340  ierr = SNESSetTolerances(snes,PETSC_DEFAULT,tol,PETSC_DEFAULT,PETSC_DEFAULT,
1341  PETSC_DEFAULT);
1342  }
1343  else if (cid == TS_CLASSID)
1344  {
1345  TS ts = (TS)obj;
1346  ierr = TSSetTolerances(ts,PETSC_DECIDE,NULL,tol,NULL);
1347  }
1348  else
1349  {
1350  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1351  }
1352  PCHKERRQ(obj,ierr);
1353 }
1354 
1355 void PetscSolver::SetAbsTol(double tol)
1356 {
1357  if (cid == KSP_CLASSID)
1358  {
1359  KSP ksp = (KSP)obj;
1360  ierr = KSPSetTolerances(ksp,PETSC_DEFAULT,tol,PETSC_DEFAULT,PETSC_DEFAULT);
1361  }
1362  else if (cid == SNES_CLASSID)
1363  {
1364  SNES snes = (SNES)obj;
1365  ierr = SNESSetTolerances(snes,tol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,
1366  PETSC_DEFAULT);
1367  }
1368  else if (cid == TS_CLASSID)
1369  {
1370  TS ts = (TS)obj;
1371  ierr = TSSetTolerances(ts,tol,NULL,PETSC_DECIDE,NULL);
1372  }
1373  else
1374  {
1375  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1376  }
1377  PCHKERRQ(obj,ierr);
1378 }
1379 
1380 void PetscSolver::SetMaxIter(int max_iter)
1381 {
1382  if (cid == KSP_CLASSID)
1383  {
1384  KSP ksp = (KSP)obj;
1385  ierr = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,
1386  max_iter);
1387  }
1388  else if (cid == SNES_CLASSID)
1389  {
1390  SNES snes = (SNES)obj;
1391  ierr = SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,
1392  max_iter,PETSC_DEFAULT);
1393  }
1394  else if (cid == TS_CLASSID)
1395  {
1396  TS ts = (TS)obj;
1397  ierr = TSSetMaxSteps(ts,max_iter);
1398  }
1399  else
1400  {
1401  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1402  }
1403  PCHKERRQ(obj,ierr);
1404 }
1405 
1406 
1407 void PetscSolver::SetPrintLevel(int plev)
1408 {
1409  typedef PetscErrorCode (*myPetscFunc)(void**);
1410  PetscViewerAndFormat *vf = NULL;
1411  PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm(obj));
1412 
1413  if (plev > 0)
1414  {
1415  ierr = PetscViewerAndFormatCreate(viewer,PETSC_VIEWER_DEFAULT,&vf);
1416  PCHKERRQ(obj,ierr);
1417  }
1418  if (cid == KSP_CLASSID)
1419  {
1420  // there are many other options, see the function KSPSetFromOptions() in
1421  // src/ksp/ksp/interface/itcl.c
1422  typedef PetscErrorCode (*myMonitor)(KSP,PetscInt,PetscReal,void*);
1423  KSP ksp = (KSP)obj;
1424  if (plev >= 0)
1425  {
1426  ierr = KSPMonitorCancel(ksp); PCHKERRQ(ksp,ierr);
1427  }
1428  if (plev == 1)
1429  {
1430  ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorDefault,vf,
1431  (myPetscFunc)PetscViewerAndFormatDestroy);
1432  PCHKERRQ(ksp,ierr);
1433  }
1434  else if (plev > 1)
1435  {
1436  ierr = KSPSetComputeSingularValues(ksp,PETSC_TRUE); PCHKERRQ(ksp,ierr);
1437  ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorSingularValue,vf,
1438  (myPetscFunc)PetscViewerAndFormatDestroy);
1439  PCHKERRQ(ksp,ierr);
1440  if (plev > 2)
1441  {
1442  ierr = PetscViewerAndFormatCreate(viewer,PETSC_VIEWER_DEFAULT,&vf);
1443  PCHKERRQ(viewer,ierr);
1444  ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidualNorm,vf,
1445  (myPetscFunc)PetscViewerAndFormatDestroy);
1446  PCHKERRQ(ksp,ierr);
1447  }
1448  }
1449  }
1450  else if (cid == SNES_CLASSID)
1451  {
1452  typedef PetscErrorCode (*myMonitor)(SNES,PetscInt,PetscReal,void*);
1453  SNES snes = (SNES)obj;
1454  if (plev >= 0)
1455  {
1456  ierr = SNESMonitorCancel(snes); PCHKERRQ(snes,ierr);
1457  }
1458  if (plev > 0)
1459  {
1460  ierr = SNESMonitorSet(snes,(myMonitor)SNESMonitorDefault,vf,
1461  (myPetscFunc)PetscViewerAndFormatDestroy);
1462  PCHKERRQ(snes,ierr);
1463  }
1464  }
1465  else if (cid == TS_CLASSID)
1466  {
1467  TS ts = (TS)obj;
1468  if (plev >= 0)
1469  {
1470  ierr = TSMonitorCancel(ts); PCHKERRQ(ts,ierr);
1471  }
1472  }
1473  else
1474  {
1475  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1476  }
1477 }
1478 
1479 void PetscSolver::SetMonitor(PetscSolverMonitor *ctx)
1480 {
1481  if (cid == KSP_CLASSID)
1482  {
1483  ierr = KSPMonitorSet((KSP)obj,__mfem_ksp_monitor,ctx,NULL);
1484  PCHKERRQ(obj,ierr);
1485  }
1486  else if (cid == SNES_CLASSID)
1487  {
1488  ierr = SNESMonitorSet((SNES)obj,__mfem_snes_monitor,ctx,NULL);
1489  PCHKERRQ(obj,ierr);
1490  }
1491  else if (cid == TS_CLASSID)
1492  {
1493  ierr = TSMonitorSet((TS)obj,__mfem_ts_monitor,ctx,NULL);
1494  PCHKERRQ(obj,ierr);
1495  }
1496  else
1497  {
1498  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1499  }
1500 }
1501 
1502 void PetscSolver::SetBCHandler(PetscBCHandler *bch)
1503 {
1504  bchandler = bch;
1505  if (cid == SNES_CLASSID)
1506  {
1507  __mfem_snes_ctx* snes_ctx = (__mfem_snes_ctx*)private_ctx;
1508  snes_ctx->bchandler = bchandler;
1509  }
1510  else if (cid == TS_CLASSID)
1511  {
1512  __mfem_ts_ctx* ts_ctx = (__mfem_ts_ctx*)private_ctx;
1513  ts_ctx->bchandler = bchandler;
1514  }
1515  else
1516  {
1517  MFEM_ABORT("Handling of essential bc only implemented for nonlinear and time-dependent solvers");
1518  }
1519 }
1520 
1521 void PetscSolver::SetPreconditionerFactory(PetscPreconditionerFactory *factory)
1522 {
1523  PC pc = NULL;
1524  if (cid == TS_CLASSID)
1525  {
1526  SNES snes;
1527  KSP ksp;
1528 
1529  ierr = TSGetSNES((TS)obj,&snes); PCHKERRQ(obj,ierr);
1530  ierr = SNESGetKSP(snes,&ksp); PCHKERRQ(obj,ierr);
1531  ierr = KSPGetPC(ksp,&pc); PCHKERRQ(obj,ierr);
1532  }
1533  else if (cid == SNES_CLASSID)
1534  {
1535  KSP ksp;
1536 
1537  ierr = SNESGetKSP((SNES)obj,&ksp); PCHKERRQ(obj,ierr);
1538  ierr = KSPGetPC(ksp,&pc); PCHKERRQ(obj,ierr);
1539  }
1540  else if (cid == KSP_CLASSID)
1541  {
1542  ierr = KSPGetPC((KSP)obj,&pc); PCHKERRQ(obj,ierr);
1543  }
1544  else if (cid == PC_CLASSID)
1545  {
1546  pc = (PC)obj;
1547  }
1548  else
1549  {
1550  MFEM_ABORT("No support for PetscPreconditionerFactory for this object");
1551  }
1552  ierr = MakeShellPCWithFactory(pc,factory); PCHKERRQ(pc,ierr);
1553 }
1554 
1555 void PetscSolver::Customize(bool customize) const
1556 {
1557  if (!customize) { clcustom = true; }
1558  if (!clcustom)
1559  {
1560  if (cid == PC_CLASSID)
1561  {
1562  PC pc = (PC)obj;
1563  ierr = PCSetFromOptions(pc); PCHKERRQ(pc, ierr);
1564  }
1565  else if (cid == KSP_CLASSID)
1566  {
1567  KSP ksp = (KSP)obj;
1568  ierr = KSPSetFromOptions(ksp); PCHKERRQ(ksp, ierr);
1569  }
1570  else if (cid == SNES_CLASSID)
1571  {
1572  SNES snes = (SNES)obj;
1573  ierr = SNESSetFromOptions(snes); PCHKERRQ(snes, ierr);
1574  }
1575  else if (cid == TS_CLASSID)
1576  {
1577  TS ts = (TS)obj;
1578  ierr = TSSetFromOptions(ts); PCHKERRQ(ts, ierr);
1579  }
1580  else
1581  {
1582  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1583  }
1584  }
1585  clcustom = true;
1586 }
1587 
1588 int PetscSolver::GetConverged()
1589 {
1590  if (cid == KSP_CLASSID)
1591  {
1592  KSP ksp = (KSP)obj;
1593  KSPConvergedReason reason;
1594  ierr = KSPGetConvergedReason(ksp,&reason);
1595  PCHKERRQ(ksp,ierr);
1596  return reason > 0 ? 1 : 0;
1597  }
1598  else if (cid == SNES_CLASSID)
1599  {
1600  SNES snes = (SNES)obj;
1601  SNESConvergedReason reason;
1602  ierr = SNESGetConvergedReason(snes,&reason);
1603  PCHKERRQ(snes,ierr);
1604  return reason > 0 ? 1 : 0;
1605  }
1606  else if (cid == TS_CLASSID)
1607  {
1608  TS ts = (TS)obj;
1609  TSConvergedReason reason;
1610  ierr = TSGetConvergedReason(ts,&reason);
1611  PCHKERRQ(ts,ierr);
1612  return reason > 0 ? 1 : 0;
1613  }
1614  else
1615  {
1616  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1617  return -1;
1618  }
1619 }
1620 
1621 int PetscSolver::GetNumIterations()
1622 {
1623  if (cid == KSP_CLASSID)
1624  {
1625  KSP ksp = (KSP)obj;
1626  PetscInt its;
1627  ierr = KSPGetIterationNumber(ksp,&its);
1628  PCHKERRQ(ksp,ierr);
1629  return its;
1630  }
1631  else if (cid == SNES_CLASSID)
1632  {
1633  SNES snes = (SNES)obj;
1634  PetscInt its;
1635  ierr = SNESGetIterationNumber(snes,&its);
1636  PCHKERRQ(snes,ierr);
1637  return its;
1638  }
1639  else if (cid == TS_CLASSID)
1640  {
1641  TS ts = (TS)obj;
1642  PetscInt its;
1643  ierr = TSGetStepNumber(ts,&its);
1644  PCHKERRQ(ts,ierr);
1645  return its;
1646  }
1647  else
1648  {
1649  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1650  return -1;
1651  }
1652 }
1653 
1654 double PetscSolver::GetFinalNorm()
1655 {
1656  if (cid == KSP_CLASSID)
1657  {
1658  KSP ksp = (KSP)obj;
1659  PetscReal norm;
1660  ierr = KSPGetResidualNorm(ksp,&norm);
1661  PCHKERRQ(ksp,ierr);
1662  return norm;
1663  }
1664  if (cid == SNES_CLASSID)
1665  {
1666  SNES snes = (SNES)obj;
1667  PetscReal norm;
1668  ierr = SNESGetFunctionNorm(snes,&norm);
1669  PCHKERRQ(snes,ierr);
1670  return norm;
1671  }
1672  else
1673  {
1674  MFEM_ABORT("CLASSID = " << cid << " is not implemented!");
1675  return PETSC_MAX_REAL;
1676  }
1677 }
1678 
1679 void PetscSolver::CreatePrivateContext()
1680 {
1681  FreePrivateContext();
1682  if (cid == SNES_CLASSID)
1683  {
1684  __mfem_snes_ctx *snes_ctx;
1685  ierr = PetscNew(&snes_ctx); CCHKERRQ(PETSC_COMM_SELF,ierr);
1686  snes_ctx->op = NULL;
1687  snes_ctx->bchandler = NULL;
1688  snes_ctx->work = NULL;
1689  snes_ctx->jacType = Operator::PETSC_MATAIJ;
1690  private_ctx = (void*) snes_ctx;
1691  }
1692  else if (cid == TS_CLASSID)
1693  {
1694  __mfem_ts_ctx *ts_ctx;
1695  ierr = PetscNew(&ts_ctx); CCHKERRQ(PETSC_COMM_SELF,ierr);
1696  ts_ctx->op = NULL;
1697  ts_ctx->bchandler = NULL;
1698  ts_ctx->work = NULL;
1699  ts_ctx->cached_shift = std::numeric_limits<PetscReal>::min();
1700  ts_ctx->cached_ijacstate = -1;
1701  ts_ctx->cached_rhsjacstate = -1;
1702  ts_ctx->type = PetscODESolver::ODE_SOLVER_GENERAL;
1703  ts_ctx->jacType = Operator::PETSC_MATAIJ;
1704  private_ctx = (void*) ts_ctx;
1705  }
1706 }
1707 
1708 void PetscSolver::FreePrivateContext()
1709 {
1710  if (!private_ctx) { return; }
1711  // free private context's owned objects
1712  if (cid == SNES_CLASSID)
1713  {
1714  __mfem_snes_ctx *snes_ctx = (__mfem_snes_ctx *)private_ctx;
1715  delete snes_ctx->work;
1716  }
1717  else if (cid == TS_CLASSID)
1718  {
1719  __mfem_ts_ctx *ts_ctx = (__mfem_ts_ctx *)private_ctx;
1720  delete ts_ctx->work;
1721  }
1722  ierr = PetscFree(private_ctx); CCHKERRQ(PETSC_COMM_SELF,ierr);
1723 }
1724 
1725 // PetscBCHandler methods
1726 
1727 PetscBCHandler::PetscBCHandler(Array<int>& ess_tdof_list,
1728  enum PetscBCHandler::Type _type)
1729  : bctype(_type), setup(false), eval_t(0.0),
1730  eval_t_cached(std::numeric_limits<double>::min())
1731 {
1732  SetTDofs(ess_tdof_list);
1733 }
1734 
1735 void PetscBCHandler::SetTDofs(Array<int>& list)
1736 {
1737  ess_tdof_list.SetSize(list.Size());
1738  ess_tdof_list.Assign(list);
1739  setup = false;
1740 }
1741 
1742 void PetscBCHandler::SetUp(PetscInt n)
1743 {
1744  if (setup) { return; }
1745  if (bctype == CONSTANT)
1746  {
1747  eval_g.SetSize(n);
1748  this->Eval(eval_t,eval_g);
1749  eval_t_cached = eval_t;
1750  }
1751  else if (bctype == TIME_DEPENDENT)
1752  {
1753  eval_g.SetSize(n);
1754  }
1755  setup = true;
1756 }
1757 
1758 void PetscBCHandler::ApplyBC(const Vector &x, Vector &y)
1759 {
1760  if (!setup) { MFEM_ABORT("PetscBCHandler not yet setup"); }
1761  y = x;
1762  if (bctype == ZERO)
1763  {
1764  for (PetscInt i = 0; i < ess_tdof_list.Size(); ++i)
1765  {
1766  y[ess_tdof_list[i]] = 0.0;
1767  }
1768  }
1769  else
1770  {
1771  if (bctype != CONSTANT && eval_t != eval_t_cached)
1772  {
1773  Eval(eval_t,eval_g);
1774  eval_t_cached = eval_t;
1775  }
1776  for (PetscInt i = 0; i < ess_tdof_list.Size(); ++i)
1777  {
1778  y[ess_tdof_list[i]] = eval_g[ess_tdof_list[i]];
1779  }
1780  }
1781 }
1782 
1783 void PetscBCHandler::FixResidualBC(const Vector& x, Vector& y)
1784 {
1785  if (!setup) { MFEM_ABORT("PetscBCHandler not yet setup"); }
1786  if (bctype == ZERO)
1787  {
1788  for (PetscInt i = 0; i < ess_tdof_list.Size(); ++i)
1789  {
1790  y[ess_tdof_list[i]] = x[ess_tdof_list[i]];
1791  }
1792  }
1793  else
1794  {
1795  for (PetscInt i = 0; i < ess_tdof_list.Size(); ++i)
1796  {
1797  y[ess_tdof_list[i]] = x[ess_tdof_list[i]] - eval_g[ess_tdof_list[i]];
1798  }
1799  }
1800 }
1801 
1802 // PetscLinearSolver methods
1803 
1804 PetscLinearSolver::PetscLinearSolver(MPI_Comm comm, const std::string &prefix,
1805  bool wrapin)
1806  : PetscSolver(), Solver(), wrap(wrapin)
1807 {
1808  KSP ksp;
1809  ierr = KSPCreate(comm,&ksp); CCHKERRQ(comm,ierr);
1810  obj = (PetscObject)ksp;
1811  ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
1812  ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
1813 }
1814 
1815 PetscLinearSolver::PetscLinearSolver(const PetscParMatrix &A,
1816  const std::string &prefix)
1817  : PetscSolver(), Solver(), wrap(false)
1818 {
1819  KSP ksp;
1820  ierr = KSPCreate(A.GetComm(),&ksp); CCHKERRQ(A.GetComm(),ierr);
1821  obj = (PetscObject)ksp;
1822  ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
1823  ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
1824  SetOperator(A);
1825 }
1826 
1827 PetscLinearSolver::PetscLinearSolver(const HypreParMatrix &A, bool wrapin,
1828  const std::string &prefix)
1829  : PetscSolver(), Solver(), wrap(wrapin)
1830 {
1831  KSP ksp;
1832  ierr = KSPCreate(A.GetComm(),&ksp); CCHKERRQ(A.GetComm(),ierr);
1833  obj = (PetscObject)ksp;
1834  ierr = PetscObjectGetClassId(obj, &cid); PCHKERRQ(obj, ierr);
1835  ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
1836  SetOperator(A);
1837 }
1838 
1839 void PetscLinearSolver::SetOperator(const Operator &op)
1840 {
1841  const HypreParMatrix *hA = dynamic_cast<const HypreParMatrix *>(&op);
1842  PetscParMatrix *pA = const_cast<PetscParMatrix *>
1843  (dynamic_cast<const PetscParMatrix *>(&op));
1844  const Operator *oA = dynamic_cast<const Operator *>(&op);
1845 
1846  // Preserve Pmat if already set
1847  KSP ksp = (KSP)obj;
1848  Mat P = NULL;
1849  PetscBool pmat;
1850  ierr = KSPGetOperatorsSet(ksp,NULL,&pmat); PCHKERRQ(ksp,ierr);
1851  if (pmat)
1852  {
1853  ierr = KSPGetOperators(ksp,NULL,&P); PCHKERRQ(ksp,ierr);
1854  ierr = PetscObjectReference((PetscObject)P); PCHKERRQ(ksp,ierr);
1855  }
1856 
1857  // update base classes: Operator, Solver, PetscLinearSolver
1858  bool delete_pA = false;
1859  if (!pA)
1860  {
1861  if (hA)
1862  {
1863  // Create MATSHELL object or convert into a format suitable to construct preconditioners
1864  pA = new PetscParMatrix(hA, wrap ? PETSC_MATSHELL : PETSC_MATAIJ);
1865  delete_pA = true;
1866  }
1867  else if (oA) // fallback to general operator
1868  {
1869  // Create MATSHELL or MATNEST (if oA is a BlockOperator) object
1870  // If oA is a BlockOperator, Operator::Type is relevant to the subblocks
1871  pA = new PetscParMatrix(PetscObjectComm(obj),oA,
1872  wrap ? PETSC_MATSHELL : PETSC_MATAIJ);
1873  delete_pA = true;
1874  }
1875  }
1876  MFEM_VERIFY(pA, "Unsupported operation!");
1877 
1878  // Set operators into PETSc KSP
1879  Mat A = pA->A;
1880  if (operatorset)
1881  {
1882  Mat C;
1883  PetscInt nheight,nwidth,oheight,owidth;
1884 
1885  ierr = KSPGetOperators(ksp,&C,NULL); PCHKERRQ(ksp,ierr);
1886  ierr = MatGetSize(A,&nheight,&nwidth); PCHKERRQ(A,ierr);
1887  ierr = MatGetSize(C,&oheight,&owidth); PCHKERRQ(A,ierr);
1888  if (nheight != oheight || nwidth != owidth)
1889  {
1890  // reinit without destroying the KSP
1891  // communicator remains the same
1892  ierr = KSPReset(ksp); PCHKERRQ(ksp,ierr);
1893  delete X;
1894  delete B;
1895  X = B = NULL;
1896  wrap = false;
1897  }
1898  }
1899  if (P)
1900  {
1901  ierr = KSPSetOperators(ksp,A,P); PCHKERRQ(ksp,ierr);
1902  ierr = MatDestroy(&P); PCHKERRQ(ksp,ierr);
1903  }
1904  else
1905  {
1906  ierr = KSPSetOperators(ksp,A,A); PCHKERRQ(ksp,ierr);
1907  }
1908 
1909  // Update PetscSolver
1910  operatorset = true;
1911 
1912  // Update the Operator fields.
1913  height = pA->Height();
1914  width = pA->Width();
1915 
1916  if (delete_pA) { delete pA; }
1917 }
1918 
1919 void PetscLinearSolver::SetOperator(const Operator &op, const Operator &pop)
1920 {
1921  const HypreParMatrix *hA = dynamic_cast<const HypreParMatrix *>(&op);
1922  PetscParMatrix *pA = const_cast<PetscParMatrix *>
1923  (dynamic_cast<const PetscParMatrix *>(&op));
1924  const Operator *oA = dynamic_cast<const Operator *>(&op);
1925 
1926  PetscParMatrix *ppA = const_cast<PetscParMatrix *>
1927  (dynamic_cast<const PetscParMatrix *>(&pop));
1928  const Operator *poA = dynamic_cast<const Operator *>(&pop);
1929 
1930  // Convert Operator for linear system
1931  bool delete_pA = false;
1932  if (!pA)
1933  {
1934  if (hA)
1935  {
1936  // Create MATSHELL object or convert into a format suitable to construct preconditioners
1937  pA = new PetscParMatrix(hA, wrap ? PETSC_MATSHELL : PETSC_MATAIJ);
1938  delete_pA = true;
1939  }
1940  else if (oA) // fallback to general operator
1941  {
1942  // Create MATSHELL or MATNEST (if oA is a BlockOperator) object
1943  // If oA is a BlockOperator, Operator::Type is relevant to the subblocks
1944  pA = new PetscParMatrix(PetscObjectComm(obj),oA,
1945  wrap ? PETSC_MATSHELL : PETSC_MATAIJ);
1946  delete_pA = true;
1947  }
1948  }
1949  MFEM_VERIFY(pA, "Unsupported operation!");
1950 
1951  // Convert Operator to be preconditioned
1952  bool delete_ppA = false;
1953  if (!ppA)
1954  {
1955  if (oA == poA && !wrap) // Same operator, already converted
1956  {
1957  ppA = pA;
1958  }
1959  else
1960  {
1961  ppA = new PetscParMatrix(PetscObjectComm(obj), poA, PETSC_MATAIJ);
1962  delete_ppA = true;
1963  }
1964  }
1965  MFEM_VERIFY(ppA, "Unsupported operation!");
1966 
1967  // Set operators into PETSc KSP
1968  KSP ksp = (KSP)obj;
1969  Mat A = pA->A;
1970  Mat P = ppA->A;
1971  if (operatorset)
1972  {
1973  Mat C;
1974  PetscInt nheight,nwidth,oheight,owidth;
1975 
1976  ierr = KSPGetOperators(ksp,&C,NULL); PCHKERRQ(ksp,ierr);
1977  ierr = MatGetSize(A,&nheight,&nwidth); PCHKERRQ(A,ierr);
1978  ierr = MatGetSize(C,&oheight,&owidth); PCHKERRQ(A,ierr);
1979  if (nheight != oheight || nwidth != owidth)
1980  {
1981  // reinit without destroying the KSP
1982  // communicator remains the same
1983  ierr = KSPReset(ksp); PCHKERRQ(ksp,ierr);
1984  delete X;
1985  delete B;
1986  X = B = NULL;
1987  wrap = false;
1988  }
1989  }
1990  ierr = KSPSetOperators(ksp,A,P); PCHKERRQ(ksp,ierr);
1991 
1992  // Update PetscSolver
1993  operatorset = true;
1994 
1995  // Update the Operator fields.
1996  height = pA->Height();
1997  width = pA->Width();
1998 
1999  if (delete_pA) { delete pA; }
2000  if (delete_ppA) { delete ppA; }
2001 }
2002 
2003 void PetscLinearSolver::SetPreconditioner(Solver &precond)
2004 {
2005  KSP ksp = (KSP)obj;
2006 
2007  // Preserve Amat if already set
2008  Mat A = NULL;
2009  PetscBool amat;
2010  ierr = KSPGetOperatorsSet(ksp,&amat,NULL); PCHKERRQ(ksp,ierr);
2011  if (amat)
2012  {
2013  ierr = KSPGetOperators(ksp,&A,NULL); PCHKERRQ(ksp,ierr);
2014  ierr = PetscObjectReference((PetscObject)A); PCHKERRQ(ksp,ierr);
2015  }
2016  PetscPreconditioner *ppc = dynamic_cast<PetscPreconditioner *>(&precond);
2017  if (ppc)
2018  {
2019  ierr = KSPSetPC(ksp,*ppc); PCHKERRQ(ksp,ierr);
2020  }
2021  else
2022  {
2023  // wrap the Solver action
2024  // Solver is assumed to be already setup
2025  // ownership of precond is not tranferred,
2026  // consistently with other MFEM's linear solvers
2027  PC pc;
2028  ierr = KSPGetPC(ksp,&pc); PCHKERRQ(ksp,ierr);
2029  ierr = MakeShellPC(pc,precond,false); PCHKERRQ(ksp,ierr);
2030  }
2031  if (A)
2032  {
2033  Mat P;
2034 
2035  ierr = KSPGetOperators(ksp,NULL,&P); PCHKERRQ(ksp,ierr);
2036  ierr = PetscObjectReference((PetscObject)P); PCHKERRQ(ksp,ierr);
2037  ierr = KSPSetOperators(ksp,A,P); PCHKERRQ(ksp,ierr);
2038  ierr = MatDestroy(&A); PCHKERRQ(ksp,ierr);
2039  ierr = MatDestroy(&P); PCHKERRQ(ksp,ierr);
2040  }
2041 }
2042 
2043 void PetscLinearSolver::Mult(const Vector &b, Vector &x) const
2044 {
2045  KSP ksp = (KSP)obj;
2046 
2047  if (!B || !X)
2048  {
2049  Mat pA = NULL;
2050  ierr = KSPGetOperators(ksp, &pA, NULL); PCHKERRQ(obj, ierr);
2051  if (!B)
2052  {
2053  PetscParMatrix A = PetscParMatrix(pA, true);
2054  B = new PetscParVector(A, true, false);
2055  }
2056  if (!X)
2057  {
2058  PetscParMatrix A = PetscParMatrix(pA, true);
2059  X = new PetscParVector(A, false, false);
2060  }
2061  }
2062  B->PlaceArray(b.GetData());
2063  X->PlaceArray(x.GetData());
2064 
2065  Customize();
2066 
2067  ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
2068  PCHKERRQ(ksp, ierr);
2069 
2070  // Solve the system.
2071  ierr = KSPSolve(ksp, B->x, X->x); PCHKERRQ(ksp,ierr);
2072  B->ResetArray();
2073  X->ResetArray();
2074 }
2075 
2076 PetscLinearSolver::~PetscLinearSolver()
2077 {
2078  MPI_Comm comm;
2079  KSP ksp = (KSP)obj;
2080  ierr = PetscObjectGetComm((PetscObject)ksp,&comm); PCHKERRQ(ksp,ierr);
2081  ierr = KSPDestroy(&ksp); CCHKERRQ(comm,ierr);
2082 }
2083 
2084 // PetscPCGSolver methods
2085 
2086 PetscPCGSolver::PetscPCGSolver(MPI_Comm comm, const std::string &prefix)
2087  : PetscLinearSolver(comm,prefix)
2088 {
2089  KSP ksp = (KSP)obj;
2090  ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
2091  // this is to obtain a textbook PCG
2092  ierr = KSPSetNormType(ksp,KSP_NORM_NATURAL); PCHKERRQ(ksp,ierr);
2093 }
2094 
2095 PetscPCGSolver::PetscPCGSolver(PetscParMatrix& A, const std::string &prefix)
2096  : PetscLinearSolver(A,prefix)
2097 {
2098  KSP ksp = (KSP)obj;
2099  ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
2100  // this is to obtain a textbook PCG
2101  ierr = KSPSetNormType(ksp,KSP_NORM_NATURAL); PCHKERRQ(ksp,ierr);
2102 }
2103 
2104 PetscPCGSolver::PetscPCGSolver(HypreParMatrix& A, bool wrap,
2105  const std::string &prefix)
2106  : PetscLinearSolver(A,wrap,prefix)
2107 {
2108  KSP ksp = (KSP)obj;
2109  ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
2110  // this is to obtain a textbook PCG
2111  ierr = KSPSetNormType(ksp,KSP_NORM_NATURAL); PCHKERRQ(ksp,ierr);
2112 }
2113 
2114 // PetscPreconditioner methods
2115 
2116 PetscPreconditioner::PetscPreconditioner(MPI_Comm comm,
2117  const std::string &prefix)
2118  : PetscSolver(), Solver()
2119 {
2120  PC pc;
2121  ierr = PCCreate(comm,&pc); CCHKERRQ(comm,ierr);
2122  obj = (PetscObject)pc;
2123  ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
2124  ierr = PCSetOptionsPrefix(pc, prefix.c_str()); PCHKERRQ(pc, ierr);
2125 }
2126 
2127 PetscPreconditioner::PetscPreconditioner(PetscParMatrix &A,
2128  const string &prefix)
2129  : PetscSolver(), Solver()
2130 {
2131  PC pc;
2132  ierr = PCCreate(A.GetComm(),&pc); CCHKERRQ(A.GetComm(),ierr);
2133  obj = (PetscObject)pc;
2134  ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
2135  ierr = PCSetOptionsPrefix(pc, prefix.c_str()); PCHKERRQ(pc, ierr);
2136  SetOperator(A);
2137 }
2138 
2139 PetscPreconditioner::PetscPreconditioner(MPI_Comm comm, Operator &op,
2140  const string &prefix)
2141  : PetscSolver(), Solver()
2142 {
2143  PC pc;
2144  ierr = PCCreate(comm,&pc); CCHKERRQ(comm,ierr);
2145  obj = (PetscObject)pc;
2146  ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
2147  ierr = PCSetOptionsPrefix(pc, prefix.c_str()); PCHKERRQ(pc, ierr);
2148  SetOperator(op);
2149 }
2150 
2151 void PetscPreconditioner::SetOperator(const Operator &op)
2152 {
2153  bool delete_pA = false;
2154  PetscParMatrix *pA = const_cast<PetscParMatrix *>
2155  (dynamic_cast<const PetscParMatrix *>(&op));
2156 
2157  if (!pA)
2158  {
2159  const Operator *cop = dynamic_cast<const Operator *>(&op);
2160  pA = new PetscParMatrix(PetscObjectComm(obj),cop,PETSC_MATAIJ);
2161  delete_pA = true;
2162  }
2163 
2164  // Set operators into PETSc PC
2165  PC pc = (PC)obj;
2166  Mat A = pA->A;
2167  if (operatorset)
2168  {
2169  Mat C;
2170  PetscInt nheight,nwidth,oheight,owidth;
2171 
2172  ierr = PCGetOperators(pc,&C,NULL); PCHKERRQ(pc,ierr);
2173  ierr = MatGetSize(A,&nheight,&nwidth); PCHKERRQ(A,ierr);
2174  ierr = MatGetSize(C,&oheight,&owidth); PCHKERRQ(A,ierr);
2175  if (nheight != oheight || nwidth != owidth)
2176  {
2177  // reinit without destroying the PC
2178  // communicator remains the same
2179  ierr = PCReset(pc); PCHKERRQ(pc,ierr);
2180  delete X;
2181  delete B;
2182  X = B = NULL;
2183  }
2184  }
2185  ierr = PCSetOperators(pc,pA->A,pA->A); PCHKERRQ(obj,ierr);
2186 
2187  // Update PetscSolver
2188  operatorset = true;
2189 
2190  // Update the Operator fields.
2191  height = pA->Height();
2192  width = pA->Width();
2193 
2194  if (delete_pA) { delete pA; };
2195 }
2196 
2197 void PetscPreconditioner::Mult(const Vector &b, Vector &x) const
2198 {
2199  PC pc = (PC)obj;
2200 
2201  if (!B || !X)
2202  {
2203  Mat pA = NULL;
2204  ierr = PCGetOperators(pc, NULL, &pA); PCHKERRQ(obj, ierr);
2205  if (!B)
2206  {
2207  PetscParMatrix A(pA, true);
2208  B = new PetscParVector(A, true, false);
2209  }
2210  if (!X)
2211  {
2212  PetscParMatrix A(pA, true);
2213  X = new PetscParVector(A, false, false);
2214  }
2215  }
2216  B->PlaceArray(b.GetData());
2217  X->PlaceArray(x.GetData());
2218 
2219  Customize();
2220 
2221  // Apply the preconditioner.
2222  ierr = PCApply(pc, B->x, X->x); PCHKERRQ(pc, ierr);
2223  B->ResetArray();
2224  X->ResetArray();
2225 }
2226 
2227 PetscPreconditioner::~PetscPreconditioner()
2228 {
2229  MPI_Comm comm;
2230  PC pc = (PC)obj;
2231  ierr = PetscObjectGetComm((PetscObject)pc,&comm); PCHKERRQ(pc,ierr);
2232  ierr = PCDestroy(&pc); CCHKERRQ(comm,ierr);
2233 }
2234 
2235 // PetscBDDCSolver methods
2236 
2237 void PetscBDDCSolver::BDDCSolverConstructor(const PetscBDDCSolverParams &opts)
2238 {
2239  MPI_Comm comm = PetscObjectComm(obj);
2240 
2241  // get PETSc object
2242  PC pc = (PC)obj;
2243  Mat pA;
2244  ierr = PCGetOperators(pc,NULL,&pA); PCHKERRQ(pc,ierr);
2245 
2246  // matrix type should be of type MATIS
2247  PetscBool ismatis;
2248  ierr = PetscObjectTypeCompare((PetscObject)pA,MATIS,&ismatis);
2249  PCHKERRQ(pA,ierr);
2250  MFEM_VERIFY(ismatis,"PetscBDDCSolver needs the matrix in unassembled format");
2251 
2252  // set PETSc PC type to PCBDDC
2253  ierr = PCSetType(pc,PCBDDC); PCHKERRQ(obj,ierr);
2254 
2255  // index sets for fields splitting
2256  IS *fields = NULL;
2257  PetscInt nf = 0;
2258 
2259  // index sets for boundary dofs specification (Essential = dir, Natural = neu)
2260  IS dir = NULL, neu = NULL;
2261  PetscInt rst;
2262 
2263  // Extract l2l matrices
2264  Array<Mat> *l2l = NULL;
2265  if (opts.ess_dof_local || opts.nat_dof_local)
2266  {
2267  PetscContainer c;
2268 
2269  ierr = PetscObjectQuery((PetscObject)pA,"__mfem_l2l",(PetscObject*)&c);
2270  MFEM_VERIFY(c,"Local-to-local PETSc container not present");
2271  ierr = PetscContainerGetPointer(c,(void**)&l2l); PCHKERRQ(c,ierr);
2272  }
2273 
2274  // check information about index sets (essential dofs, fields, etc.)
2275 #ifdef MFEM_DEBUG
2276  {
2277  // make sure ess/nat_dof have been collectively set
2278  PetscBool lpr = PETSC_FALSE,pr;
2279  if (opts.ess_dof) { lpr = PETSC_TRUE; }
2280  ierr = MPI_Allreduce(&lpr,&pr,1,MPIU_BOOL,MPI_LOR,comm);
2281  PCHKERRQ(pA,ierr);
2282  MFEM_VERIFY(lpr == pr,"ess_dof should be collectively set");
2283  lpr = PETSC_FALSE;
2284  if (opts.nat_dof) { lpr = PETSC_TRUE; }
2285  ierr = MPI_Allreduce(&lpr,&pr,1,MPIU_BOOL,MPI_LOR,comm);
2286  PCHKERRQ(pA,ierr);
2287  MFEM_VERIFY(lpr == pr,"nat_dof should be collectively set");
2288  // make sure fields have been collectively set
2289  PetscInt ms[2],Ms[2];
2290  ms[0] = -nf; ms[1] = nf;
2291  ierr = MPI_Allreduce(&ms,&Ms,2,MPIU_INT,MPI_MAX,comm);
2292  PCHKERRQ(pA,ierr);
2293  MFEM_VERIFY(-Ms[0] == Ms[1],
2294  "number of fields should be the same across processes");
2295  }
2296 #endif
2297 
2298  // boundary sets
2299  ierr = MatGetOwnershipRange(pA,&rst,NULL); PCHKERRQ(pA,ierr);
2300  if (opts.ess_dof)
2301  {
2302  PetscInt st = opts.ess_dof_local ? 0 : rst;
2303  if (!opts.ess_dof_local)
2304  {
2305  // need to compute the boundary dofs in global ordering
2306  ierr = Convert_Array_IS(comm,true,opts.ess_dof,st,&dir);
2307  CCHKERRQ(comm,ierr);
2308  ierr = PCBDDCSetDirichletBoundaries(pc,dir); PCHKERRQ(pc,ierr);
2309  }
2310  else
2311  {
2312  // need to compute a list for the marked boundary dofs in local ordering
2313  ierr = Convert_Vmarks_IS(comm,*l2l,opts.ess_dof,st,&dir);
2314  CCHKERRQ(comm,ierr);
2315  ierr = PCBDDCSetDirichletBoundariesLocal(pc,dir); PCHKERRQ(pc,ierr);
2316  }
2317  }
2318  if (opts.nat_dof)
2319  {
2320  PetscInt st = opts.nat_dof_local ? 0 : rst;
2321  if (!opts.nat_dof_local)
2322  {
2323  // need to compute the boundary dofs in global ordering
2324  ierr = Convert_Array_IS(comm,true,opts.nat_dof,st,&neu);
2325  CCHKERRQ(comm,ierr);
2326  ierr = PCBDDCSetNeumannBoundaries(pc,neu); PCHKERRQ(pc,ierr);
2327  }
2328  else
2329  {
2330  // need to compute a list for the marked boundary dofs in local ordering
2331  ierr = Convert_Vmarks_IS(comm,*l2l,opts.nat_dof,st,&neu);
2332  CCHKERRQ(comm,ierr);
2333  ierr = PCBDDCSetNeumannBoundariesLocal(pc,neu); PCHKERRQ(pc,ierr);
2334  }
2335  }
2336 
2337  // field splitting
2338  if (nf)
2339  {
2340  ierr = PCBDDCSetDofsSplitting(pc,nf,fields); PCHKERRQ(pc,ierr);
2341  for (int i = 0; i < nf; i++)
2342  {
2343  ierr = ISDestroy(&fields[i]); CCHKERRQ(comm,ierr);
2344  }
2345  ierr = PetscFree(fields); PCHKERRQ(pc,ierr);
2346  }
2347 
2348  // code for block size is disabled since we cannot change the matrix
2349  // block size after it has been setup
2350  // int bs = 1;
2351 
2352  // Customize using the finite element space (if any)
2353  ParFiniteElementSpace *fespace = opts.fespace;
2354  if (fespace)
2355  {
2356  const FiniteElementCollection *fec = fespace->FEColl();
2357  bool edgespace, rtspace;
2358  bool needint = false;
2359  bool tracespace, rt_tracespace, edge_tracespace;
2360  int dim , p;
2361  PetscBool B_is_Trans = PETSC_FALSE;
2362 
2363  ParMesh *pmesh = (ParMesh *) fespace->GetMesh();
2364  dim = pmesh->Dimension();
2365  // bs = fec->DofForGeometry(Geometry::POINT);
2366  // bs = bs ? bs : 1;
2367  rtspace = dynamic_cast<const RT_FECollection*>(fec);
2368  edgespace = dynamic_cast<const ND_FECollection*>(fec);
2369  edge_tracespace = dynamic_cast<const ND_Trace_FECollection*>(fec);
2370  rt_tracespace = dynamic_cast<const RT_Trace_FECollection*>(fec);
2371  tracespace = edge_tracespace || rt_tracespace;
2372 
2373  p = 1;
2374  if (fespace->GetNE() > 0)
2375  {
2376  if (!tracespace)
2377  {
2378  p = fespace->GetOrder(0);
2379  }
2380  else
2381  {
2382  p = fespace->GetFaceOrder(0);
2383  if (dim == 2) { p++; }
2384  }
2385  }
2386 
2387  if (edgespace) // H(curl)
2388  {
2389  if (dim == 2)
2390  {
2391  needint = true;
2392  if (tracespace)
2393  {
2394  MFEM_WARNING("Tracespace case doesn't work for H(curl) and p=2,"
2395  " not using auxiliary quadrature");
2396  needint = false;
2397  }
2398  }
2399  else
2400  {
2401  FiniteElementCollection *vfec;
2402  if (tracespace)
2403  {
2404  vfec = new H1_Trace_FECollection(p,dim);
2405  }
2406  else
2407  {
2408  vfec = new H1_FECollection(p,dim);
2409  }
2410  ParFiniteElementSpace *vfespace = new ParFiniteElementSpace(pmesh,vfec);
2411  ParDiscreteLinearOperator *grad;
2412  grad = new ParDiscreteLinearOperator(vfespace,fespace);
2413  if (tracespace)
2414  {
2415  grad->AddTraceFaceInterpolator(new GradientInterpolator);
2416  }
2417  else
2418  {
2419  grad->AddDomainInterpolator(new GradientInterpolator);
2420  }
2421  grad->Assemble();
2422  grad->Finalize();
2423  HypreParMatrix *hG = grad->ParallelAssemble();
2424  PetscParMatrix *G = new PetscParMatrix(hG,PETSC_MATAIJ);
2425  delete hG;
2426  delete grad;
2427 
2428  PetscBool conforming = PETSC_TRUE;
2429  if (pmesh->Nonconforming()) { conforming = PETSC_FALSE; }
2430  ierr = PCBDDCSetDiscreteGradient(pc,*G,p,0,PETSC_TRUE,conforming);
2431  PCHKERRQ(pc,ierr);
2432  delete vfec;
2433  delete vfespace;
2434  delete G;
2435  }
2436  }
2437  else if (rtspace) // H(div)
2438  {
2439  needint = true;
2440  if (tracespace)
2441  {
2442  MFEM_WARNING("Tracespace case doesn't work for H(div), not using"
2443  " auxiliary quadrature");
2444  needint = false;
2445  }
2446  }
2447  //else if (bs == dim) // Elasticity?
2448  //{
2449  // needint = true;
2450  //}
2451 
2452  PetscParMatrix *B = NULL;
2453  if (needint)
2454  {
2455  // Generate bilinear form in unassembled format which is used to
2456  // compute the net-flux across subdomain boundaries for H(div) and
2457  // Elasticity, and the line integral \int u x n of 2D H(curl) fields
2458  FiniteElementCollection *auxcoll;
2459  if (tracespace) { auxcoll = new RT_Trace_FECollection(p,dim); }
2460  else { auxcoll = new L2_FECollection(p,dim); };
2461  ParFiniteElementSpace *pspace = new ParFiniteElementSpace(pmesh,auxcoll);
2462  ParMixedBilinearForm *b = new ParMixedBilinearForm(fespace,pspace);
2463 
2464  if (edgespace)
2465  {
2466  if (tracespace)
2467  {
2468  b->AddTraceFaceIntegrator(new VectorFECurlIntegrator);
2469  }
2470  else
2471  {
2472  b->AddDomainIntegrator(new VectorFECurlIntegrator);
2473  }
2474  }
2475  else
2476  {
2477  if (tracespace)
2478  {
2479  b->AddTraceFaceIntegrator(new VectorFEDivergenceIntegrator);
2480  }
2481  else
2482  {
2483  b->AddDomainIntegrator(new VectorFEDivergenceIntegrator);
2484  }
2485  }
2486  b->Assemble();
2487  b->Finalize();
2488  OperatorHandle Bh(Operator::PETSC_MATIS);
2489  b->ParallelAssemble(Bh);
2490  Bh.Get(B);
2491  Bh.SetOperatorOwner(false);
2492 
2493  if (dir) // if essential dofs are present, we need to zero the columns
2494  {
2495  Mat pB = *B;
2496  ierr = MatTranspose(pB,MAT_INPLACE_MATRIX,&pB); PCHKERRQ(pA,ierr);
2497  if (!opts.ess_dof_local)
2498  {
2499  ierr = MatZeroRowsIS(pB,dir,0.,NULL,NULL); PCHKERRQ(pA,ierr);
2500  }
2501  else
2502  {
2503  ierr = MatZeroRowsLocalIS(pB,dir,0.,NULL,NULL); PCHKERRQ(pA,ierr);
2504  }
2505  B_is_Trans = PETSC_TRUE;
2506  }
2507  delete b;
2508  delete pspace;
2509  delete auxcoll;
2510  }
2511 
2512  if (B)
2513  {
2514  ierr = PCBDDCSetDivergenceMat(pc,*B,B_is_Trans,NULL); PCHKERRQ(pc,ierr);
2515  }
2516  delete B;
2517  }
2518  ierr = ISDestroy(&dir); PCHKERRQ(pc,ierr);
2519  ierr = ISDestroy(&neu); PCHKERRQ(pc,ierr);
2520 }
2521 
2522 PetscBDDCSolver::PetscBDDCSolver(PetscParMatrix &A,
2523  const PetscBDDCSolverParams &opts,
2524  const std::string &prefix)
2525  : PetscPreconditioner(A,prefix)
2526 {
2527  BDDCSolverConstructor(opts);
2528  Customize();
2529 }
2530 
2531 PetscBDDCSolver::PetscBDDCSolver(MPI_Comm comm, Operator &op,
2532  const PetscBDDCSolverParams &opts,
2533  const std::string &prefix)
2534  : PetscPreconditioner(comm,op,prefix)
2535 {
2536  BDDCSolverConstructor(opts);
2537  Customize();
2538 }
2539 
2540 PetscFieldSplitSolver::PetscFieldSplitSolver(MPI_Comm comm, Operator &op,
2541  const string &prefix)
2542  : PetscPreconditioner(comm,op,prefix)
2543 {
2544  PC pc = (PC)obj;
2545 
2546  Mat pA;
2547  ierr = PCGetOperators(pc,&pA,NULL); PCHKERRQ(pc,ierr);
2548 
2549  // Check if pA is of type MATNEST
2550  // (this requirement can be removed when we can pass fields).
2551  PetscBool isnest;
2552  ierr = PetscObjectTypeCompare((PetscObject)pA,MATNEST,&isnest);
2553  PCHKERRQ(pA,ierr);
2554  MFEM_VERIFY(isnest,
2555  "PetscFieldSplitSolver needs the matrix in nested format.");
2556 
2557  PetscInt nr;
2558  IS *isrow;
2559  ierr = PCSetType(pc,PCFIELDSPLIT); PCHKERRQ(pc,ierr);
2560  ierr = MatNestGetSize(pA,&nr,NULL); PCHKERRQ(pc,ierr);
2561  ierr = PetscCalloc1(nr,&isrow); CCHKERRQ(PETSC_COMM_SELF,ierr);
2562  ierr = MatNestGetISs(pA,isrow,NULL); PCHKERRQ(pc,ierr);
2563 
2564  // We need to customize here, before setting the index sets.
2565  Customize();
2566 
2567  for (PetscInt i=0; i<nr; i++)
2568  {
2569  ierr = PCFieldSplitSetIS(pc,NULL,isrow[i]); PCHKERRQ(pc,ierr);
2570  }
2571  ierr = PetscFree(isrow); CCHKERRQ(PETSC_COMM_SELF,ierr);
2572 }
2573 
2574 // PetscNonlinearSolver methods
2575 
2576 PetscNonlinearSolver::PetscNonlinearSolver(MPI_Comm comm,
2577  const std::string &prefix)
2578  : PetscSolver(), Solver()
2579 {
2580  // Create the actual solver object
2581  SNES snes;
2582  ierr = SNESCreate(comm, &snes); CCHKERRQ(comm, ierr);
2583  obj = (PetscObject)snes;
2584  ierr = PetscObjectGetClassId(obj, &cid); PCHKERRQ(obj, ierr);
2585  ierr = SNESSetOptionsPrefix(snes, prefix.c_str()); PCHKERRQ(snes, ierr);
2586 
2587  // Allocate private solver context
2588  CreatePrivateContext();
2589 }
2590 
2591 PetscNonlinearSolver::PetscNonlinearSolver(MPI_Comm comm, Operator &op,
2592  const std::string &prefix)
2593  : PetscSolver(), Solver()
2594 {
2595  // Create the actual solver object
2596  SNES snes;
2597  ierr = SNESCreate(comm, &snes); CCHKERRQ(comm, ierr);
2598  obj = (PetscObject)snes;
2599  ierr = PetscObjectGetClassId(obj, &cid); PCHKERRQ(obj, ierr);
2600  ierr = SNESSetOptionsPrefix(snes, prefix.c_str()); PCHKERRQ(snes, ierr);
2601 
2602  // Allocate private solver context
2603  CreatePrivateContext();
2604 
2605  SetOperator(op);
2606 }
2607 
2608 PetscNonlinearSolver::~PetscNonlinearSolver()
2609 {
2610  MPI_Comm comm;
2611  SNES snes = (SNES)obj;
2612  ierr = PetscObjectGetComm(obj,&comm); PCHKERRQ(obj, ierr);
2613  ierr = SNESDestroy(&snes); CCHKERRQ(comm, ierr);
2614 }
2615 
2616 void PetscNonlinearSolver::SetOperator(const Operator &op)
2617 {
2618  SNES snes = (SNES)obj;
2619 
2620  if (operatorset)
2621  {
2622  PetscBool ls,gs;
2623  void *fctx,*jctx;
2624 
2625  ierr = SNESGetFunction(snes, NULL, NULL, &fctx);
2626  PCHKERRQ(snes, ierr);
2627  ierr = SNESGetJacobian(snes, NULL, NULL, NULL, &jctx);
2628  PCHKERRQ(snes, ierr);
2629 
2630  ls = (PetscBool)(height == op.Height() && width == op.Width() &&
2631  (void*)&op == fctx &&
2632  (void*)&op == jctx);
2633  ierr = MPI_Allreduce(&ls,&gs,1,MPIU_BOOL,MPI_LAND,
2634  PetscObjectComm((PetscObject)snes));
2635  PCHKERRQ(snes,ierr);
2636  if (!gs)
2637  {
2638  ierr = SNESReset(snes); PCHKERRQ(snes,ierr);
2639  delete X;
2640  delete B;
2641  X = B = NULL;
2642  }
2643  }
2644 
2645  __mfem_snes_ctx *snes_ctx = (__mfem_snes_ctx*)private_ctx;
2646  snes_ctx->op = (mfem::Operator*)&op;
2647  ierr = SNESSetFunction(snes, NULL, __mfem_snes_function, (void *)snes_ctx);
2648  PCHKERRQ(snes, ierr);
2649  ierr = SNESSetJacobian(snes, NULL, NULL, __mfem_snes_jacobian,
2650  (void *)snes_ctx);
2651  PCHKERRQ(snes, ierr);
2652 
2653  // Update PetscSolver
2654  operatorset = true;
2655 
2656  // Update the Operator fields.
2657  height = op.Height();
2658  width = op.Width();
2659 }
2660 
2661 void PetscNonlinearSolver::SetJacobianType(Operator::Type jacType)
2662 {
2663  __mfem_snes_ctx *snes_ctx = (__mfem_snes_ctx*)private_ctx;
2664  snes_ctx->jacType = jacType;
2665 }
2666 
2667 void PetscNonlinearSolver::Mult(const Vector &b, Vector &x) const
2668 {
2669  SNES snes = (SNES)obj;
2670 
2671  bool b_nonempty = b.Size();
2672  if (!B) { B = new PetscParVector(PetscObjectComm(obj), *this, true); }
2673  if (!X) { X = new PetscParVector(PetscObjectComm(obj), *this, false, false); }
2674  X->PlaceArray(x.GetData());
2675  if (b_nonempty) { B->PlaceArray(b.GetData()); }
2676  else { *B = 0.0; }
2677 
2678  Customize();
2679 
2680  if (!iterative_mode) { *X = 0.; }
2681 
2682  if (bchandler) { bchandler->SetUp(X->Size()); }
2683 
2684  // Solve the system.
2685  ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
2686  X->ResetArray();
2687  if (b_nonempty) { B->ResetArray(); }
2688 }
2689 
2690 // PetscODESolver methods
2691 
2692 PetscODESolver::PetscODESolver(MPI_Comm comm, const string &prefix)
2693  : PetscSolver(), ODESolver()
2694 {
2695  // Create the actual solver object
2696  TS ts;
2697  ierr = TSCreate(comm,&ts); CCHKERRQ(comm,ierr);
2698  obj = (PetscObject)ts;
2699  ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
2700  ierr = TSSetOptionsPrefix(ts, prefix.c_str()); PCHKERRQ(ts, ierr);
2701 
2702  // Allocate private solver context
2703  CreatePrivateContext();
2704 
2705  // Default options, to comply with the current interface to ODESolver.
2706  ierr = TSSetMaxSteps(ts,PETSC_MAX_INT-1);
2707  PCHKERRQ(ts,ierr);
2708  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
2709  PCHKERRQ(ts,ierr);
2710  TSAdapt tsad;
2711  ierr = TSGetAdapt(ts,&tsad);
2712  PCHKERRQ(ts,ierr);
2713  ierr = TSAdaptSetType(tsad,TSADAPTNONE);
2714  PCHKERRQ(ts,ierr);
2715 }
2716 
2717 PetscODESolver::~PetscODESolver()
2718 {
2719  MPI_Comm comm;
2720  TS ts = (TS)obj;
2721  ierr = PetscObjectGetComm(obj,&comm); PCHKERRQ(obj,ierr);
2722  ierr = TSDestroy(&ts); CCHKERRQ(comm,ierr);
2723 }
2724 
2725 void PetscODESolver::Init(TimeDependentOperator &f_,
2726  enum PetscODESolver::Type type)
2727 {
2728  TS ts = (TS)obj;
2729 
2730  __mfem_ts_ctx *ts_ctx = (__mfem_ts_ctx*)private_ctx;
2731  if (operatorset)
2732  {
2733  PetscBool ls,gs;
2734  void *fctx = NULL,*jctx = NULL,*rfctx = NULL,*rjctx = NULL;
2735 
2736  if (f->isImplicit())
2737  {
2738  ierr = TSGetIFunction(ts, NULL, NULL, &fctx);
2739  PCHKERRQ(ts, ierr);
2740  ierr = TSGetIJacobian(ts, NULL, NULL, NULL, &jctx);
2741  PCHKERRQ(ts, ierr);
2742  }
2743  if (!f->isHomogeneous())
2744  {
2745  ierr = TSGetRHSFunction(ts, NULL, NULL, &rfctx);
2746  PCHKERRQ(ts, ierr);
2747  ierr = TSGetRHSJacobian(ts, NULL, NULL, NULL, &rjctx);
2748  PCHKERRQ(ts, ierr);
2749  }
2750  ls = (PetscBool)(f->Height() == f_.Height() &&
2751  f->Width() == f_.Width() &&
2752  f->isExplicit() == f_.isExplicit() &&
2753  f->isImplicit() == f_.isImplicit() &&
2754  f->isHomogeneous() == f_.isHomogeneous());
2755  if (ls && f_.isImplicit())
2756  {
2757  ls = (PetscBool)(ls && (void*)&f_ == fctx && (void*)&f_ == jctx);
2758  }
2759  if (ls && !f_.isHomogeneous())
2760  {
2761  ls = (PetscBool)(ls && (void*)&f_ == rfctx && (void*)&f_ == rjctx);
2762  }
2763  ierr = MPI_Allreduce(&ls,&gs,1,MPIU_BOOL,MPI_LAND,
2764  PetscObjectComm((PetscObject)ts));
2765  PCHKERRQ(ts,ierr);
2766  if (!gs)
2767  {
2768  ierr = TSReset(ts); PCHKERRQ(ts,ierr);
2769  delete X;
2770  X = NULL;
2771  ts_ctx->cached_shift = std::numeric_limits<PetscReal>::min();
2772  ts_ctx->cached_ijacstate = -1;
2773  ts_ctx->cached_rhsjacstate = -1;
2774  }
2775  }
2776  f = &f_;
2777 
2778  // Set functions in TS
2779  ts_ctx->op = &f_;
2780  if (f_.isImplicit())
2781  {
2782  ierr = TSSetIFunction(ts, NULL, __mfem_ts_ifunction, (void *)ts_ctx);
2783  PCHKERRQ(ts, ierr);
2784  ierr = TSSetIJacobian(ts, NULL, NULL, __mfem_ts_ijacobian, (void *)ts_ctx);
2785  PCHKERRQ(ts, ierr);
2786  ierr = TSSetEquationType(ts, TS_EQ_IMPLICIT);
2787  PCHKERRQ(ts, ierr);
2788  }
2789  if (!f_.isHomogeneous())
2790  {
2791  if (!f_.isImplicit())
2792  {
2793  ierr = TSSetEquationType(ts, TS_EQ_EXPLICIT);
2794  PCHKERRQ(ts, ierr);
2795  }
2796  ierr = TSSetRHSFunction(ts, NULL, __mfem_ts_rhsfunction, (void *)ts_ctx);
2797  PCHKERRQ(ts, ierr);
2798  ierr = TSSetRHSJacobian(ts, NULL, NULL, __mfem_ts_rhsjacobian, (void *)ts_ctx);
2799  PCHKERRQ(ts, ierr);
2800  }
2801  operatorset = true;
2802 
2803  ts_ctx->type = type;
2804  if (type == ODE_SOLVER_LINEAR)
2805  {
2806  ierr = TSSetProblemType(ts, TS_LINEAR);
2807  PCHKERRQ(ts, ierr);
2808  }
2809  else
2810  {
2811  ierr = TSSetProblemType(ts, TS_NONLINEAR);
2812  PCHKERRQ(ts, ierr);
2813  }
2814 }
2815 
2816 void PetscODESolver::SetJacobianType(Operator::Type jacType)
2817 {
2818  __mfem_ts_ctx *ts_ctx = (__mfem_ts_ctx*)private_ctx;
2819  ts_ctx->jacType = jacType;
2820 }
2821 
2822 void PetscODESolver::Step(Vector &x, double &t, double &dt)
2823 {
2824  // Pass the parameters to PETSc.
2825  TS ts = (TS)obj;
2826  ierr = TSSetTime(ts, t); PCHKERRQ(ts, ierr);
2827  ierr = TSSetTimeStep(ts, dt); PCHKERRQ(ts, ierr);
2828 
2829  if (!X) { X = new PetscParVector(PetscObjectComm(obj), *f, false, false); }
2830  X->PlaceArray(x.GetData());
2831 
2832  Customize();
2833 
2834  if (bchandler) { bchandler->SetUp(x.Size()); }
2835 
2836  // Take the step.
2837  ierr = TSSetSolution(ts, *X); PCHKERRQ(ts, ierr);
2838  ierr = TSStep(ts); PCHKERRQ(ts, ierr);
2839  X->ResetArray();
2840 
2841  // Get back current time and the time step used to caller.
2842  // We cannot use TSGetTimeStep() as it returns the next candidate step
2843  PetscReal pt;
2844  ierr = TSGetTime(ts,&pt); PCHKERRQ(ts,ierr);
2845  dt = pt - (PetscReal)t;
2846  t = pt;
2847 }
2848 
2849 void PetscODESolver::Run(Vector &x, double &t, double &dt, double t_final)
2850 {
2851  // Give the parameters to PETSc.
2852  TS ts = (TS)obj;
2853  ierr = TSSetTime(ts, t); PCHKERRQ(ts, ierr);
2854  ierr = TSSetTimeStep(ts, dt); PCHKERRQ(ts, ierr);
2855  ierr = TSSetMaxTime(ts, t_final); PCHKERRQ(ts, ierr);
2856  ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
2857  PCHKERRQ(ts, ierr);
2858 
2859  if (!X) { X = new PetscParVector(PetscObjectComm(obj), *f, false, false); }
2860  X->PlaceArray(x.GetData());
2861 
2862  Customize();
2863 
2864  if (bchandler) { bchandler->SetUp(x.Size()); }
2865 
2866  // Take the steps.
2867  ierr = TSSolve(ts, X->x); PCHKERRQ(ts, ierr);
2868  X->ResetArray();
2869 
2870  // Get back final time and time step to caller.
2871  PetscReal pt;
2872  ierr = TSGetTime(ts, &pt); PCHKERRQ(ts,ierr);
2873  t = pt;
2874  ierr = TSGetTimeStep(ts,&pt); PCHKERRQ(ts,ierr);
2875  dt = pt;
2876 }
2877 
2878 } // namespace mfem
2879 
2880 #include "petsc/private/petscimpl.h"
2881 
2882 // auxiliary functions
2883 static PetscErrorCode __mfem_ts_monitor(TS ts, PetscInt it, PetscReal t, Vec x,
2884  void* ctx)
2885 {
2886  mfem::PetscSolverMonitor *monitor_ctx = (mfem::PetscSolverMonitor *)ctx;
2887 
2888  PetscFunctionBeginUser;
2889  if (!ctx)
2890  {
2891  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER, "No monitor context provided");
2892  }
2893  if (monitor_ctx->mon_sol)
2894  {
2895  mfem::PetscParVector V(x,true);
2896  monitor_ctx->MonitorSolution(it,t,V);
2897  }
2898  if (monitor_ctx->mon_res)
2899  {
2900  SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,
2901  "Cannot monitor the residual with TS");
2902  }
2903  PetscFunctionReturn(0);
2904 }
2905 
2906 static PetscErrorCode __mfem_ksp_monitor(KSP ksp, PetscInt it, PetscReal res,
2907  void* ctx)
2908 {
2909  mfem::PetscSolverMonitor *monitor_ctx = (mfem::PetscSolverMonitor *)ctx;
2910  Vec x;
2911  PetscErrorCode ierr;
2912 
2913  PetscFunctionBeginUser;
2914  if (!ctx)
2915  {
2916  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"No monitor context provided");
2917  }
2918  if (monitor_ctx->mon_sol)
2919  {
2920  ierr = KSPBuildSolution(ksp,NULL,&x); CHKERRQ(ierr);
2921  mfem::PetscParVector V(x,true);
2922  monitor_ctx->MonitorSolution(it,res,V);
2923  }
2924  if (monitor_ctx->mon_res)
2925  {
2926  ierr = KSPBuildResidual(ksp,NULL,NULL,&x); CHKERRQ(ierr);
2927  mfem::PetscParVector V(x,true);
2928  monitor_ctx->MonitorResidual(it,res,V);
2929  }
2930  PetscFunctionReturn(0);
2931 }
2932 
2933 static PetscErrorCode __mfem_ts_ifunction(TS ts, PetscReal t, Vec x, Vec xp,
2934  Vec f,void *ctx)
2935 {
2936  __mfem_ts_ctx* ts_ctx = (__mfem_ts_ctx*)ctx;
2937  PetscErrorCode ierr;
2938 
2939  PetscFunctionBeginUser;
2940  mfem::PetscParVector xx(x,true);
2941  mfem::PetscParVector yy(xp,true);
2942  mfem::PetscParVector ff(f,true);
2943 
2944  mfem::TimeDependentOperator *op = ts_ctx->op;
2945  op->SetTime(t);
2946 
2947  if (ts_ctx->bchandler)
2948  {
2949  // we evaluate the ImplicitMult method with the correct bc
2950  if (!ts_ctx->work) { ts_ctx->work = new mfem::Vector(xx.Size()); }
2951  mfem::PetscBCHandler *bchandler = ts_ctx->bchandler;
2952  mfem::Vector* txx = ts_ctx->work;
2953  bchandler->SetTime(t);
2954  bchandler->ApplyBC(xx,*txx);
2955  op->ImplicitMult(*txx,yy,ff);
2956  // and fix the residual (i.e. f_\partial\Omega = u - g(t))
2957  bchandler->FixResidualBC(xx,ff);
2958  }
2959  else
2960  {
2961  // use the ImplicitMult method of the class
2962  op->ImplicitMult(xx,yy,ff);
2963  }
2964 
2965  // need to tell PETSc the Vec has been updated
2966  ierr = PetscObjectStateIncrease((PetscObject)f); CHKERRQ(ierr);
2967  PetscFunctionReturn(0);
2968 }
2969 
2970 static PetscErrorCode __mfem_ts_rhsfunction(TS ts, PetscReal t, Vec x, Vec f,
2971  void *ctx)
2972 {
2973  __mfem_ts_ctx* ts_ctx = (__mfem_ts_ctx*)ctx;
2974  PetscErrorCode ierr;
2975 
2976  PetscFunctionBeginUser;
2977  if (ts_ctx->bchandler) { MFEM_ABORT("RHS evaluation with bc not implemented"); } // TODO
2978  mfem::PetscParVector xx(x,true);
2979  mfem::PetscParVector ff(f,true);
2980  mfem::TimeDependentOperator *top = ts_ctx->op;
2981  top->SetTime(t);
2982 
2983  // use the ExplicitMult method - compute the RHS function
2984  top->ExplicitMult(xx,ff);
2985 
2986  // need to tell PETSc the Vec has been updated
2987  ierr = PetscObjectStateIncrease((PetscObject)f); CHKERRQ(ierr);
2988  PetscFunctionReturn(0);
2989 }
2990 
2991 static PetscErrorCode __mfem_ts_ijacobian(TS ts, PetscReal t, Vec x,
2992  Vec xp, PetscReal shift, Mat A, Mat P,
2993  void *ctx)
2994 {
2995  __mfem_ts_ctx* ts_ctx = (__mfem_ts_ctx*)ctx;
2996  mfem::Vector *xx;
2997  PetscScalar *array;
2998  PetscReal eps = 0.001; /* 0.1% difference */
2999  PetscInt n;
3000  PetscObjectState state;
3001  PetscErrorCode ierr;
3002 
3003  PetscFunctionBeginUser;
3004  // update time
3005  mfem::TimeDependentOperator *op = ts_ctx->op;
3006  op->SetTime(t);
3007 
3008  // prevent to recompute a Jacobian if we already did so
3009  // the relative tolerance comparison should be fine given the fact
3010  // that two consecutive shifts should have similar magnitude
3011  ierr = PetscObjectStateGet((PetscObject)A,&state); CHKERRQ(ierr);
3012  if (ts_ctx->type == mfem::PetscODESolver::ODE_SOLVER_LINEAR &&
3013  std::abs(ts_ctx->cached_shift/shift - 1.0) < eps &&
3014  state == ts_ctx->cached_ijacstate) { PetscFunctionReturn(0); }
3015 
3016  // wrap Vecs with Vectors
3017  ierr = VecGetLocalSize(x,&n); CHKERRQ(ierr);
3018  ierr = VecGetArrayRead(xp,(const PetscScalar**)&array); CHKERRQ(ierr);
3019  mfem::Vector yy(array,n);
3020  ierr = VecRestoreArrayRead(xp,(const PetscScalar**)&array); CHKERRQ(ierr);
3021  ierr = VecGetArrayRead(x,(const PetscScalar**)&array); CHKERRQ(ierr);
3022  if (!ts_ctx->bchandler)
3023  {
3024  xx = new mfem::Vector(array,n);
3025  }
3026  else
3027  {
3028  // make sure we compute a Jacobian with the correct boundary values
3029  if (!ts_ctx->work) { ts_ctx->work = new mfem::Vector(n); }
3030  mfem::Vector txx(array,n);
3031  mfem::PetscBCHandler *bchandler = ts_ctx->bchandler;
3032  xx = ts_ctx->work;
3033  bchandler->SetTime(t);
3034  bchandler->ApplyBC(txx,*xx);
3035  }
3036  ierr = VecRestoreArrayRead(x,(const PetscScalar**)&array); CHKERRQ(ierr);
3037 
3038  // Use TimeDependentOperator::GetImplicitGradient(x,y,s)
3039  mfem::Operator& J = op->GetImplicitGradient(*xx,yy,shift);
3040  if (!ts_ctx->bchandler) { delete xx; }
3041  ts_ctx->cached_shift = shift;
3042 
3043  // Convert to the operator type requested if needed
3044  bool delete_pA = false;
3045  mfem::PetscParMatrix *pA = const_cast<mfem::PetscParMatrix *>
3046  (dynamic_cast<const mfem::PetscParMatrix *>(&J));
3047  if (!pA || pA->GetType() != ts_ctx->jacType)
3048  {
3049  pA = new mfem::PetscParMatrix(PetscObjectComm((PetscObject)ts),&J,
3050  ts_ctx->jacType);
3051  delete_pA = true;
3052  }
3053 
3054  // Eliminate essential dofs
3055  if (ts_ctx->bchandler)
3056  {
3057  mfem::PetscBCHandler *bchandler = ts_ctx->bchandler;
3058  mfem::PetscParVector dummy(PetscObjectComm((PetscObject)ts),0);
3059  pA->EliminateRowsCols(bchandler->GetTDofs(),dummy,dummy);
3060  }
3061 
3062  // Avoid unneeded copy of the matrix by hacking
3063  Mat B;
3064  B = pA->ReleaseMat(false);
3065  ierr = MatHeaderReplace(A,&B); CHKERRQ(ierr);
3066  if (delete_pA) { delete pA; }
3067 
3068  // Jacobian reusage
3069  ierr = PetscObjectStateGet((PetscObject)A,&ts_ctx->cached_ijacstate);
3070  CHKERRQ(ierr);
3071  PetscFunctionReturn(0);
3072 }
3073 
3074 static PetscErrorCode __mfem_ts_rhsjacobian(TS ts, PetscReal t, Vec x,
3075  Mat A, Mat P, void *ctx)
3076 {
3077  __mfem_ts_ctx* ts_ctx = (__mfem_ts_ctx*)ctx;
3078  mfem::Vector *xx;
3079  PetscScalar *array;
3080  PetscInt n;
3081  PetscObjectState state;
3082  PetscErrorCode ierr;
3083 
3084  PetscFunctionBeginUser;
3085  // update time
3086  mfem::TimeDependentOperator *op = ts_ctx->op;
3087  op->SetTime(t);
3088 
3089  // prevent to recompute a Jacobian if we already did so
3090  ierr = PetscObjectStateGet((PetscObject)A,&state); CHKERRQ(ierr);
3091  if (ts_ctx->type == mfem::PetscODESolver::ODE_SOLVER_LINEAR &&
3092  state == ts_ctx->cached_rhsjacstate) { PetscFunctionReturn(0); }
3093 
3094  // wrap Vec with Vector
3095  ierr = VecGetLocalSize(x,&n); CHKERRQ(ierr);
3096  ierr = VecGetArrayRead(x,(const PetscScalar**)&array); CHKERRQ(ierr);
3097  if (!ts_ctx->bchandler)
3098  {
3099  xx = new mfem::Vector(array,n);
3100  }
3101  else
3102  {
3103  // make sure we compute a Jacobian with the correct boundary values
3104  if (!ts_ctx->work) { ts_ctx->work = new mfem::Vector(n); }
3105  mfem::Vector txx(array,n);
3106  mfem::PetscBCHandler *bchandler = ts_ctx->bchandler;
3107  xx = ts_ctx->work;
3108  bchandler->SetTime(t);
3109  bchandler->ApplyBC(txx,*xx);
3110  }
3111  ierr = VecRestoreArrayRead(x,(const PetscScalar**)&array); CHKERRQ(ierr);
3112 
3113  // Use TimeDependentOperator::GetExplicitGradient(x)
3114  mfem::Operator& J = op->GetExplicitGradient(*xx);
3115  if (!ts_ctx->bchandler) { delete xx; }
3116 
3117  // Convert to the operator type requested if needed
3118  bool delete_pA = false;
3119  mfem::PetscParMatrix *pA = const_cast<mfem::PetscParMatrix *>
3120  (dynamic_cast<const mfem::PetscParMatrix *>(&J));
3121  if (!pA || pA->GetType() != ts_ctx->jacType)
3122  {
3123  pA = new mfem::PetscParMatrix(PetscObjectComm((PetscObject)ts),&J,
3124  ts_ctx->jacType);
3125  delete_pA = true;
3126  }
3127 
3128  // Eliminate essential dofs
3129  if (ts_ctx->bchandler)
3130  {
3131  mfem::PetscBCHandler *bchandler = ts_ctx->bchandler;
3132  mfem::PetscParVector dummy(PetscObjectComm((PetscObject)ts),0);
3133  pA->EliminateRowsCols(bchandler->GetTDofs(),dummy,dummy);
3134  }
3135 
3136  // Avoid unneeded copy of the matrix by hacking
3137  Mat B;
3138  B = pA->ReleaseMat(false);
3139  ierr = MatHeaderReplace(A,&B); CHKERRQ(ierr);
3140  if (delete_pA) { delete pA; }
3141 
3142  // Jacobian reusage
3143  if (ts_ctx->type == mfem::PetscODESolver::ODE_SOLVER_LINEAR)
3144  {
3145  ierr = TSRHSJacobianSetReuse(ts,PETSC_TRUE); PCHKERRQ(ts,ierr);
3146  }
3147  ierr = PetscObjectStateGet((PetscObject)A,&ts_ctx->cached_rhsjacstate);
3148  CHKERRQ(ierr);
3149  PetscFunctionReturn(0);
3150 }
3151 
3152 static PetscErrorCode __mfem_snes_monitor(SNES snes, PetscInt it, PetscReal res,
3153  void* ctx)
3154 {
3155  mfem::PetscSolverMonitor *monitor_ctx = (mfem::PetscSolverMonitor *)ctx;
3156  Vec x;
3157  PetscErrorCode ierr;
3158 
3159  PetscFunctionBeginUser;
3160  if (!ctx)
3161  {
3162  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"No monitor context provided");
3163  }
3164  if (monitor_ctx->mon_sol)
3165  {
3166  ierr = SNESGetSolution(snes,&x); CHKERRQ(ierr);
3167  mfem::PetscParVector V(x,true);
3168  monitor_ctx->MonitorSolution(it,res,V);
3169  }
3170  if (monitor_ctx->mon_res)
3171  {
3172  ierr = SNESGetFunction(snes,&x,NULL,NULL); CHKERRQ(ierr);
3173  mfem::PetscParVector V(x,true);
3174  monitor_ctx->MonitorResidual(it,res,V);
3175  }
3176  PetscFunctionReturn(0);
3177 }
3178 
3179 static PetscErrorCode __mfem_snes_jacobian(SNES snes, Vec x, Mat A, Mat P,
3180  void *ctx)
3181 {
3182  PetscScalar *array;
3183  PetscInt n;
3184  PetscErrorCode ierr;
3185  mfem::Vector *xx;
3186  __mfem_snes_ctx *snes_ctx = (__mfem_snes_ctx*)ctx;
3187 
3188  PetscFunctionBeginUser;
3189  ierr = VecGetArrayRead(x,(const PetscScalar**)&array); CHKERRQ(ierr);
3190  ierr = VecGetLocalSize(x,&n); CHKERRQ(ierr);
3191  if (!snes_ctx->bchandler)
3192  {
3193  xx = new mfem::Vector(array,n);
3194  }
3195  else
3196  {
3197  // make sure we compute a Jacobian with the correct boundary values
3198  if (!snes_ctx->work) { snes_ctx->work = new mfem::Vector(n); }
3199  mfem::Vector txx(array,n);
3200  mfem::PetscBCHandler *bchandler = snes_ctx->bchandler;
3201  xx = snes_ctx->work;
3202  bchandler->ApplyBC(txx,*xx);
3203  }
3204 
3205  // Use Operator::GetGradient(x)
3206  mfem::Operator& J = snes_ctx->op->GetGradient(*xx);
3207  ierr = VecRestoreArrayRead(x,(const PetscScalar**)&array); CHKERRQ(ierr);
3208  if (!snes_ctx->bchandler) { delete xx; }
3209 
3210  // Convert to the operator type requested if needed
3211  bool delete_pA = false;
3212  mfem::PetscParMatrix *pA = const_cast<mfem::PetscParMatrix *>
3213  (dynamic_cast<const mfem::PetscParMatrix *>(&J));
3214  if (!pA || pA->GetType() != snes_ctx->jacType)
3215  {
3216  pA = new mfem::PetscParMatrix(PetscObjectComm((PetscObject)snes),&J,
3217  snes_ctx->jacType);
3218  delete_pA = true;
3219  }
3220 
3221  // Eliminate essential dofs
3222  if (snes_ctx->bchandler)
3223  {
3224  mfem::PetscBCHandler *bchandler = snes_ctx->bchandler;
3225  mfem::PetscParVector dummy(PetscObjectComm((PetscObject)snes),0);
3226  pA->EliminateRowsCols(bchandler->GetTDofs(),dummy,dummy);
3227  }
3228 
3229  // Avoid unneeded copy of the matrix by hacking
3230  Mat B = pA->ReleaseMat(false);
3231  ierr = MatHeaderReplace(A,&B); CHKERRQ(ierr);
3232  if (delete_pA) { delete pA; }
3233  PetscFunctionReturn(0);
3234 }
3235 
3236 static PetscErrorCode __mfem_snes_function(SNES snes, Vec x, Vec f, void *ctx)
3237 {
3238  __mfem_snes_ctx* snes_ctx = (__mfem_snes_ctx*)ctx;
3239 
3240  PetscFunctionBeginUser;
3241  mfem::PetscParVector xx(x,true);
3242  mfem::PetscParVector ff(f,true);
3243  if (snes_ctx->bchandler)
3244  {
3245  // we evaluate the Mult method with the correct bc
3246  if (!snes_ctx->work) { snes_ctx->work = new mfem::Vector(xx.Size()); }
3247  mfem::PetscBCHandler *bchandler = snes_ctx->bchandler;
3248  mfem::Vector* txx = snes_ctx->work;
3249  bchandler->ApplyBC(xx,*txx);
3250  snes_ctx->op->Mult(*txx,ff);
3251  // and fix the residual (i.e. f_\partial\Omega = u - g)
3252  bchandler->FixResidualBC(xx,ff);
3253  }
3254  else
3255  {
3256  // use the Mult method of the class
3257  snes_ctx->op->Mult(xx,ff);
3258  }
3259  // need to tell PETSc the Vec has been updated
3260  ierr = PetscObjectStateIncrease((PetscObject)f); CHKERRQ(ierr);
3261  PetscFunctionReturn(0);
3262 }
3263 
3264 static PetscErrorCode __mfem_mat_shell_apply(Mat A, Vec x, Vec y)
3265 {
3266  __mfem_mat_shell_ctx *ctx;
3267  PetscErrorCode ierr;
3268 
3269  PetscFunctionBeginUser;
3270  ierr = MatShellGetContext(A,(void **)&ctx); CHKERRQ(ierr);
3271  mfem::PetscParVector xx(x,true);
3272  mfem::PetscParVector yy(y,true);
3273  ctx->op->Mult(xx,yy);
3274  // need to tell PETSc the Vec has been updated
3275  ierr = PetscObjectStateIncrease((PetscObject)y); CHKERRQ(ierr);
3276  PetscFunctionReturn(0);
3277 }
3278 
3279 static PetscErrorCode __mfem_mat_shell_apply_transpose(Mat A, Vec x, Vec y)
3280 {
3281  __mfem_mat_shell_ctx *ctx;
3282  PetscErrorCode ierr;
3283 
3284  PetscFunctionBeginUser;
3285  ierr = MatShellGetContext(A,(void **)&ctx); CHKERRQ(ierr);
3286  mfem::PetscParVector xx(x,true);
3287  mfem::PetscParVector yy(y,true);
3288  ctx->op->MultTranspose(xx,yy);
3289  // need to tell PETSc the Vec has been updated
3290  ierr = PetscObjectStateIncrease((PetscObject)y); CHKERRQ(ierr);
3291  PetscFunctionReturn(0);
3292 }
3293 
3294 static PetscErrorCode __mfem_mat_shell_destroy(Mat A)
3295 {
3296  __mfem_mat_shell_ctx *ctx;
3297  PetscErrorCode ierr;
3298 
3299  PetscFunctionBeginUser;
3300  ierr = MatShellGetContext(A,(void **)&ctx); CHKERRQ(ierr);
3301  delete ctx;
3302  PetscFunctionReturn(0);
3303 }
3304 
3305 static PetscErrorCode __mfem_pc_shell_view(PC pc, PetscViewer viewer)
3306 {
3307  __mfem_pc_shell_ctx *ctx;
3308  PetscErrorCode ierr;
3309 
3310  PetscFunctionBeginUser;
3311  ierr = PCShellGetContext(pc,(void **)&ctx); CHKERRQ(ierr);
3312  if (ctx->op)
3313  {
3314  PetscBool isascii;
3315  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
3316  CHKERRQ(ierr);
3317 
3318  mfem::PetscPreconditioner *ppc = dynamic_cast<mfem::PetscPreconditioner *>
3319  (ctx->op);
3320  if (ppc)
3321  {
3322  ierr = PCView(*ppc,viewer); CHKERRQ(ierr);
3323  }
3324  else
3325  {
3326  if (isascii)
3327  {
3328  ierr = PetscViewerASCIIPrintf(viewer,
3329  "No information available on the mfem::Solver\n");
3330  CHKERRQ(ierr);
3331  }
3332  }
3333  if (isascii && ctx->factory)
3334  {
3335  ierr = PetscViewerASCIIPrintf(viewer,
3336  "Number of preconditioners created by the factory %lu\n",ctx->numprec);
3337  CHKERRQ(ierr);
3338  }
3339  }
3340  PetscFunctionReturn(0);
3341 }
3342 
3343 static PetscErrorCode __mfem_pc_shell_apply(PC pc, Vec x, Vec y)
3344 {
3345  __mfem_pc_shell_ctx *ctx;
3346  PetscErrorCode ierr;
3347 
3348  PetscFunctionBeginUser;
3349  mfem::PetscParVector xx(x,true);
3350  mfem::PetscParVector yy(y,true);
3351  ierr = PCShellGetContext(pc,(void **)&ctx); CHKERRQ(ierr);
3352  if (ctx->op)
3353  {
3354  ctx->op->Mult(xx,yy);
3355  // need to tell PETSc the Vec has been updated
3356  ierr = PetscObjectStateIncrease((PetscObject)y); CHKERRQ(ierr);
3357  }
3358  else // operator is not present, copy x
3359  {
3360  yy = xx;
3361  }
3362  PetscFunctionReturn(0);
3363 }
3364 
3365 static PetscErrorCode __mfem_pc_shell_apply_transpose(PC pc, Vec x, Vec y)
3366 {
3367  __mfem_pc_shell_ctx *ctx;
3368  PetscErrorCode ierr;
3369 
3370  PetscFunctionBeginUser;
3371  mfem::PetscParVector xx(x,true);
3372  mfem::PetscParVector yy(y,true);
3373  ierr = PCShellGetContext(pc,(void **)&ctx); CHKERRQ(ierr);
3374  if (ctx->op)
3375  {
3376  ctx->op->MultTranspose(xx,yy);
3377  // need to tell PETSc the Vec has been updated
3378  ierr = PetscObjectStateIncrease((PetscObject)y); CHKERRQ(ierr);
3379  }
3380  else // operator is not present, copy x
3381  {
3382  yy = xx;
3383  }
3384  PetscFunctionReturn(0);
3385 }
3386 
3387 static PetscErrorCode __mfem_pc_shell_setup(PC pc)
3388 {
3389  __mfem_pc_shell_ctx *ctx;
3390 
3391  PetscFunctionBeginUser;
3392  ierr = PCShellGetContext(pc,(void **)&ctx); CHKERRQ(ierr);
3393  if (ctx->factory)
3394  {
3395  // Delete any owned operator
3396  if (ctx->ownsop)
3397  {
3398  delete ctx->op;
3399  }
3400 
3401  // Get current preconditioning Mat
3402  Mat B;
3403  ierr = PCGetOperators(pc,NULL,&B); CHKERRQ(ierr);
3404 
3405  // Call user-defined setup
3406  mfem::OperatorHandle hB(new mfem::PetscParMatrix(B,true),true);
3407  mfem::PetscPreconditionerFactory *factory = ctx->factory;
3408  ctx->op = factory->NewPreconditioner(hB);
3409  ctx->ownsop = true;
3410  ctx->numprec++;
3411  }
3412  PetscFunctionReturn(0);
3413 }
3414 
3415 static PetscErrorCode __mfem_pc_shell_destroy(PC pc)
3416 {
3417  __mfem_pc_shell_ctx *ctx;
3418  PetscErrorCode ierr;
3419 
3420  PetscFunctionBeginUser;
3421  ierr = PCShellGetContext(pc,(void **)&ctx); CHKERRQ(ierr);
3422  if (ctx->ownsop)
3423  {
3424  delete ctx->op;
3425  }
3426  delete ctx;
3427  PetscFunctionReturn(0);
3428 }
3429 
3430 static PetscErrorCode __mfem_array_container_destroy(void *ptr)
3431 {
3432  PetscErrorCode ierr;
3433 
3434  PetscFunctionBeginUser;
3435  ierr = PetscFree(ptr); CHKERRQ(ierr);
3436  PetscFunctionReturn(0);
3437 }
3438 
3439 static PetscErrorCode __mfem_matarray_container_destroy(void *ptr)
3440 {
3441  mfem::Array<Mat> *a = (mfem::Array<Mat>*)ptr;
3442  PetscErrorCode ierr;
3443 
3444  PetscFunctionBeginUser;
3445  for (int i=0; i<a->Size(); i++)
3446  {
3447  Mat M = (*a)[i];
3448  MPI_Comm comm = PetscObjectComm((PetscObject)M);
3449  ierr = MatDestroy(&M); CCHKERRQ(comm,ierr);
3450  }
3451  delete a;
3452  PetscFunctionReturn(0);
3453 }
3454 
3455 // Sets the type of PC to PCSHELL and wraps the solver action
3456 // if ownsop is true, ownership of precond is transferred to the PETSc object
3457 PetscErrorCode MakeShellPC(PC pc, mfem::Solver &precond, bool ownsop)
3458 {
3459  PetscFunctionBeginUser;
3460  __mfem_pc_shell_ctx *ctx = new __mfem_pc_shell_ctx;
3461  ctx->op = &precond;
3462  ctx->ownsop = ownsop;
3463  ctx->factory = NULL;
3464  ctx->numprec = 0;
3465 
3466  ierr = PCSetType(pc,PCSHELL); CHKERRQ(ierr);
3467  ierr = PCShellSetName(pc,"MFEM Solver"); CHKERRQ(ierr);
3468  ierr = PCShellSetContext(pc,(void *)ctx); CHKERRQ(ierr);
3469  ierr = PCShellSetApply(pc,__mfem_pc_shell_apply); CHKERRQ(ierr);
3470  ierr = PCShellSetApplyTranspose(pc,__mfem_pc_shell_apply_transpose);
3471  CHKERRQ(ierr);
3472  ierr = PCShellSetSetUp(pc,__mfem_pc_shell_setup); CHKERRQ(ierr);
3473  ierr = PCShellSetView(pc,__mfem_pc_shell_view); CHKERRQ(ierr);
3474  ierr = PCShellSetDestroy(pc,__mfem_pc_shell_destroy); CHKERRQ(ierr);
3475  PetscFunctionReturn(0);
3476 }
3477 
3478 // Sets the type of PC to PCSHELL. Uses a PetscPreconditionerFactory to construct the solver
3479 // Takes ownership of the solver created by the factory
3480 PetscErrorCode MakeShellPCWithFactory(PC pc,
3481  mfem::PetscPreconditionerFactory *factory)
3482 {
3483  PetscFunctionBeginUser;
3484  __mfem_pc_shell_ctx *ctx = new __mfem_pc_shell_ctx;
3485  ctx->op = NULL;
3486  ctx->ownsop = true;
3487  ctx->factory = factory;
3488  ctx->numprec = 0;
3489 
3490  ierr = PCSetType(pc,PCSHELL); CHKERRQ(ierr);
3491  ierr = PCShellSetName(pc,factory->GetName()); CHKERRQ(ierr);
3492  ierr = PCShellSetContext(pc,(void *)ctx); CHKERRQ(ierr);
3493  ierr = PCShellSetApply(pc,__mfem_pc_shell_apply); CHKERRQ(ierr);
3494  ierr = PCShellSetApplyTranspose(pc,__mfem_pc_shell_apply_transpose);
3495  CHKERRQ(ierr);
3496  ierr = PCShellSetSetUp(pc,__mfem_pc_shell_setup); CHKERRQ(ierr);
3497  ierr = PCShellSetView(pc,__mfem_pc_shell_view); CHKERRQ(ierr);
3498  ierr = PCShellSetDestroy(pc,__mfem_pc_shell_destroy); CHKERRQ(ierr);
3499  PetscFunctionReturn(0);
3500 }
3501 
3502 // Converts from a list (or a marked Array if islist is false) to an IS
3503 // st indicates the offset where to start numbering
3504 static PetscErrorCode Convert_Array_IS(MPI_Comm comm, bool islist,
3505  const mfem::Array<int> *list,
3506  PetscInt st, IS* is)
3507 {
3508  PetscInt n = list->Size(),*idxs;
3509  const int *data = list->GetData();
3510  PetscErrorCode ierr;
3511 
3512  PetscFunctionBeginUser;
3513  ierr = PetscMalloc1(n,&idxs); CHKERRQ(ierr);
3514  if (islist)
3515  {
3516  for (PetscInt i=0; i<n; i++) { idxs[i] = data[i] + st; }
3517  }
3518  else
3519  {
3520  PetscInt cum = 0;
3521  for (PetscInt i=0; i<n; i++)
3522  {
3523  if (data[i]) { idxs[cum++] = i+st; }
3524  }
3525  n = cum;
3526  }
3527  ierr = ISCreateGeneral(comm,n,idxs,PETSC_OWN_POINTER,is);
3528  CHKERRQ(ierr);
3529  PetscFunctionReturn(0);
3530 }
3531 
3532 // Converts from a marked Array of Vdofs to an IS
3533 // st indicates the offset where to start numbering
3534 // l2l is a vector of matrices generated during RAP
3535 static PetscErrorCode Convert_Vmarks_IS(MPI_Comm comm,
3536  mfem::Array<Mat> &pl2l,
3537  const mfem::Array<int> *mark,
3538  PetscInt st, IS* is)
3539 {
3540  mfem::Array<int> sub_dof_marker;
3541  mfem::Array<mfem::SparseMatrix*> l2l(pl2l.Size());
3542  PetscInt nl;
3543  PetscErrorCode ierr;
3544 
3545  PetscFunctionBeginUser;
3546  for (int i = 0; i < pl2l.Size(); i++)
3547  {
3548  PetscInt m,n,*ii,*jj;
3549  PetscBool done;
3550  ierr = MatGetRowIJ(pl2l[i],0,PETSC_FALSE,PETSC_FALSE,&m,(const PetscInt**)&ii,
3551  (const PetscInt**)&jj,&done); CHKERRQ(ierr);
3552  MFEM_VERIFY(done,"Unable to perform MatGetRowIJ on " << i << " l2l matrix");
3553  ierr = MatGetSize(pl2l[i],NULL,&n); CHKERRQ(ierr);
3554  l2l[i] = new mfem::SparseMatrix(ii,jj,NULL,m,n,false,true,true);
3555  }
3556  nl = 0;
3557  for (int i = 0; i < l2l.Size(); i++) { nl += l2l[i]->Width(); }
3558  sub_dof_marker.SetSize(nl);
3559  const int* vdata = mark->GetData();
3560  int* sdata = sub_dof_marker.GetData();
3561  int cumh = 0, cumw = 0;
3562  for (int i = 0; i < l2l.Size(); i++)
3563  {
3564  const mfem::Array<int> vf_marker(const_cast<int*>(vdata)+cumh,
3565  l2l[i]->Height());
3566  mfem::Array<int> sf_marker(sdata+cumw,l2l[i]->Width());
3567  l2l[i]->BooleanMultTranspose(vf_marker,sf_marker);
3568  cumh += l2l[i]->Height();
3569  cumw += l2l[i]->Width();
3570  }
3571  ierr = Convert_Array_IS(comm,false,&sub_dof_marker,st,is); CCHKERRQ(comm,ierr);
3572  for (int i = 0; i < pl2l.Size(); i++)
3573  {
3574  PetscInt m = l2l[i]->Height();
3575  PetscInt *ii = l2l[i]->GetI(),*jj = l2l[i]->GetJ();
3576  PetscBool done;
3577  ierr = MatRestoreRowIJ(pl2l[i],0,PETSC_FALSE,PETSC_FALSE,&m,
3578  (const PetscInt**)&ii,
3579  (const PetscInt**)&jj,&done); CHKERRQ(ierr);
3580  MFEM_VERIFY(done,"Unable to perform MatRestoreRowIJ on "
3581  << i << " l2l matrix");
3582  delete l2l[i];
3583  }
3584  PetscFunctionReturn(0);
3585 }
3586 
3587 #if !defined(PETSC_HAVE_HYPRE)
3588 
3589 #include "_hypre_parcsr_mv.h"
3590 static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix* hA,Mat* pA)
3591 {
3592  MPI_Comm comm = hypre_ParCSRMatrixComm(hA);
3593  hypre_CSRMatrix *hdiag,*hoffd;
3594  PetscScalar *da,*oa,*aptr;
3595  PetscInt *dii,*djj,*oii,*ojj,*iptr;
3596  PetscInt i,dnnz,onnz,m,n;
3597  PetscMPIInt size;
3598  PetscErrorCode ierr;
3599 
3600  PetscFunctionBeginUser;
3601  hdiag = hypre_ParCSRMatrixDiag(hA);
3602  hoffd = hypre_ParCSRMatrixOffd(hA);
3603  m = hypre_CSRMatrixNumRows(hdiag);
3604  n = hypre_CSRMatrixNumCols(hdiag);
3605  dnnz = hypre_CSRMatrixNumNonzeros(hdiag);
3606  onnz = hypre_CSRMatrixNumNonzeros(hoffd);
3607  ierr = PetscMalloc1(m+1,&dii); CHKERRQ(ierr);
3608  ierr = PetscMalloc1(dnnz,&djj); CHKERRQ(ierr);
3609  ierr = PetscMalloc1(dnnz,&da); CHKERRQ(ierr);
3610  ierr = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));
3611  CHKERRQ(ierr);
3612  ierr = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));
3613  CHKERRQ(ierr);
3614  ierr = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));
3615  CHKERRQ(ierr);
3616  iptr = djj;
3617  aptr = da;
3618  for (i=0; i<m; i++)
3619  {
3620  PetscInt nc = dii[i+1]-dii[i];
3621  ierr = PetscSortIntWithScalarArray(nc,iptr,aptr); CHKERRQ(ierr);
3622  iptr += nc;
3623  aptr += nc;
3624  }
3625  ierr = MPI_Comm_size(comm,&size); CHKERRQ(ierr);
3626  if (size > 1)
3627  {
3628  PetscInt *offdj,*coffd;
3629 
3630  ierr = PetscMalloc1(m+1,&oii); CHKERRQ(ierr);
3631  ierr = PetscMalloc1(onnz,&ojj); CHKERRQ(ierr);
3632  ierr = PetscMalloc1(onnz,&oa); CHKERRQ(ierr);
3633  ierr = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));
3634  CHKERRQ(ierr);
3635  offdj = hypre_CSRMatrixJ(hoffd);
3636  coffd = hypre_ParCSRMatrixColMapOffd(hA);
3637  for (i=0; i<onnz; i++) { ojj[i] = coffd[offdj[i]]; }
3638  ierr = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));
3639  CHKERRQ(ierr);
3640  iptr = ojj;
3641  aptr = oa;
3642  for (i=0; i<m; i++)
3643  {
3644  PetscInt nc = oii[i+1]-oii[i];
3645  ierr = PetscSortIntWithScalarArray(nc,iptr,aptr); CHKERRQ(ierr);
3646  iptr += nc;
3647  aptr += nc;
3648  }
3649  ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,
3650  djj,da,oii,ojj,oa,pA); CHKERRQ(ierr);
3651  }
3652  else
3653  {
3654  oii = ojj = NULL;
3655  oa = NULL;
3656  ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,pA); CHKERRQ(ierr);
3657  }
3658  /* We are responsible to free the CSR arrays. However, since we can take
3659  references of a PetscParMatrix but we cannot take reference of PETSc
3660  arrays, we need to create a PetscContainer object to take reference of
3661  these arrays in reference objects */
3662  void *ptrs[6] = {dii,djj,da,oii,ojj,oa};
3663  const char *names[6] = {"_mfem_csr_dii",
3664  "_mfem_csr_djj",
3665  "_mfem_csr_da",
3666  "_mfem_csr_oii",
3667  "_mfem_csr_ojj",
3668  "_mfem_csr_oa"
3669  };
3670  for (i=0; i<6; i++)
3671  {
3672  PetscContainer c;
3673 
3674  ierr = PetscContainerCreate(comm,&c); CHKERRQ(ierr);
3675  ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
3676  ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
3677  CHKERRQ(ierr);
3678  ierr = PetscObjectCompose((PetscObject)(*pA),names[i],(PetscObject)c);
3679  CHKERRQ(ierr);
3680  ierr = PetscContainerDestroy(&c); CHKERRQ(ierr);
3681  }
3682  PetscFunctionReturn(0);
3683 }
3684 
3685 static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix* hA,Mat* pA)
3686 {
3687  Mat lA;
3688  ISLocalToGlobalMapping rl2g,cl2g;
3689  IS is;
3690  hypre_CSRMatrix *hdiag,*hoffd;
3691  MPI_Comm comm = hypre_ParCSRMatrixComm(hA);
3692  void *ptrs[2];
3693  const char *names[2] = {"_mfem_csr_aux",
3694  "_mfem_csr_data"
3695  };
3696  PetscScalar *hdd,*hod,*aa,*data;
3697  PetscInt *col_map_offd,*hdi,*hdj,*hoi,*hoj;
3698  PetscInt *aux,*ii,*jj;
3699  PetscInt cum,dr,dc,oc,str,stc,nnz,i,jd,jo;
3700  PetscErrorCode ierr;
3701 
3702  PetscFunctionBeginUser;
3703  /* access relevant information in ParCSR */
3704  str = hypre_ParCSRMatrixFirstRowIndex(hA);
3705  stc = hypre_ParCSRMatrixFirstColDiag(hA);
3706  hdiag = hypre_ParCSRMatrixDiag(hA);
3707  hoffd = hypre_ParCSRMatrixOffd(hA);
3708  dr = hypre_CSRMatrixNumRows(hdiag);
3709  dc = hypre_CSRMatrixNumCols(hdiag);
3710  nnz = hypre_CSRMatrixNumNonzeros(hdiag);
3711  hdi = hypre_CSRMatrixI(hdiag);
3712  hdj = hypre_CSRMatrixJ(hdiag);
3713  hdd = hypre_CSRMatrixData(hdiag);
3714  oc = hypre_CSRMatrixNumCols(hoffd);
3715  nnz += hypre_CSRMatrixNumNonzeros(hoffd);
3716  hoi = hypre_CSRMatrixI(hoffd);
3717  hoj = hypre_CSRMatrixJ(hoffd);
3718  hod = hypre_CSRMatrixData(hoffd);
3719 
3720  /* generate l2g maps for rows and cols */
3721  ierr = ISCreateStride(comm,dr,str,1,&is); CHKERRQ(ierr);
3722  ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g); CHKERRQ(ierr);
3723  ierr = ISDestroy(&is); CHKERRQ(ierr);
3724  col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
3725  ierr = PetscMalloc1(dc+oc,&aux); CHKERRQ(ierr);
3726  for (i=0; i<dc; i++) { aux[i] = i+stc; }
3727  for (i=0; i<oc; i++) { aux[i+dc] = col_map_offd[i]; }
3728  ierr = ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is); CHKERRQ(ierr);
3729  ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g); CHKERRQ(ierr);
3730  ierr = ISDestroy(&is); CHKERRQ(ierr);
3731 
3732  /* create MATIS object */
3733  ierr = MatCreate(comm,pA); CHKERRQ(ierr);
3734  ierr = MatSetSizes(*pA,dr,dc,PETSC_DECIDE,PETSC_DECIDE); CHKERRQ(ierr);
3735  ierr = MatSetType(*pA,MATIS); CHKERRQ(ierr);
3736  ierr = MatSetLocalToGlobalMapping(*pA,rl2g,cl2g); CHKERRQ(ierr);
3737  ierr = ISLocalToGlobalMappingDestroy(&rl2g); CHKERRQ(ierr);
3738  ierr = ISLocalToGlobalMappingDestroy(&cl2g); CHKERRQ(ierr);
3739 
3740  /* merge local matrices */
3741  ierr = PetscMalloc1(nnz+dr+1,&aux); CHKERRQ(ierr);
3742  ierr = PetscMalloc1(nnz,&data); CHKERRQ(ierr);
3743  ii = aux;
3744  jj = aux+dr+1;
3745  aa = data;
3746  *ii = *(hdi++) + *(hoi++);
3747  for (jd=0,jo=0,cum=0; *ii<nnz; cum++)
3748  {
3749  PetscScalar *aold = aa;
3750  PetscInt *jold = jj,nc = jd+jo;
3751  for (; jd<*hdi; jd++) { *jj++ = *hdj++; *aa++ = *hdd++; }
3752  for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
3753  *(++ii) = *(hdi++) + *(hoi++);
3754  ierr = PetscSortIntWithScalarArray(jd+jo-nc,jold,aold); CHKERRQ(ierr);
3755  }
3756  for (; cum<dr; cum++) { *(++ii) = nnz; }
3757  ii = aux;
3758  jj = aux+dr+1;
3759  aa = data;
3760  ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,ii,jj,aa,&lA);
3761  CHKERRQ(ierr);
3762  ptrs[0] = aux;
3763  ptrs[1] = data;
3764  for (i=0; i<2; i++)
3765  {
3766  PetscContainer c;
3767 
3768  ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); CHKERRQ(ierr);
3769  ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
3770  ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
3771  CHKERRQ(ierr);
3772  ierr = PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
3773  CHKERRQ(ierr);
3774  ierr = PetscContainerDestroy(&c); CHKERRQ(ierr);
3775  }
3776  ierr = MatISSetLocalMat(*pA,lA); CHKERRQ(ierr);
3777  ierr = MatDestroy(&lA); CHKERRQ(ierr);
3778  ierr = MatAssemblyBegin(*pA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
3779  ierr = MatAssemblyEnd(*pA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
3780  PetscFunctionReturn(0);
3781 }
3782 #endif
3783 
3784 #endif // MFEM_USE_PETSC
3785 #endif // MFEM_USE_MPI
mfem::Array::Size
int Size() const
Logical size of the array.
Definition: array.hpp:110
mfem::PetscLinearSolver::SetPreconditioner
void SetPreconditioner(Solver &precond)
Sets the solver to perform preconditioning.
Definition: petsc.cpp:2003
mfem::PetscSolver::SetPrintLevel
void SetPrintLevel(int plev)
Definition: petsc.cpp:1407
mfem::SparseMatrix::NumNonZeroElems
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:972
mfem::PetscSolver::B
PetscParVector * B
Right-hand side and solution vector.
Definition: petsc.hpp:416
mfem::PetscSolver::CreatePrivateContext
void CreatePrivateContext()
Definition: petsc.cpp:1679
mfem::PetscParMatrix::Print
void Print(const char *fname=NULL, bool binary=false) const
Prints the matrix (to stdout if fname is NULL)
Definition: petsc.cpp:1072
mfem::Vector::Vector
Vector()
Default constructor for Vector. Sets size = 0 and data = NULL.
Definition: vector.hpp:51
mfem::HypreParMatrix::GetComm
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:298
mfem::TimeDependentOperator::isImplicit
bool isImplicit() const
True if type is IMPLICIT or HOMOGENEOUS.
Definition: operator.hpp:177
mfem::TimeDependentOperator::isHomogeneous
bool isHomogeneous() const
True if type is HOMOGENEOUS.
Definition: operator.hpp:179
mfem::PetscNonlinearSolver::SetOperator
virtual void SetOperator(const Operator &op)
Specification of the nonlinear operator.
Definition: petsc.cpp:2616
mfem::PetscPreconditioner
Abstract class for PETSc&#39;s preconditioners.
Definition: petsc.hpp:521
mfem::PetscParMatrix::operator+=
PetscParMatrix & operator+=(const PetscParMatrix &B)
Definition: petsc.cpp:484
mfem::PetscLinearSolver::SetOperator
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: petsc.cpp:1839
mfem::RAP
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1476
mfem::PetscSolver::clcustom
bool clcustom
Boolean to handle SetFromOptions calls.
Definition: petsc.hpp:407
mfem::PetscParMatrix::ConvertOperator
void ConvertOperator(MPI_Comm comm, const Operator &op, Mat *B, Operator::Type tid)
Definition: petsc.cpp:656
mfem::PetscSolver::GetFinalNorm
double GetFinalNorm()
Definition: petsc.cpp:1654
mfem::PetscParMatrix
Wrapper for PETSc&#39;s matrix class.
Definition: petsc.hpp:133
mfem::PetscBDDCSolverParams::nat_dof
const Array< int > * nat_dof
Definition: petsc.hpp:549
mfem::PetscBCHandler::ApplyBC
void ApplyBC(const Vector &x, Vector &y)
y = x on ess_tdof_list_c and y = g (internally evaluated) on ess_tdof_list
Definition: petsc.cpp:1758
mfem::PetscLinearSolver
Abstract class for PETSc&#39;s linear solvers.
Definition: petsc.hpp:474
mfem::PetscParVector::x
Vec x
The actual PETSc object.
Definition: petsc.hpp:37
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
mfem::PetscSolver
Abstract class for PETSc&#39;s solvers.
Definition: petsc.hpp:403
mfem::PetscODESolver::Run
virtual void Run(Vector &x, double &t, double &dt, double t_final)
Perform time integration from time t [in] to time tf [in].
Definition: petsc.cpp:2849
mfem::PetscBDDCSolverParams::fespace
ParFiniteElementSpace * fespace
Definition: petsc.hpp:546
mfem::TimeDependentOperator
Base abstract class for time dependent operators.
Definition: operator.hpp:142
mfem::PetscParMatrix::GetNumCols
PetscInt GetNumCols() const
Returns the local number of columns.
Definition: petsc.cpp:366
mfem::Operator::GetGradient
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
Definition: operator.hpp:57
mfem::Vector::MakeDataOwner
void MakeDataOwner()
Definition: vector.hpp:107
mfem::Operator::Width
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
mfem::OperatorHandle
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
mfem::PetscParMatrix::Mult
void Mult(double a, const Vector &x, double b, Vector &y) const
Matvec: y = a A x + b y.
Definition: petsc.cpp:1039
mfem::PetscBDDCSolver::PetscBDDCSolver
PetscBDDCSolver(MPI_Comm comm, Operator &op, const PetscBDDCSolverParams &opts=PetscBDDCSolverParams(), const std::string &prefix=std::string())
Definition: petsc.cpp:2531
mfem::PetscSolver::~PetscSolver
virtual ~PetscSolver()
Destroy the PetscParVectors allocated (if any).
Definition: petsc.cpp:1318
mfem::PetscNonlinearSolver::PetscNonlinearSolver
PetscNonlinearSolver(MPI_Comm comm, const std::string &prefix=std::string())
Definition: petsc.cpp:2576
mfem::Array::GetData
T * GetData()
Returns the data.
Definition: array.hpp:92
mfem::PetscParMatrix::GetType
Type GetType() const
Definition: petsc.cpp:1260
mfem::TimeDependentOperator::SetTime
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:172
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
mfem::PetscODESolver::Step
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: petsc.cpp:2822
mfem::ODESolver
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
mfem::PetscParMatrix::Init
void Init()
Initialize with defaults. Does not initialize inherited members.
Definition: petsc.cpp:394
mfem::PetscPreconditioner::PetscPreconditioner
PetscPreconditioner(MPI_Comm comm, const std::string &prefix=std::string())
Definition: petsc.cpp:2116
mfem::PetscSolver::cid
PetscClassId cid
The class id of the actual PETSc object.
Definition: petsc.hpp:413
mfem::BlockOperator::IsZeroBlock
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
Definition: blockoperator.hpp:75
mfem::ParFiniteElementSpace
Abstract parallel finite element space.
Definition: pfespace.hpp:31
mfem::PetscBCHandler::CONSTANT
Constant in time b.c.
Definition: petsc.hpp:334
mfem::PetscSolver::SetRelTol
void SetRelTol(double tol)
Definition: petsc.cpp:1330
mfem::Solver::iterative_mode
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:263
mfem::PetscPCGSolver::PetscPCGSolver
PetscPCGSolver(MPI_Comm comm, const std::string &prefix=std::string())
Definition: petsc.cpp:2086
mfem::PetscParMatrix::Destroy
void Destroy()
Delete all owned data. Does not perform re-initialization with defaults.
Definition: petsc.cpp:907
mfem::TimeDependentOperator::GetExplicitGradient
virtual Operator & GetExplicitGradient(const Vector &x) const
Return an Operator representing dG/dx at the given point x and the currently set time.
Definition: operator.hpp:248
mfem::PetscSolver::SetMonitor
void SetMonitor(PetscSolverMonitor *ctx)
Sets user-defined monitoring routine.
Definition: petsc.cpp:1479
mfem::PetscSolver::SetMaxIter
void SetMaxIter(int max_iter)
Definition: petsc.cpp:1380
mfem::Vector::GetData
double * GetData() const
Definition: vector.hpp:121
mfem::SparseMatrix::SortColumnIndices
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:326
mfem::PetscParVector
Wrapper for PETSc&#39;s vector class.
Definition: petsc.hpp:33
mfem::PetscSolver::Customize
void Customize(bool customize=true) const
Customize object with options set.
Definition: petsc.cpp:1555
mfem::PetscParMatrix::M
PetscInt M() const
Returns the global number of rows.
Definition: petsc.cpp:373
mfem::BlockOperator::GetBlock
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
Definition: blockoperator.hpp:77
mfem::PetscNonlinearSolver::Mult
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
Definition: petsc.cpp:2667
mfem::PetscParMatrix::Y
PetscParVector * Y
Definition: petsc.hpp:140
mfem::PetscLinearSolver::~PetscLinearSolver
virtual ~PetscLinearSolver()
Definition: petsc.cpp:2076
mfem::PetscParMatrix::GetNumRows
PetscInt GetNumRows() const
Returns the local number of rows.
Definition: petsc.cpp:359
mfem::PetscParVector::Randomize
void Randomize(PetscInt seed)
Set random values.
Definition: petsc.cpp:320
mfem::PetscParVector::Print
void Print(const char *fname=NULL, bool binary=false) const
Prints the vector (to stdout if fname is NULL)
Definition: petsc.cpp:332
mfem::EliminateBC
void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Definition: hypre.cpp:1530
mfem::PetscParVector::ResetArray
void ResetArray()
Reset the PETSc Vec object to use its default data. Call this method after the use of PlaceArray()...
Definition: petsc.cpp:315
mfem::TimeDependentOperator::GetImplicitGradient
virtual Operator & GetImplicitGradient(const Vector &x, const Vector &k, double shift) const
Return an Operator representing (dF/dk shift + dF/dx) at the given x, k, and the currently set time...
Definition: operator.hpp:235
mfem::PetscODESolver::Init
virtual void Init(TimeDependentOperator &f_, enum PetscODESolver::Type type)
Initialize the ODE solver.
Definition: petsc.cpp:2725
mfem::PetscParMatrix::Transpose
PetscParMatrix * Transpose(bool action=false)
Returns the transpose of the PetscParMatrix.
Definition: petsc.cpp:1020
dim
int dim
Definition: ex3.cpp:47
mfem::ParFiniteElementSpace::GetComm
MPI_Comm GetComm() const
Definition: pfespace.hpp:166
mfem::PetscSolver::obj
PetscObject obj
The actual PETSc object (KSP, PC, SNES or TS).
Definition: petsc.hpp:410
mfem::PetscSolverMonitor
Abstract class for monitoring PETSc&#39;s solvers.
Definition: petsc.hpp:659
mfem::SparseMatrix
Data type sparse matrix.
Definition: sparsemat.hpp:38
mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
mfem::PetscNonlinearSolver::SetJacobianType
void SetJacobianType(Operator::Type type)
Definition: petsc.cpp:2661
mfem::to_int
int to_int(const std::string &str)
Definition: text.hpp:68
mfem::PetscBDDCSolverParams
Auxiliary class for BDDC customization.
Definition: petsc.hpp:543
mfem::PetscPreconditioner::~PetscPreconditioner
virtual ~PetscPreconditioner()
Definition: petsc.cpp:2227
mfem::PetscParMatrix::GetComm
MPI_Comm GetComm() const
Get the associated MPI communicator.
Definition: petsc.hpp:242
mfem::PetscBCHandler::FixResidualBC
void FixResidualBC(const Vector &x, Vector &y)
y = x-g on ess_tdof_list, the rest of y is unchanged
Definition: petsc.cpp:1783
mfem::Array::LoseData
void LoseData()
NULL-ifies the data.
Definition: array.hpp:104
mfem::PetscParMatrix::operator=
PetscParMatrix & operator=(const PetscParMatrix &B)
Definition: petsc.cpp:468
mfem::PetscParMatrix::EliminateRowsCols
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
mfem::PetscParVector::~PetscParVector
virtual ~PetscParVector()
Calls PETSc&#39;s destroy function.
Definition: petsc.cpp:180
mfem::PetscFieldSplitSolver::PetscFieldSplitSolver
PetscFieldSplitSolver(MPI_Comm comm, Operator &op, const std::string &prefix=std::string())
Definition: petsc.cpp:2540
mfem::PetscSolver::bchandler
PetscBCHandler * bchandler
Handler for boundary conditions.
Definition: petsc.hpp:419
mfem::PetscPreconditionerFactory::NewPreconditioner
virtual Solver * NewPreconditioner(const OperatorHandle &oh)=0
mfem::Array::Assign
void Assign(const T *)
Copy data from a pointer. Size() elements are copied.
Definition: array.hpp:688
mfem::PetscLinearSolver::PetscLinearSolver
PetscLinearSolver(MPI_Comm comm, const std::string &prefix=std::string(), bool wrap=false)
Definition: petsc.cpp:1804
mfem::TimeDependentOperator::isExplicit
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:175
mfem::PetscBCHandler::SetTime
void SetTime(double t)
Sets the current time.
Definition: petsc.hpp:364
mfem::SparseMatrix::GetData
double * GetData() const
Return element data, i.e. array A.
Definition: sparsemat.hpp:135
mfem::SparseMatrix::GetI
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:131
mfem::PetscParMatrix::PetscParMatrix
PetscParMatrix()
Create an empty matrix to be used as a reference to an existing matrix.
Definition: petsc.cpp:401
mfem::PetscSolver::GetNumIterations
int GetNumIterations()
Definition: petsc.cpp:1621
mfem::PetscSolver::SetBCHandler
void SetBCHandler(PetscBCHandler *bch)
Sets the object to handle essential boundary conditions.
Definition: petsc.cpp:1502
mfem::HypreParVector
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:73
mfem::PetscBCHandler::Eval
virtual void Eval(double t, Vector &g)
Boundary conditions evaluation.
Definition: petsc.hpp:354
mfem::Operator::Type
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:123
mfem::PetscBCHandler::PetscBCHandler
PetscBCHandler(Type _type=ZERO)
Definition: petsc.hpp:338
mfem::PetscParMatrix::NNZ
PetscInt NNZ() const
Returns the number of nonzeros.
Definition: petsc.cpp:387
mfem::TimeDependentOperator::ExplicitMult
virtual void ExplicitMult(const Vector &x, Vector &y) const
Perform the action of the explicit part of the operator, G: y = G(x, t) where t is the current time...
Definition: operator.hpp:186
mfem::PetscPreconditioner::Mult
virtual void Mult(const Vector &b, Vector &x) const
Application of the preconditioner.
Definition: petsc.cpp:2197
mfem::Array::SetSize
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:503
mfem::PetscParMatrix::A
Mat A
The actual PETSc object.
Definition: petsc.hpp:137
mfem::PetscParMatrix::N
PetscInt N() const
Returns the global number of columns.
Definition: petsc.cpp:380
mfem::PetscParMatrix::GetY
PetscParVector * GetY() const
Returns the inner vector in the range of A (it creates it if needed)
Definition: petsc.cpp:1010
mfem::PetscODESolver::SetJacobianType
void SetJacobianType(Operator::Type type)
Definition: petsc.cpp:2816
mfem::PetscParMatrix::MakeWrapper
void MakeWrapper(MPI_Comm comm, const Operator *op, Mat *B)
Creates a wrapper around a mfem::Operator op using PETSc&#39;s MATSHELL object and returns the Mat in B...
Definition: petsc.cpp:635
mfem::PetscSolver::operatorset
bool operatorset
Boolean to handle SetOperator calls.
Definition: petsc.hpp:425
mfem::PetscParMatrix::X
PetscParVector * X
Auxiliary vectors for typecasting.
Definition: petsc.hpp:140
mfem::PetscParMatrix::MultTranspose
void MultTranspose(double a, const Vector &x, double b, Vector &y) const
Matvec transpose: y = a A^T x + b y.
Definition: petsc.cpp:1055
mfem::Operator::height
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
mfem::PetscNonlinearSolver::~PetscNonlinearSolver
virtual ~PetscNonlinearSolver()
Definition: petsc.cpp:2608
mfem::PetscODESolver::~PetscODESolver
virtual ~PetscODESolver()
Definition: petsc.cpp:2717
mfem::PetscSolver::GetConverged
int GetConverged()
Definition: petsc.cpp:1588
mfem::PetscSolver::PetscSolver
PetscSolver()
Construct an empty PetscSolver. Initialize protected objects to NULL.
Definition: petsc.cpp:1308
mfem::PetscSolver::SetPreconditionerFactory
void SetPreconditionerFactory(PetscPreconditionerFactory *factory)
Sets the object for the creation of the preconditioner.
Definition: petsc.cpp:1521
mfem::PetscParMatrix::ReleaseMat
Mat ReleaseMat(bool dereference)
Release the PETSc Mat object. If dereference is true, decrement the refcount of the Mat object...
Definition: petsc.cpp:1247
mfem::PetscParVector::PetscParVector
PetscParVector(MPI_Comm comm, PetscInt glob_size, PetscInt *col=NULL)
Creates vector with given global size and partitioning of the columns.
Definition: petsc.cpp:162
mfem::PetscSolverMonitor::MonitorSolution
virtual void MonitorSolution(PetscInt it, PetscReal norm, const Vector &x)
Monitor the solution vector x.
Definition: petsc.hpp:669
mfem::PetscPreconditioner::SetOperator
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: petsc.cpp:2151
mfem::PetscParVector::GlobalSize
PetscInt GlobalSize() const
Returns the global number of rows.
Definition: petsc.cpp:147
mfem::PetscSolver::X
PetscParVector * X
Definition: petsc.hpp:416
mfem::PetscBDDCSolverParams::ess_dof
const Array< int > * ess_dof
Definition: petsc.hpp:547
mfem::PetscParVector::_SetDataAndSize_
void _SetDataAndSize_()
Definition: petsc.cpp:136
mfem::PetscSolver::SetAbsTol
void SetAbsTol(double tol)
Definition: petsc.cpp:1355
mfem::PetscBCHandler
Helper class for handling essential boundary conditions.
Definition: petsc.hpp:328
mfem::PetscParMatrix::GetX
PetscParVector * GetX() const
Returns the inner vector in the domain of A (it creates it if needed)
Definition: petsc.cpp:1000
mfem::PetscBCHandler::GetTDofs
Array< int > & GetTDofs()
Gets essential dofs (local, per-process numbering)
Definition: petsc.hpp:361
mfem::PetscSolver::SetTol
void SetTol(double tol)
Definition: petsc.cpp:1325
mfem::PetscParVector::operator=
PetscParVector & operator=(PetscScalar d)
Set constant values.
Definition: petsc.cpp:298
mfem::ParFiniteElementSpace::GetTrueDofOffsets
HYPRE_Int * GetTrueDofOffsets() const
Definition: pfespace.hpp:175
mfem::PetscSolverMonitor::MonitorResidual
virtual void MonitorResidual(PetscInt it, PetscReal norm, const Vector &r)
Monitor the residual vector r.
Definition: petsc.hpp:675
mfem::BlockOperator::NumRowBlocks
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockoperator.hpp:70
mfem::Vector
Vector data type.
Definition: vector.hpp:41
mfem::FiniteElementSpace::FEColl
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:177
mfem::IdentityOperator
Identity Operator I: x -&gt; x.
Definition: operator.hpp:282
mfem::PetscParVector::PlaceArray
void PlaceArray(PetscScalar *temp_data)
Temporarily replace the data of the PETSc Vec object. To return to the original data array...
Definition: petsc.cpp:310
mfem::PetscSolver::FreePrivateContext
void FreePrivateContext()
Definition: petsc.cpp:1708
mfem::ODESolver::f
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:26
mfem::PetscParVector::GlobalVector
Vector * GlobalVector() const
Returns the global vector in each processor.
Definition: petsc.cpp:273
mfem::PetscODESolver::PetscODESolver
PetscODESolver(MPI_Comm comm, const std::string &prefix=std::string())
Definition: petsc.cpp:2692
mfem::Solver
Base class for solvers.
Definition: operator.hpp:259
mfem::Vector::data
double * data
Definition: vector.hpp:46
mfem::out
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
mfem::Operator
Abstract operator.
Definition: operator.hpp:21
mfem::HypreParMatrix
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
mfem::BlockOperator
A class to handle Block systems in a matrix-free implementation.
Definition: blockoperator.hpp:34
mfem::SparseMatrix::GetJ
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:133
mfem::PetscLinearSolver::Mult
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
Definition: petsc.cpp:2043
mfem::PetscParMatrix::operator*=
void operator*=(double s)
Scale all entries by s: A_scaled = s*A.
Definition: petsc.cpp:1034
mfem::BlockOperator::NumColBlocks
int NumColBlocks() const
Return the number of column blocks.
Definition: blockoperator.hpp:72
mfem::PetscBCHandler::SetTDofs
void SetTDofs(Array< int > &list)
Sets essential dofs (local, per-process numbering)
Definition: petsc.cpp:1735
mfem::Operator::width
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
mfem::PetscSolver::private_ctx
void * private_ctx
Private context for solver.
Definition: petsc.hpp:422
mfem::PetscBCHandler::SetUp
void SetUp(PetscInt n)
SetUp the helper object, where n is the size of the solution vector.
Definition: petsc.cpp:1742
mfem::PetscParMatrix::MakeRef
void MakeRef(const PetscParMatrix &master)
Makes this object a reference to another PetscParMatrix.
Definition: petsc.cpp:990
mfem::TimeDependentOperator::ImplicitMult
virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const
Perform the action of the implicit part of the operator, F: y = F(x, k, t) where t is the current tim...
Definition: operator.hpp:196
Generated on Fri Feb 18 2022 19:42:04 for MFEM by   doxygen 1.8.5