petsc_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local Includes
23 #include "libmesh/petsc_vector.h"
24 #include "libmesh/petsc_matrix.h"
25 
26 #ifdef LIBMESH_HAVE_PETSC
27 
29 #include "libmesh/dense_vector.h"
30 #include "libmesh/parallel.h"
31 #include "libmesh/petsc_macro.h"
32 #include "libmesh/utility.h"
33 
34 namespace libMesh
35 {
36 
37 //-----------------------------------------------------------------------
38 // PetscVector members
39 
40 template <typename T>
42 {
43  this->_restore_array();
44  libmesh_assert(this->closed());
45 
46  PetscErrorCode ierr=0;
47  PetscScalar value=0.;
48 
49  ierr = VecSum (_vec, &value);
50  LIBMESH_CHKERR(ierr);
51 
52  return static_cast<T>(value);
53 }
54 
55 
56 template <typename T>
58 {
59  this->_restore_array();
60  libmesh_assert(this->closed());
61 
62  PetscErrorCode ierr=0;
63  PetscReal value=0.;
64 
65  ierr = VecNorm (_vec, NORM_1, &value);
66  LIBMESH_CHKERR(ierr);
67 
68  return static_cast<Real>(value);
69 }
70 
71 
72 
73 template <typename T>
75 {
76  this->_restore_array();
77  libmesh_assert(this->closed());
78 
79  PetscErrorCode ierr=0;
80  PetscReal value=0.;
81 
82  ierr = VecNorm (_vec, NORM_2, &value);
83  LIBMESH_CHKERR(ierr);
84 
85  return static_cast<Real>(value);
86 }
87 
88 
89 
90 
91 template <typename T>
93 {
94  this->_restore_array();
95  libmesh_assert(this->closed());
96 
97  PetscErrorCode ierr=0;
98  PetscReal value=0.;
99 
100  ierr = VecNorm (_vec, NORM_INFINITY, &value);
101  LIBMESH_CHKERR(ierr);
102 
103  return static_cast<Real>(value);
104 }
105 
106 
107 
108 
109 template <typename T>
112 {
113  this->_restore_array();
114  libmesh_assert(this->closed());
115 
116  this->add(1., v);
117 
118  return *this;
119 }
120 
121 
122 
123 template <typename T>
126 {
127  this->_restore_array();
128  libmesh_assert(this->closed());
129 
130  this->add(-1., v);
131 
132  return *this;
133 }
134 
135 
136 
137 template <typename T>
138 void PetscVector<T>::set (const numeric_index_type i, const T value)
139 {
140  this->_restore_array();
141  libmesh_assert_less (i, size());
142 
143  PetscErrorCode ierr=0;
144  PetscInt i_val = static_cast<PetscInt>(i);
145  PetscScalar petsc_value = static_cast<PetscScalar>(value);
146 
147  ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES);
148  LIBMESH_CHKERR(ierr);
149 
150  this->_is_closed = false;
151 }
152 
153 
154 
155 template <typename T>
157 {
158  PetscErrorCode ierr = 0;
159 
160  // VecReciprocal has been in PETSc since at least 2.3.3 days
161  ierr = VecReciprocal(_vec);
162  LIBMESH_CHKERR(ierr);
163 }
164 
165 
166 
167 template <typename T>
169 {
170  PetscErrorCode ierr = 0;
171 
172  // We just call the PETSc VecConjugate
173  ierr = VecConjugate(_vec);
174  LIBMESH_CHKERR(ierr);
175 }
176 
177 
178 
179 template <typename T>
180 void PetscVector<T>::add (const numeric_index_type i, const T value)
181 {
182  this->_restore_array();
183  libmesh_assert_less (i, size());
184 
185  PetscErrorCode ierr=0;
186  PetscInt i_val = static_cast<PetscInt>(i);
187  PetscScalar petsc_value = static_cast<PetscScalar>(value);
188 
189  ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES);
190  LIBMESH_CHKERR(ierr);
191 
192  this->_is_closed = false;
193 }
194 
195 
196 
197 template <typename T>
198 void PetscVector<T>::add_vector (const T * v,
199  const std::vector<numeric_index_type> & dof_indices)
200 {
201  // If we aren't adding anything just return
202  if(dof_indices.empty())
203  return;
204 
205  this->_restore_array();
206 
207  PetscErrorCode ierr=0;
208  const PetscInt * i_val = reinterpret_cast<const PetscInt *>(&dof_indices[0]);
209  const PetscScalar * petsc_value = static_cast<const PetscScalar *>(v);
210 
211  ierr = VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
212  i_val, petsc_value, ADD_VALUES);
213  LIBMESH_CHKERR(ierr);
214 
215  this->_is_closed = false;
216 }
217 
218 
219 
220 template <typename T>
222  const SparseMatrix<T> & A_in)
223 {
224  this->_restore_array();
225  // Make sure the data passed in are really of Petsc types
226  const PetscVector<T> * V = cast_ptr<const PetscVector<T> *>(&V_in);
227  const PetscMatrix<T> * A = cast_ptr<const PetscMatrix<T> *>(&A_in);
228 
229  PetscErrorCode ierr=0;
230 
231  A->close();
232 
233  // The const_cast<> is not elegant, but it is required since PETSc
234  // is not const-correct.
235  ierr = MatMultAdd(const_cast<PetscMatrix<T> *>(A)->mat(), V->_vec, _vec, _vec);
236  LIBMESH_CHKERR(ierr);
237 }
238 
239 
240 
241 template <typename T>
243  const SparseMatrix<T> & A_in)
244 {
245  this->_restore_array();
246  // Make sure the data passed in are really of Petsc types
247  const PetscVector<T> * V = cast_ptr<const PetscVector<T> *>(&V_in);
248  const PetscMatrix<T> * A = cast_ptr<const PetscMatrix<T> *>(&A_in);
249 
250  PetscErrorCode ierr=0;
251 
252  A->close();
253 
254  // The const_cast<> is not elegant, but it is required since PETSc
255  // is not const-correct.
256  ierr = MatMultTransposeAdd(const_cast<PetscMatrix<T> *>(A)->mat(), V->_vec, _vec, _vec);
257  LIBMESH_CHKERR(ierr);
258 }
259 
260 
261 #if PETSC_VERSION_LESS_THAN(3,1,0)
262 template <typename T>
264  const SparseMatrix<T> &)
265 {
266  libmesh_error_msg("MatMultHermitianTranspose was introduced in PETSc 3.1.0," \
267  << "No one has made it backwards compatible with older " \
268  << "versions of PETSc so far.");
269 }
270 
271 #else
272 
273 template <typename T>
275  const SparseMatrix<T> & A_in)
276 {
277  this->_restore_array();
278  // Make sure the data passed in are really of Petsc types
279  const PetscVector<T> * V = cast_ptr<const PetscVector<T> *>(&V_in);
280  const PetscMatrix<T> * A = cast_ptr<const PetscMatrix<T> *>(&A_in);
281 
282  A->close();
283 
284  // Store a temporary copy since MatMultHermitianTransposeAdd doesn't seem to work
285  // TODO: Find out why MatMultHermitianTransposeAdd doesn't work, might be a PETSc bug?
286  UniquePtr< NumericVector<Number> > this_clone = this->clone();
287 
288  // The const_cast<> is not elegant, but it is required since PETSc
289  // is not const-correct.
290  PetscErrorCode ierr = MatMultHermitianTranspose(const_cast<PetscMatrix<T> *>(A)->mat(), V->_vec, _vec);
291  LIBMESH_CHKERR(ierr);
292 
293  // Add the temporary copy to the matvec result
294  this->add(1., *this_clone);
295 }
296 #endif
297 
298 
299 template <typename T>
300 void PetscVector<T>::add (const T v_in)
301 {
302  this->_get_array(false);
303 
304  for (numeric_index_type i=0; i<_local_size; i++)
305  _values[i] += v_in;
306 }
307 
308 
309 
310 template <typename T>
312 {
313  this->add (1., v);
314 }
315 
316 
317 
318 template <typename T>
319 void PetscVector<T>::add (const T a_in, const NumericVector<T> & v_in)
320 {
321  this->_restore_array();
322 
323  PetscErrorCode ierr = 0;
324  PetscScalar a = static_cast<PetscScalar>(a_in);
325 
326  // Make sure the NumericVector passed in is really a PetscVector
327  const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
328  v->_restore_array();
329 
330  libmesh_assert_equal_to (this->size(), v->size());
331 
332  if(this->type() != GHOSTED)
333  {
334  ierr = VecAXPY(_vec, a, v->_vec);
335  LIBMESH_CHKERR(ierr);
336  }
337  else
338  {
339  Vec loc_vec;
340  Vec v_loc_vec;
341  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
342  LIBMESH_CHKERR(ierr);
343  ierr = VecGhostGetLocalForm (v->_vec,&v_loc_vec);
344  LIBMESH_CHKERR(ierr);
345 
346  ierr = VecAXPY(loc_vec, a, v_loc_vec);
347  LIBMESH_CHKERR(ierr);
348 
349  ierr = VecGhostRestoreLocalForm (v->_vec,&v_loc_vec);
350  LIBMESH_CHKERR(ierr);
351  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
352  LIBMESH_CHKERR(ierr);
353  }
354 }
355 
356 
357 
358 template <typename T>
359 void PetscVector<T>::insert (const T * v,
360  const std::vector<numeric_index_type> & dof_indices)
361 {
362  if (dof_indices.empty())
363  return;
364 
365  this->_restore_array();
366 
367  PetscErrorCode ierr=0;
368  PetscInt * idx_values = numeric_petsc_cast(&dof_indices[0]);
369  ierr = VecSetValues (_vec, dof_indices.size(), idx_values, v, INSERT_VALUES);
370  LIBMESH_CHKERR(ierr);
371 
372  this->_is_closed = false;
373 }
374 
375 
376 
377 template <typename T>
378 void PetscVector<T>::scale (const T factor_in)
379 {
380  this->_restore_array();
381 
382  PetscErrorCode ierr = 0;
383  PetscScalar factor = static_cast<PetscScalar>(factor_in);
384 
385  if(this->type() != GHOSTED)
386  {
387  ierr = VecScale(_vec, factor);
388  LIBMESH_CHKERR(ierr);
389  }
390  else
391  {
392  Vec loc_vec;
393  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
394  LIBMESH_CHKERR(ierr);
395 
396  ierr = VecScale(loc_vec, factor);
397  LIBMESH_CHKERR(ierr);
398 
399  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
400  LIBMESH_CHKERR(ierr);
401  }
402 }
403 
404 template <typename T>
406 {
407  PetscErrorCode ierr = 0;
408 
409  const PetscVector<T> * v_vec = cast_ptr<const PetscVector<T> *>(&v);
410 
411  ierr = VecPointwiseDivide(_vec, _vec, v_vec->_vec);
412  LIBMESH_CHKERR(ierr);
413 
414  return *this;
415 }
416 
417 template <typename T>
419 {
420  this->_restore_array();
421 
422  PetscErrorCode ierr = 0;
423 
424  if(this->type() != GHOSTED)
425  {
426  ierr = VecAbs(_vec);
427  LIBMESH_CHKERR(ierr);
428  }
429  else
430  {
431  Vec loc_vec;
432  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
433  LIBMESH_CHKERR(ierr);
434 
435  ierr = VecAbs(loc_vec);
436  LIBMESH_CHKERR(ierr);
437 
438  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
439  LIBMESH_CHKERR(ierr);
440  }
441 }
442 
443 template <typename T>
445 {
446  this->_restore_array();
447 
448  // Error flag
449  PetscErrorCode ierr = 0;
450 
451  // Return value
452  PetscScalar value=0.;
453 
454  // Make sure the NumericVector passed in is really a PetscVector
455  const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&V);
456 
457  // 2.3.x (at least) style. Untested for previous versions.
458  ierr = VecDot(this->_vec, v->_vec, &value);
459  LIBMESH_CHKERR(ierr);
460 
461  return static_cast<T>(value);
462 }
463 
464 template <typename T>
466 {
467  this->_restore_array();
468 
469  // Error flag
470  PetscErrorCode ierr = 0;
471 
472  // Return value
473  PetscScalar value=0.;
474 
475  // Make sure the NumericVector passed in is really a PetscVector
476  const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&V);
477 
478  // 2.3.x (at least) style. Untested for previous versions.
479  ierr = VecTDot(this->_vec, v->_vec, &value);
480  LIBMESH_CHKERR(ierr);
481 
482  return static_cast<T>(value);
483 }
484 
485 
486 template <typename T>
489 {
490  this->_restore_array();
491  libmesh_assert(this->closed());
492 
493  PetscErrorCode ierr = 0;
494  PetscScalar s = static_cast<PetscScalar>(s_in);
495 
496  if (this->size() != 0)
497  {
498  if(this->type() != GHOSTED)
499  {
500  ierr = VecSet(_vec, s);
501  LIBMESH_CHKERR(ierr);
502  }
503  else
504  {
505  Vec loc_vec;
506  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
507  LIBMESH_CHKERR(ierr);
508 
509  ierr = VecSet(loc_vec, s);
510  LIBMESH_CHKERR(ierr);
511 
512  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
513  LIBMESH_CHKERR(ierr);
514  }
515  }
516 
517  return *this;
518 }
519 
520 
521 
522 template <typename T>
525 {
526  // Make sure the NumericVector passed in is really a PetscVector
527  const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
528 
529  *this = *v;
530 
531  return *this;
532 }
533 
534 
535 
536 template <typename T>
539 {
540  this->_restore_array();
541  v._restore_array();
542 
543  libmesh_assert_equal_to (this->size(), v.size());
544  libmesh_assert_equal_to (this->local_size(), v.local_size());
545  libmesh_assert (v.closed());
546 
547  PetscErrorCode ierr = 0;
548 
549  if (((this->type()==PARALLEL) && (v.type()==GHOSTED)) ||
550  ((this->type()==GHOSTED) && (v.type()==PARALLEL)) ||
551  ((this->type()==GHOSTED) && (v.type()==SERIAL)) ||
552  ((this->type()==SERIAL) && (v.type()==GHOSTED)))
553  {
554  /* Allow assignment of a ghosted to a parallel vector since this
555  causes no difficulty. See discussion in libmesh-devel of
556  June 24, 2010. */
557  ierr = VecCopy (v._vec, this->_vec);
558  LIBMESH_CHKERR(ierr);
559  }
560  else
561  {
562  /* In all other cases, we assert that both vectors are of equal
563  type. */
564  libmesh_assert_equal_to (this->_type, v._type);
565 
566  if (v.size() != 0)
567  {
568  if(this->type() != GHOSTED)
569  {
570  ierr = VecCopy (v._vec, this->_vec);
571  LIBMESH_CHKERR(ierr);
572  }
573  else
574  {
575  Vec loc_vec;
576  Vec v_loc_vec;
577  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
578  LIBMESH_CHKERR(ierr);
579  ierr = VecGhostGetLocalForm (v._vec,&v_loc_vec);
580  LIBMESH_CHKERR(ierr);
581 
582  ierr = VecCopy (v_loc_vec, loc_vec);
583  LIBMESH_CHKERR(ierr);
584 
585  ierr = VecGhostRestoreLocalForm (v._vec,&v_loc_vec);
586  LIBMESH_CHKERR(ierr);
587  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
588  LIBMESH_CHKERR(ierr);
589  }
590  }
591  }
592 
593  close();
594 
595  return *this;
596 }
597 
598 
599 
600 template <typename T>
602 PetscVector<T>::operator = (const std::vector<T> & v)
603 {
604  this->_get_array(false);
605 
610  if (this->size() == v.size())
611  {
612  numeric_index_type first = first_local_index();
613  numeric_index_type last = last_local_index();
614  for (numeric_index_type i=0; i<last-first; i++)
615  _values[i] = static_cast<PetscScalar>(v[first + i]);
616  }
617 
622  else
623  {
624  for (numeric_index_type i=0; i<_local_size; i++)
625  _values[i] = static_cast<PetscScalar>(v[i]);
626  }
627 
628  // Make sure ghost dofs are up to date
629  if (this->type() == GHOSTED)
630  this->close();
631 
632  return *this;
633 }
634 
635 
636 
637 template <typename T>
639 {
640  this->_restore_array();
641 
642  // Make sure the NumericVector passed in is really a PetscVector
643  PetscVector<T> * v_local = cast_ptr<PetscVector<T> *>(&v_local_in);
644 
645  libmesh_assert(v_local);
646  libmesh_assert_equal_to (v_local->size(), this->size());
647 
648  PetscErrorCode ierr = 0;
649  const PetscInt n = this->size();
650 
651  IS is;
652  VecScatter scatter;
653 
654  // Create idx, idx[i] = i;
655  std::vector<PetscInt> idx(n); Utility::iota (idx.begin(), idx.end(), 0);
656 
657  // Create the index set & scatter object
658  ierr = ISCreateLibMesh(this->comm().get(), n, &idx[0], PETSC_USE_POINTER, &is);
659  LIBMESH_CHKERR(ierr);
660 
661  ierr = VecScatterCreate(_vec, is,
662  v_local->_vec, is,
663  &scatter);
664  LIBMESH_CHKERR(ierr);
665 
666  // Perform the scatter
667  ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
668  INSERT_VALUES, SCATTER_FORWARD);
669  LIBMESH_CHKERR(ierr);
670 
671  ierr = VecScatterEnd (scatter, _vec, v_local->_vec,
672  INSERT_VALUES, SCATTER_FORWARD);
673  LIBMESH_CHKERR(ierr);
674 
675  // Clean up
676  ierr = LibMeshISDestroy (&is);
677  LIBMESH_CHKERR(ierr);
678 
679  ierr = LibMeshVecScatterDestroy(&scatter);
680  LIBMESH_CHKERR(ierr);
681 
682  // Make sure ghost dofs are up to date
683  if (v_local->type() == GHOSTED)
684  v_local->close();
685 }
686 
687 
688 
689 template <typename T>
691  const std::vector<numeric_index_type> & send_list) const
692 {
693  // FIXME: Workaround for a strange bug at large-scale.
694  // If we have ghosting, PETSc lets us just copy the solution, and
695  // doing so avoids a segfault?
696  if (v_local_in.type() == GHOSTED &&
697  this->type() == PARALLEL)
698  {
699  v_local_in = *this;
700  return;
701  }
702 
703  // Normal code path begins here
704 
705  this->_restore_array();
706 
707  // Make sure the NumericVector passed in is really a PetscVector
708  PetscVector<T> * v_local = cast_ptr<PetscVector<T> *>(&v_local_in);
709 
710  libmesh_assert(v_local);
711  libmesh_assert_equal_to (v_local->size(), this->size());
712  libmesh_assert_less_equal (send_list.size(), v_local->size());
713 
714  PetscErrorCode ierr=0;
715  const numeric_index_type n_sl =
716  cast_int<numeric_index_type>(send_list.size());
717 
718  IS is;
719  VecScatter scatter;
720 
721  std::vector<PetscInt> idx(n_sl + this->local_size());
722 
723  for (numeric_index_type i=0; i<n_sl; i++)
724  idx[i] = static_cast<PetscInt>(send_list[i]);
725  for (numeric_index_type i = 0; i != this->local_size(); ++i)
726  idx[n_sl+i] = i + this->first_local_index();
727 
728  // Create the index set & scatter object
729  if (idx.empty())
730  ierr = ISCreateLibMesh(this->comm().get(),
731  n_sl+this->local_size(), PETSC_NULL, PETSC_USE_POINTER, &is);
732  else
733  ierr = ISCreateLibMesh(this->comm().get(),
734  n_sl+this->local_size(), &idx[0], PETSC_USE_POINTER, &is);
735  LIBMESH_CHKERR(ierr);
736 
737  ierr = VecScatterCreate(_vec, is,
738  v_local->_vec, is,
739  &scatter);
740  LIBMESH_CHKERR(ierr);
741 
742 
743  // Perform the scatter
744  ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
745  INSERT_VALUES, SCATTER_FORWARD);
746  LIBMESH_CHKERR(ierr);
747 
748  ierr = VecScatterEnd (scatter, _vec, v_local->_vec,
749  INSERT_VALUES, SCATTER_FORWARD);
750  LIBMESH_CHKERR(ierr);
751 
752  // Clean up
753  ierr = LibMeshISDestroy (&is);
754  LIBMESH_CHKERR(ierr);
755 
756  ierr = LibMeshVecScatterDestroy(&scatter);
757  LIBMESH_CHKERR(ierr);
758 
759  // Make sure ghost dofs are up to date
760  if (v_local->type() == GHOSTED)
761  v_local->close();
762 }
763 
764 
765 
766 template <typename T>
767 void PetscVector<T>::localize (std::vector<T> & v_local,
768  const std::vector<numeric_index_type> & indices) const
769 {
770  // Error code used to check the status of all PETSc function calls.
771  PetscErrorCode ierr = 0;
772 
773  // Create a sequential destination Vec with the right number of entries on each proc.
774  Vec dest;
775  ierr = VecCreateSeq(PETSC_COMM_SELF, indices.size(), &dest);
776  LIBMESH_CHKERR(ierr);
777 
778  // Create an IS using the libmesh routine. PETSc does not own the
779  // IS memory in this case, it is automatically cleaned up by the
780  // std::vector destructor.
781  IS is;
782  ierr = ISCreateLibMesh(this->comm().get(),
783  indices.size(),
784  numeric_petsc_cast(&indices[0]),
786  &is);
787  LIBMESH_CHKERR(ierr);
788 
789  // Create the VecScatter object. "NULL" means "use the identity IS".
790  VecScatter scatter;
791  ierr = VecScatterCreate(_vec,
792  /*src is=*/is,
793  /*dest vec=*/dest,
794  /*dest is=*/NULL,
795  &scatter);
796  LIBMESH_CHKERR(ierr);
797 
798  // Do the scatter
799  ierr = VecScatterBegin(scatter, _vec, dest, INSERT_VALUES, SCATTER_FORWARD);
800  LIBMESH_CHKERR(ierr);
801 
802  ierr = VecScatterEnd(scatter, _vec, dest, INSERT_VALUES, SCATTER_FORWARD);
803  LIBMESH_CHKERR(ierr);
804 
805  // Get access to the values stored in dest.
806  PetscScalar * values;
807  ierr = VecGetArray (dest, &values);
808  LIBMESH_CHKERR(ierr);
809 
810  // Store values into the provided v_local. Make sure there is enough
811  // space reserved and then clear out any existing entries before
812  // inserting.
813  v_local.reserve(indices.size());
814  v_local.clear();
815  v_local.insert(v_local.begin(), values, values+indices.size());
816 
817  // We are done using it, so restore the array.
818  ierr = VecRestoreArray (dest, &values);
819  LIBMESH_CHKERR(ierr);
820 
821  // Clean up PETSc data structures.
822  ierr = LibMeshVecScatterDestroy(&scatter);
823  LIBMESH_CHKERR(ierr);
824 
825  ierr = LibMeshISDestroy (&is);
826  LIBMESH_CHKERR(ierr);
827 
828  ierr = LibMeshVecDestroy(&dest);
829  LIBMESH_CHKERR(ierr);
830 }
831 
832 
833 
834 template <typename T>
835 void PetscVector<T>::localize (const numeric_index_type first_local_idx,
836  const numeric_index_type last_local_idx,
837  const std::vector<numeric_index_type> & send_list)
838 {
839  this->_restore_array();
840 
841  libmesh_assert_less_equal (send_list.size(), this->size());
842  libmesh_assert_less_equal (last_local_idx+1, this->size());
843 
844  const numeric_index_type my_size = this->size();
845  const numeric_index_type my_local_size = (last_local_idx - first_local_idx + 1);
846  PetscErrorCode ierr=0;
847 
848  // Don't bother for serial cases
849  // if ((first_local_idx == 0) &&
850  // (my_local_size == my_size))
851  // But we do need to stay in sync for degenerate cases
852  if (this->n_processors() == 1)
853  return;
854 
855 
856  // Build a parallel vector, initialize it with the local
857  // parts of (*this)
858  PetscVector<T> parallel_vec(this->comm(), PARALLEL);
859 
860  parallel_vec.init (my_size, my_local_size, true, PARALLEL);
861 
862 
863  // Copy part of *this into the parallel_vec
864  {
865  IS is;
866  VecScatter scatter;
867 
868  // Create idx, idx[i] = i+first_local_idx;
869  std::vector<PetscInt> idx(my_local_size);
870  Utility::iota (idx.begin(), idx.end(), first_local_idx);
871 
872  // Create the index set & scatter object
873  ierr = ISCreateLibMesh(this->comm().get(), my_local_size,
874  my_local_size ? &idx[0] : libmesh_nullptr, PETSC_USE_POINTER, &is);
875  LIBMESH_CHKERR(ierr);
876 
877  ierr = VecScatterCreate(_vec, is,
878  parallel_vec._vec, is,
879  &scatter);
880  LIBMESH_CHKERR(ierr);
881 
882  ierr = VecScatterBegin(scatter, _vec, parallel_vec._vec,
883  INSERT_VALUES, SCATTER_FORWARD);
884  LIBMESH_CHKERR(ierr);
885 
886  ierr = VecScatterEnd (scatter, _vec, parallel_vec._vec,
887  INSERT_VALUES, SCATTER_FORWARD);
888  LIBMESH_CHKERR(ierr);
889 
890  // Clean up
891  ierr = LibMeshISDestroy (&is);
892  LIBMESH_CHKERR(ierr);
893 
894  ierr = LibMeshVecScatterDestroy(&scatter);
895  LIBMESH_CHKERR(ierr);
896  }
897 
898  // localize like normal
899  parallel_vec.close();
900  parallel_vec.localize (*this, send_list);
901  this->close();
902 }
903 
904 
905 
906 template <typename T>
907 void PetscVector<T>::localize (std::vector<T> & v_local) const
908 {
909  this->_restore_array();
910 
911  // This function must be run on all processors at once
912  parallel_object_only();
913 
914  PetscErrorCode ierr=0;
915  const PetscInt n = this->size();
916  const PetscInt nl = this->local_size();
917  PetscScalar * values;
918 
919  v_local.clear();
920  v_local.resize(n, 0.);
921 
922  ierr = VecGetArray (_vec, &values);
923  LIBMESH_CHKERR(ierr);
924 
925  numeric_index_type ioff = first_local_index();
926 
927  for (PetscInt i=0; i<nl; i++)
928  v_local[i+ioff] = static_cast<T>(values[i]);
929 
930  ierr = VecRestoreArray (_vec, &values);
931  LIBMESH_CHKERR(ierr);
932 
933  this->comm().sum(v_local);
934 }
935 
936 
937 
938 // Full specialization for Real datatypes
939 #ifdef LIBMESH_USE_REAL_NUMBERS
940 
941 template <>
942 void PetscVector<Real>::localize_to_one (std::vector<Real> & v_local,
943  const processor_id_type pid) const
944 {
945  this->_restore_array();
946 
947  PetscErrorCode ierr=0;
948  const PetscInt n = size();
949  const PetscInt nl = local_size();
950  PetscScalar * values;
951 
952 
953  // only one processor
954  if (n_processors() == 1)
955  {
956  v_local.resize(n);
957 
958  ierr = VecGetArray (_vec, &values);
959  LIBMESH_CHKERR(ierr);
960 
961  for (PetscInt i=0; i<n; i++)
962  v_local[i] = static_cast<Real>(values[i]);
963 
964  ierr = VecRestoreArray (_vec, &values);
965  LIBMESH_CHKERR(ierr);
966  }
967 
968  // otherwise multiple processors
969  else
970  {
971  if(pid == 0) // optimized version for localizing to 0
972  {
973  Vec vout;
974  VecScatter ctx;
975 
976  ierr = VecScatterCreateToZero(_vec, &ctx, &vout);
977  LIBMESH_CHKERR(ierr);
978 
979  ierr = VecScatterBegin(ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
980  LIBMESH_CHKERR(ierr);
981  ierr = VecScatterEnd(ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
982  LIBMESH_CHKERR(ierr);
983 
984  if(processor_id() == 0)
985  {
986  v_local.resize(n);
987 
988  ierr = VecGetArray (vout, &values);
989  LIBMESH_CHKERR(ierr);
990 
991  for (PetscInt i=0; i<n; i++)
992  v_local[i] = static_cast<Real>(values[i]);
993 
994  ierr = VecRestoreArray (vout, &values);
995  LIBMESH_CHKERR(ierr);
996  }
997 
998  ierr = LibMeshVecScatterDestroy(&ctx);
999  LIBMESH_CHKERR(ierr);
1000  ierr = LibMeshVecDestroy(&vout);
1001  LIBMESH_CHKERR(ierr);
1002 
1003  }
1004  else
1005  {
1006  v_local.resize(n);
1007 
1008  numeric_index_type ioff = this->first_local_index();
1009  std::vector<Real> local_values (n, 0.);
1010 
1011  {
1012  ierr = VecGetArray (_vec, &values);
1013  LIBMESH_CHKERR(ierr);
1014 
1015  for (PetscInt i=0; i<nl; i++)
1016  local_values[i+ioff] = static_cast<Real>(values[i]);
1017 
1018  ierr = VecRestoreArray (_vec, &values);
1019  LIBMESH_CHKERR(ierr);
1020  }
1021 
1022 
1023  MPI_Reduce (&local_values[0], &v_local[0], n, MPI_REAL, MPI_SUM,
1024  pid, this->comm().get());
1025  }
1026  }
1027 }
1028 
1029 #endif
1030 
1031 
1032 // Full specialization for Complex datatypes
1033 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1034 
1035 template <>
1036 void PetscVector<Complex>::localize_to_one (std::vector<Complex> & v_local,
1037  const processor_id_type pid) const
1038 {
1039  this->_restore_array();
1040 
1041  PetscErrorCode ierr=0;
1042  const PetscInt n = size();
1043  const PetscInt nl = local_size();
1044  PetscScalar * values;
1045 
1046 
1047  v_local.resize(n);
1048 
1049 
1050  for (PetscInt i=0; i<n; i++)
1051  v_local[i] = 0.;
1052 
1053  // only one processor
1054  if (n == nl)
1055  {
1056  ierr = VecGetArray (_vec, &values);
1057  LIBMESH_CHKERR(ierr);
1058 
1059  for (PetscInt i=0; i<n; i++)
1060  v_local[i] = static_cast<Complex>(values[i]);
1061 
1062  ierr = VecRestoreArray (_vec, &values);
1063  LIBMESH_CHKERR(ierr);
1064  }
1065 
1066  // otherwise multiple processors
1067  else
1068  {
1069  numeric_index_type ioff = this->first_local_index();
1070 
1071  /* in here the local values are stored, acting as send buffer for MPI
1072  * initialize to zero, since we collect using MPI_SUM
1073  */
1074  std::vector<Real> real_local_values(n, 0.);
1075  std::vector<Real> imag_local_values(n, 0.);
1076 
1077  {
1078  ierr = VecGetArray (_vec, &values);
1079  LIBMESH_CHKERR(ierr);
1080 
1081  // provide my local share to the real and imag buffers
1082  for (PetscInt i=0; i<nl; i++)
1083  {
1084  real_local_values[i+ioff] = static_cast<Complex>(values[i]).real();
1085  imag_local_values[i+ioff] = static_cast<Complex>(values[i]).imag();
1086  }
1087 
1088  ierr = VecRestoreArray (_vec, &values);
1089  LIBMESH_CHKERR(ierr);
1090  }
1091 
1092  /* have buffers of the real and imaginary part of v_local.
1093  * Once MPI_Reduce() collected all the real and imaginary
1094  * parts in these std::vector<Real>, the values can be
1095  * copied to v_local
1096  */
1097  std::vector<Real> real_v_local(n);
1098  std::vector<Real> imag_v_local(n);
1099 
1100  // collect entries from other proc's in real_v_local, imag_v_local
1101  MPI_Reduce (&real_local_values[0], &real_v_local[0], n,
1102  MPI_REAL, MPI_SUM,
1103  pid, this->comm().get());
1104 
1105  MPI_Reduce (&imag_local_values[0], &imag_v_local[0], n,
1106  MPI_REAL, MPI_SUM,
1107  pid, this->comm().get());
1108 
1109  // copy real_v_local and imag_v_local to v_local
1110  for (PetscInt i=0; i<n; i++)
1111  v_local[i] = Complex(real_v_local[i], imag_v_local[i]);
1112  }
1113 }
1114 
1115 #endif
1116 
1117 
1118 
1119 template <typename T>
1121  const NumericVector<T> & vec2)
1122 {
1123  this->_restore_array();
1124 
1125  PetscErrorCode ierr = 0;
1126 
1127  // Convert arguments to PetscVector*.
1128  const PetscVector<T> * vec1_petsc = cast_ptr<const PetscVector<T> *>(&vec1);
1129  const PetscVector<T> * vec2_petsc = cast_ptr<const PetscVector<T> *>(&vec2);
1130 
1131  // Call PETSc function.
1132 
1133  if(this->type() != GHOSTED)
1134  {
1135  ierr = VecPointwiseMult(this->vec(),
1136  const_cast<PetscVector<T> *>(vec1_petsc)->vec(),
1137  const_cast<PetscVector<T> *>(vec2_petsc)->vec());
1138  LIBMESH_CHKERR(ierr);
1139  }
1140  else
1141  {
1142  Vec loc_vec;
1143  Vec v1_loc_vec;
1144  Vec v2_loc_vec;
1145  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
1146  LIBMESH_CHKERR(ierr);
1147  ierr = VecGhostGetLocalForm (const_cast<PetscVector<T> *>(vec1_petsc)->vec(),&v1_loc_vec);
1148  LIBMESH_CHKERR(ierr);
1149  ierr = VecGhostGetLocalForm (const_cast<PetscVector<T> *>(vec2_petsc)->vec(),&v2_loc_vec);
1150  LIBMESH_CHKERR(ierr);
1151 
1152  ierr = VecPointwiseMult(loc_vec,v1_loc_vec,v2_loc_vec);
1153  LIBMESH_CHKERR(ierr);
1154 
1155  ierr = VecGhostRestoreLocalForm (const_cast<PetscVector<T> *>(vec1_petsc)->vec(),&v1_loc_vec);
1156  LIBMESH_CHKERR(ierr);
1157  ierr = VecGhostRestoreLocalForm (const_cast<PetscVector<T> *>(vec2_petsc)->vec(),&v2_loc_vec);
1158  LIBMESH_CHKERR(ierr);
1159  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
1160  LIBMESH_CHKERR(ierr);
1161  }
1162 }
1163 
1164 
1165 
1166 template <typename T>
1167 void PetscVector<T>::print_matlab (const std::string & name) const
1168 {
1169  this->_restore_array();
1170  libmesh_assert (this->closed());
1171 
1172  PetscErrorCode ierr=0;
1173  PetscViewer petsc_viewer;
1174 
1175 
1176  ierr = PetscViewerCreate (this->comm().get(),
1177  &petsc_viewer);
1178  LIBMESH_CHKERR(ierr);
1179 
1184  if (name != "")
1185  {
1186  ierr = PetscViewerASCIIOpen( this->comm().get(),
1187  name.c_str(),
1188  &petsc_viewer);
1189  LIBMESH_CHKERR(ierr);
1190 
1191 #if PETSC_VERSION_LESS_THAN(3,7,0)
1192  ierr = PetscViewerSetFormat (petsc_viewer,
1193  PETSC_VIEWER_ASCII_MATLAB);
1194 #else
1195  ierr = PetscViewerPushFormat (petsc_viewer,
1196  PETSC_VIEWER_ASCII_MATLAB);
1197 #endif
1198 
1199  LIBMESH_CHKERR(ierr);
1200 
1201  ierr = VecView (_vec, petsc_viewer);
1202  LIBMESH_CHKERR(ierr);
1203  }
1204 
1208  else
1209  {
1210 
1211 #if PETSC_VERSION_LESS_THAN(3,7,0)
1212  ierr = PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
1213  PETSC_VIEWER_ASCII_MATLAB);
1214 #else
1215  ierr = PetscViewerPushFormat (PETSC_VIEWER_STDOUT_WORLD,
1216  PETSC_VIEWER_ASCII_MATLAB);
1217 #endif
1218 
1219  LIBMESH_CHKERR(ierr);
1220 
1221  ierr = VecView (_vec, PETSC_VIEWER_STDOUT_WORLD);
1222  LIBMESH_CHKERR(ierr);
1223  }
1224 
1225 
1229  ierr = LibMeshPetscViewerDestroy (&petsc_viewer);
1230  LIBMESH_CHKERR(ierr);
1231 }
1232 
1233 
1234 
1235 
1236 
1237 template <typename T>
1239  const std::vector<numeric_index_type> & rows) const
1240 {
1241  this->_restore_array();
1242 
1243  // PETSc data structures
1244  IS parent_is, subvector_is;
1245  VecScatter scatter;
1246  PetscErrorCode ierr = 0;
1247 
1248  // Make sure the passed in subvector is really a PetscVector
1249  PetscVector<T> * petsc_subvector = cast_ptr<PetscVector<T> *>(&subvector);
1250 
1251  // If the petsc_subvector is already initialized, we assume that the
1252  // user has already allocated the *correct* amount of space for it.
1253  // If not, we use the appropriate PETSc routines to initialize it.
1254  if (!petsc_subvector->initialized())
1255  {
1256  // Initialize the petsc_subvector to have enough space to hold
1257  // the entries which will be scattered into it. Note: such an
1258  // init() function (where we let PETSc decide the number of local
1259  // entries) is not currently offered by the PetscVector
1260  // class. Should we differentiate here between sequential and
1261  // parallel vector creation based on this->n_processors() ?
1262  ierr = VecCreateMPI(this->comm().get(),
1263  PETSC_DECIDE, // n_local
1264  cast_int<PetscInt>(rows.size()), // n_global
1265  &(petsc_subvector->_vec));
1266  LIBMESH_CHKERR(ierr);
1267 
1268  ierr = VecSetFromOptions (petsc_subvector->_vec);
1269  LIBMESH_CHKERR(ierr);
1270 
1271  // Mark the subvector as initialized
1272  petsc_subvector->_is_initialized = true;
1273  }
1274  else
1275  {
1276  petsc_subvector->_restore_array();
1277  }
1278 
1279  // Use iota to fill an array with entries [0,1,2,3,4,...rows.size()]
1280  std::vector<PetscInt> idx(rows.size());
1281  Utility::iota (idx.begin(), idx.end(), 0);
1282 
1283  // Construct index sets
1284  ierr = ISCreateLibMesh(this->comm().get(),
1285  rows.size(),
1286  numeric_petsc_cast(&rows[0]),
1288  &parent_is);
1289  LIBMESH_CHKERR(ierr);
1290 
1291  ierr = ISCreateLibMesh(this->comm().get(),
1292  rows.size(),
1293  &idx[0],
1295  &subvector_is);
1296  LIBMESH_CHKERR(ierr);
1297 
1298  // Construct the scatter object
1299  ierr = VecScatterCreate(this->_vec,
1300  parent_is,
1301  petsc_subvector->_vec,
1302  subvector_is,
1303  &scatter); LIBMESH_CHKERR(ierr);
1304 
1305  // Actually perform the scatter
1306  ierr = VecScatterBegin(scatter,
1307  this->_vec,
1308  petsc_subvector->_vec,
1309  INSERT_VALUES,
1310  SCATTER_FORWARD); LIBMESH_CHKERR(ierr);
1311 
1312  ierr = VecScatterEnd(scatter,
1313  this->_vec,
1314  petsc_subvector->_vec,
1315  INSERT_VALUES,
1316  SCATTER_FORWARD); LIBMESH_CHKERR(ierr);
1317 
1318  // Clean up
1319  ierr = LibMeshISDestroy(&parent_is); LIBMESH_CHKERR(ierr);
1320  ierr = LibMeshISDestroy(&subvector_is); LIBMESH_CHKERR(ierr);
1321  ierr = LibMeshVecScatterDestroy(&scatter); LIBMESH_CHKERR(ierr);
1322 }
1323 
1324 
1325 
1326 template <typename T>
1327 void PetscVector<T>::_get_array(bool read_only) const
1328 {
1329  libmesh_assert (this->initialized());
1330 
1331 #ifdef LIBMESH_HAVE_CXX11_THREAD
1332  std::atomic_thread_fence(std::memory_order_acquire);
1333 #else
1334  Threads::spin_mutex::scoped_lock lock(_petsc_vector_mutex);
1335 #endif
1336 
1337  // If the values have already been retrieved and we're currently
1338  // trying to get a non-read only view (ie read/write) and the
1339  // values are currently read only... then we need to restore
1340  // the array first... and then retrieve it again.
1341  if (_array_is_present && !read_only && _values_read_only)
1342  _restore_array();
1343 
1344  // If we already have a read/write array - and we're trying
1345  // to get a read only array - let's just use the read write
1346  if (_array_is_present && read_only && !_values_read_only)
1347  _read_only_values = _values;
1348 
1349  if (!_array_is_present)
1350  {
1351 #ifdef LIBMESH_HAVE_CXX11_THREAD
1352  std::lock_guard<std::mutex> lock(_petsc_vector_mutex);
1353 #endif
1354  if (!_array_is_present)
1355  {
1356  PetscErrorCode ierr=0;
1357  if(this->type() != GHOSTED)
1358  {
1359 #if PETSC_VERSION_LESS_THAN(3,2,0)
1360  // Vec{Get,Restore}ArrayRead were introduced in PETSc 3.2.0. If you
1361  // have an older PETSc than that, we'll do an ugly
1362  // const_cast and call VecGetArray() instead.
1363  ierr = VecGetArray(_vec, const_cast<PetscScalar **>(&_values));
1364  _values_read_only = false;
1365 #else
1366  if (read_only)
1367  {
1368  ierr = VecGetArrayRead(_vec, &_read_only_values);
1369  _values_read_only = true;
1370  }
1371  else
1372  {
1373  ierr = VecGetArray(_vec, &_values);
1374  _values_read_only = false;
1375  }
1376 #endif
1377  LIBMESH_CHKERR(ierr);
1378  _local_size = this->local_size();
1379  }
1380  else
1381  {
1382  ierr = VecGhostGetLocalForm (_vec,&_local_form);
1383  LIBMESH_CHKERR(ierr);
1384 
1385 #if PETSC_VERSION_LESS_THAN(3,2,0)
1386  // Vec{Get,Restore}ArrayRead were introduced in PETSc 3.2.0. If you
1387  // have an older PETSc than that, we'll do an ugly
1388  // const_cast and call VecGetArray() instead.
1389  ierr = VecGetArray(_local_form, const_cast<PetscScalar **>(&_values));
1390  _values_read_only = false;
1391 #else
1392  if (read_only)
1393  {
1394  ierr = VecGetArrayRead(_local_form, &_read_only_values);
1395  _values_read_only = true;
1396  }
1397  else
1398  {
1399  ierr = VecGetArray(_local_form, &_values);
1400  _values_read_only = false;
1401  }
1402 #endif
1403  LIBMESH_CHKERR(ierr);
1404 
1405  PetscInt my_local_size = 0;
1406  ierr = VecGetLocalSize(_local_form, &my_local_size);
1407  LIBMESH_CHKERR(ierr);
1408  _local_size = static_cast<numeric_index_type>(my_local_size);
1409  }
1410 
1411  { // cache ownership range
1412  PetscInt petsc_first=0, petsc_last=0;
1413  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1414  LIBMESH_CHKERR(ierr);
1415  _first = static_cast<numeric_index_type>(petsc_first);
1416  _last = static_cast<numeric_index_type>(petsc_last);
1417  }
1418 #ifdef LIBMESH_HAVE_CXX11_THREAD
1419  std::atomic_thread_fence(std::memory_order_release);
1420  _array_is_present.store(true, std::memory_order_relaxed);
1421 #else
1422  _array_is_present = true;
1423 #endif
1424  }
1425  }
1426 }
1427 
1428 
1429 
1430 template <typename T>
1432 {
1433  if (_values_manually_retrieved)
1434  libmesh_error_msg("PetscVector values were manually retrieved but are being automatically restored!");
1435 
1436 #ifdef LIBMESH_HAVE_CXX11_THREAD
1437  std::atomic_thread_fence(std::memory_order_acquire);
1438 #else
1439  Threads::spin_mutex::scoped_lock lock(_petsc_vector_mutex);
1440 #endif
1441 
1442  libmesh_assert (this->initialized());
1443  if (_array_is_present)
1444  {
1445 #ifdef LIBMESH_HAVE_CXX11_THREAD
1446  std::lock_guard<std::mutex> lock(_petsc_vector_mutex);
1447 #endif
1448 
1449  if (_array_is_present)
1450  {
1451  PetscErrorCode ierr=0;
1452  if(this->type() != GHOSTED)
1453  {
1454 #if PETSC_VERSION_LESS_THAN(3,2,0)
1455  // Vec{Get,Restore}ArrayRead were introduced in PETSc 3.2.0. If you
1456  // have an older PETSc than that, we'll do an ugly
1457  // const_cast and call VecRestoreArray() instead.
1458  ierr = VecRestoreArray (_vec, const_cast<PetscScalar **>(&_values));
1459 #else
1460  if (_values_read_only)
1461  ierr = VecRestoreArrayRead (_vec, &_read_only_values);
1462  else
1463  ierr = VecRestoreArray (_vec, &_values);
1464 #endif
1465 
1466  LIBMESH_CHKERR(ierr);
1467  _values = libmesh_nullptr;
1468  }
1469  else
1470  {
1471 #if PETSC_VERSION_LESS_THAN(3,2,0)
1472  // Vec{Get,Restore}ArrayRead were introduced in PETSc 3.2.0. If you
1473  // have an older PETSc than that, we'll do an ugly
1474  // const_cast and call VecRestoreArray() instead.
1475  ierr = VecRestoreArray (_local_form, const_cast<PetscScalar **>(&_values));
1476 #else
1477  if (_values_read_only)
1478  ierr = VecRestoreArrayRead (_local_form, &_read_only_values);
1479  else
1480  ierr = VecRestoreArray (_local_form, &_values);
1481 #endif
1482  LIBMESH_CHKERR(ierr);
1483  _values = libmesh_nullptr;
1484  ierr = VecGhostRestoreLocalForm (_vec,&_local_form);
1485  LIBMESH_CHKERR(ierr);
1486  _local_form = libmesh_nullptr;
1487  _local_size = 0;
1488  }
1489 #ifdef LIBMESH_HAVE_CXX11_THREAD
1490  std::atomic_thread_fence(std::memory_order_release);
1491  _array_is_present.store(false, std::memory_order_relaxed);
1492 #else
1493  _array_is_present = false;
1494 #endif
1495  }
1496  }
1497 }
1498 
1499 
1500 
1501 
1502 //------------------------------------------------------------------
1503 // Explicit instantiations
1504 template class PetscVector<Number>;
1505 
1506 } // namespace libMesh
1507 
1508 
1509 
1510 #endif // #ifdef LIBMESH_HAVE_PETSC
T indefinite_dot(const NumericVector< T > &) const
Definition: petsc_vector.C:465
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
void _get_array(bool read_only) const
virtual bool closed() const
bool closed()
Definition: libmesh.C:279
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) libmesh_override
NumericVector interface to PETSc Vec.
Definition: petsc_vector.h:64
void add_vector_conjugate_transpose(const NumericVector< T > &V, const SparseMatrix< T > &A_trans)
Definition: petsc_vector.C:263
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices) libmesh_override
Definition: petsc_vector.C:198
virtual void scale(const T factor) libmesh_override
Definition: petsc_vector.C:378
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType type=AUTOMATIC) libmesh_override
Definition: petsc_vector.h:835
virtual NumericVector< T > & operator+=(const NumericVector< T > &V) libmesh_override
Definition: petsc_vector.C:111
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
virtual bool initialized() const
virtual Real l1_norm() const libmesh_override
Definition: petsc_vector.C:57
void iota(ForwardIter first, ForwardIter last, T value)
Definition: utility.h:58
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscInt * numeric_petsc_cast(const numeric_index_type *p)
virtual Real linfty_norm() const libmesh_override
Definition: petsc_vector.C:92
virtual T dot(const NumericVector< T > &) const libmesh_override
Definition: petsc_vector.C:444
virtual NumericVector< T > & operator/=(NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:405
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual NumericVector< T > & operator=(const T s) libmesh_override
Definition: petsc_vector.C:488
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const libmesh_override
void _restore_array() const
virtual void set(const numeric_index_type i, const T value) libmesh_override
Definition: petsc_vector.C:138
PetscErrorCode Vec Mat Mat void * ctx
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices) libmesh_override
Definition: petsc_vector.C:359
virtual numeric_index_type size() const libmesh_override
virtual void create_subvector(NumericVector< T > &subvector, const std::vector< numeric_index_type > &rows) const libmesh_override
virtual void abs() libmesh_override
Definition: petsc_vector.C:418
std::complex< Real > Complex
virtual void close() libmesh_override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
PetscErrorCode ierr
virtual void localize(std::vector< T > &v_local) const libmesh_override
Definition: petsc_vector.C:907
static PetscErrorCode Mat * A
virtual void add(const numeric_index_type i, const T value) libmesh_override
Definition: petsc_vector.C:180
virtual void print_matlab(const std::string &name="") const libmesh_override
virtual Real l2_norm() const libmesh_override
Definition: petsc_vector.C:74
virtual numeric_index_type local_size() const libmesh_override
ParallelType type() const
virtual NumericVector< T > & operator-=(const NumericVector< T > &V) libmesh_override
Definition: petsc_vector.C:125
SparseMatrix interface to PETSc Mat.
bool initialized()
Definition: libmesh.C:272
virtual void add_vector_transpose(const NumericVector< T > &V, const SparseMatrix< T > &A_trans) libmesh_override
Definition: petsc_vector.C:242
virtual void reciprocal() libmesh_override
Definition: petsc_vector.C:156
virtual void conjugate() libmesh_override
Definition: petsc_vector.C:168
processor_id_type processor_id()
Definition: libmesh_base.h:96
virtual T sum() const libmesh_override
Definition: petsc_vector.C:41
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
processor_id_type n_processors()
Definition: libmesh_base.h:88