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