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  // 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_deprecated();
235  libmesh_warning("Matrix A must be assembled before calling PetscVector::add_vector(v, A).\n"
236  "Please update your code, as this warning will become an error in a future release.");
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_deprecated();
263  libmesh_warning("Matrix A must be assembled before calling PetscVector::add_vector_transpose(v, A).\n"
264  "Please update your code, as this warning will become an error in a future release.");
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_deprecated();
300  libmesh_warning("Matrix A must be assembled before calling PetscVector::add_vector_conjugate_transpose(v, A).\n"
301  "Please update your code, as this warning will become an error in a future release.");
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  UniquePtr< 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[0]);
390  ierr = VecSetValues (_vec, dof_indices.size(), idx_values, v, INSERT_VALUES);
391  LIBMESH_CHKERR(ierr);
392 
393  this->_is_closed = false;
394 }
395 
396 
397 
398 template <typename T>
399 void PetscVector<T>::scale (const T factor_in)
400 {
401  this->_restore_array();
402 
403  PetscErrorCode ierr = 0;
404  PetscScalar factor = static_cast<PetscScalar>(factor_in);
405 
406  if (this->type() != GHOSTED)
407  {
408  ierr = VecScale(_vec, factor);
409  LIBMESH_CHKERR(ierr);
410  }
411  else
412  {
413  Vec loc_vec;
414  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
415  LIBMESH_CHKERR(ierr);
416 
417  ierr = VecScale(loc_vec, factor);
418  LIBMESH_CHKERR(ierr);
419 
420  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
421  LIBMESH_CHKERR(ierr);
422  }
423 }
424 
425 template <typename T>
427 {
428  PetscErrorCode ierr = 0;
429 
430  const PetscVector<T> * v_vec = cast_ptr<const PetscVector<T> *>(&v);
431 
432  ierr = VecPointwiseDivide(_vec, _vec, v_vec->_vec);
433  LIBMESH_CHKERR(ierr);
434 
435  return *this;
436 }
437 
438 template <typename T>
440 {
441  this->_restore_array();
442 
443  PetscErrorCode ierr = 0;
444 
445  if (this->type() != GHOSTED)
446  {
447  ierr = VecAbs(_vec);
448  LIBMESH_CHKERR(ierr);
449  }
450  else
451  {
452  Vec loc_vec;
453  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
454  LIBMESH_CHKERR(ierr);
455 
456  ierr = VecAbs(loc_vec);
457  LIBMESH_CHKERR(ierr);
458 
459  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
460  LIBMESH_CHKERR(ierr);
461  }
462 }
463 
464 template <typename T>
465 T PetscVector<T>::dot (const NumericVector<T> & v_in) const
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_in);
477 
478  // 2.3.x (at least) style. Untested for previous versions.
479  ierr = VecDot(this->_vec, v->_vec, &value);
480  LIBMESH_CHKERR(ierr);
481 
482  return static_cast<T>(value);
483 }
484 
485 template <typename T>
487 {
488  this->_restore_array();
489 
490  // Error flag
491  PetscErrorCode ierr = 0;
492 
493  // Return value
494  PetscScalar value=0.;
495 
496  // Make sure the NumericVector passed in is really a PetscVector
497  const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
498 
499  // 2.3.x (at least) style. Untested for previous versions.
500  ierr = VecTDot(this->_vec, v->_vec, &value);
501  LIBMESH_CHKERR(ierr);
502 
503  return static_cast<T>(value);
504 }
505 
506 
507 template <typename T>
510 {
511  this->_restore_array();
512  libmesh_assert(this->closed());
513 
514  PetscErrorCode ierr = 0;
515  PetscScalar s = static_cast<PetscScalar>(s_in);
516 
517  if (this->size() != 0)
518  {
519  if (this->type() != GHOSTED)
520  {
521  ierr = VecSet(_vec, s);
522  LIBMESH_CHKERR(ierr);
523  }
524  else
525  {
526  Vec loc_vec;
527  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
528  LIBMESH_CHKERR(ierr);
529 
530  ierr = VecSet(loc_vec, s);
531  LIBMESH_CHKERR(ierr);
532 
533  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
534  LIBMESH_CHKERR(ierr);
535  }
536  }
537 
538  return *this;
539 }
540 
541 
542 
543 template <typename T>
546 {
547  // Make sure the NumericVector passed in is really a PetscVector
548  const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
549 
550  *this = *v;
551 
552  return *this;
553 }
554 
555 
556 
557 template <typename T>
560 {
561  this->_restore_array();
562  v._restore_array();
563 
564  libmesh_assert_equal_to (this->size(), v.size());
565  libmesh_assert_equal_to (this->local_size(), v.local_size());
566  libmesh_assert (v.closed());
567 
568  PetscErrorCode ierr = 0;
569 
570  if (((this->type()==PARALLEL) && (v.type()==GHOSTED)) ||
571  ((this->type()==GHOSTED) && (v.type()==PARALLEL)) ||
572  ((this->type()==GHOSTED) && (v.type()==SERIAL)) ||
573  ((this->type()==SERIAL) && (v.type()==GHOSTED)))
574  {
575  /* Allow assignment of a ghosted to a parallel vector since this
576  causes no difficulty. See discussion in libmesh-devel of
577  June 24, 2010. */
578  ierr = VecCopy (v._vec, this->_vec);
579  LIBMESH_CHKERR(ierr);
580  }
581  else
582  {
583  /* In all other cases, we assert that both vectors are of equal
584  type. */
585  libmesh_assert_equal_to (this->_type, v._type);
586 
587  if (v.size() != 0)
588  {
589  if (this->type() != GHOSTED)
590  {
591  ierr = VecCopy (v._vec, this->_vec);
592  LIBMESH_CHKERR(ierr);
593  }
594  else
595  {
596  Vec loc_vec;
597  Vec v_loc_vec;
598  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
599  LIBMESH_CHKERR(ierr);
600  ierr = VecGhostGetLocalForm (v._vec,&v_loc_vec);
601  LIBMESH_CHKERR(ierr);
602 
603  ierr = VecCopy (v_loc_vec, loc_vec);
604  LIBMESH_CHKERR(ierr);
605 
606  ierr = VecGhostRestoreLocalForm (v._vec,&v_loc_vec);
607  LIBMESH_CHKERR(ierr);
608  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
609  LIBMESH_CHKERR(ierr);
610  }
611  }
612  }
613 
614  close();
615 
616  return *this;
617 }
618 
619 
620 
621 template <typename T>
623 PetscVector<T>::operator = (const std::vector<T> & v)
624 {
625  this->_get_array(false);
626 
631  if (this->size() == v.size())
632  {
633  numeric_index_type first = first_local_index();
634  numeric_index_type last = last_local_index();
635  for (numeric_index_type i=0; i<last-first; i++)
636  _values[i] = static_cast<PetscScalar>(v[first + i]);
637  }
638 
643  else
644  {
645  for (numeric_index_type i=0; i<_local_size; i++)
646  _values[i] = static_cast<PetscScalar>(v[i]);
647  }
648 
649  // Make sure ghost dofs are up to date
650  if (this->type() == GHOSTED)
651  this->close();
652 
653  return *this;
654 }
655 
656 
657 
658 template <typename T>
660 {
661  this->_restore_array();
662 
663  // Make sure the NumericVector passed in is really a PetscVector
664  PetscVector<T> * v_local = cast_ptr<PetscVector<T> *>(&v_local_in);
665 
666  libmesh_assert(v_local);
667  libmesh_assert_equal_to (v_local->size(), this->size());
668 
669  PetscErrorCode ierr = 0;
670  const PetscInt n = this->size();
671 
672  IS is;
673  VecScatter scatter;
674 
675  // Create idx, idx[i] = i;
676  std::vector<PetscInt> idx(n); Utility::iota (idx.begin(), idx.end(), 0);
677 
678  // Create the index set & scatter object
679  ierr = ISCreateLibMesh(this->comm().get(), n, &idx[0], PETSC_USE_POINTER, &is);
680  LIBMESH_CHKERR(ierr);
681 
682  ierr = VecScatterCreate(_vec, is,
683  v_local->_vec, is,
684  &scatter);
685  LIBMESH_CHKERR(ierr);
686 
687  // Perform the scatter
688  ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
689  INSERT_VALUES, SCATTER_FORWARD);
690  LIBMESH_CHKERR(ierr);
691 
692  ierr = VecScatterEnd (scatter, _vec, v_local->_vec,
693  INSERT_VALUES, SCATTER_FORWARD);
694  LIBMESH_CHKERR(ierr);
695 
696  // Clean up
697  ierr = LibMeshISDestroy (&is);
698  LIBMESH_CHKERR(ierr);
699 
700  ierr = LibMeshVecScatterDestroy(&scatter);
701  LIBMESH_CHKERR(ierr);
702 
703  // Make sure ghost dofs are up to date
704  if (v_local->type() == GHOSTED)
705  v_local->close();
706 }
707 
708 
709 
710 template <typename T>
712  const std::vector<numeric_index_type> & send_list) const
713 {
714  // FIXME: Workaround for a strange bug at large-scale.
715  // If we have ghosting, PETSc lets us just copy the solution, and
716  // doing so avoids a segfault?
717  if (v_local_in.type() == GHOSTED &&
718  this->type() == PARALLEL)
719  {
720  v_local_in = *this;
721  return;
722  }
723 
724  // Normal code path begins here
725 
726  this->_restore_array();
727 
728  // Make sure the NumericVector passed in is really a PetscVector
729  PetscVector<T> * v_local = cast_ptr<PetscVector<T> *>(&v_local_in);
730 
731  libmesh_assert(v_local);
732  libmesh_assert_equal_to (v_local->size(), this->size());
733  libmesh_assert_less_equal (send_list.size(), v_local->size());
734 
735  PetscErrorCode ierr=0;
736  const numeric_index_type n_sl =
737  cast_int<numeric_index_type>(send_list.size());
738 
739  IS is;
740  VecScatter scatter;
741 
742  std::vector<PetscInt> idx(n_sl + this->local_size());
743 
744  for (numeric_index_type i=0; i<n_sl; i++)
745  idx[i] = static_cast<PetscInt>(send_list[i]);
746  for (numeric_index_type i = 0; i != this->local_size(); ++i)
747  idx[n_sl+i] = i + this->first_local_index();
748 
749  // Create the index set & scatter object
750  if (idx.empty())
751  ierr = ISCreateLibMesh(this->comm().get(),
752  n_sl+this->local_size(), PETSC_NULL, PETSC_USE_POINTER, &is);
753  else
754  ierr = ISCreateLibMesh(this->comm().get(),
755  n_sl+this->local_size(), &idx[0], PETSC_USE_POINTER, &is);
756  LIBMESH_CHKERR(ierr);
757 
758  ierr = VecScatterCreate(_vec, is,
759  v_local->_vec, is,
760  &scatter);
761  LIBMESH_CHKERR(ierr);
762 
763 
764  // Perform the scatter
765  ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
766  INSERT_VALUES, SCATTER_FORWARD);
767  LIBMESH_CHKERR(ierr);
768 
769  ierr = VecScatterEnd (scatter, _vec, v_local->_vec,
770  INSERT_VALUES, SCATTER_FORWARD);
771  LIBMESH_CHKERR(ierr);
772 
773  // Clean up
774  ierr = LibMeshISDestroy (&is);
775  LIBMESH_CHKERR(ierr);
776 
777  ierr = LibMeshVecScatterDestroy(&scatter);
778  LIBMESH_CHKERR(ierr);
779 
780  // Make sure ghost dofs are up to date
781  if (v_local->type() == GHOSTED)
782  v_local->close();
783 }
784 
785 
786 
787 template <typename T>
788 void PetscVector<T>::localize (std::vector<T> & v_local,
789  const std::vector<numeric_index_type> & indices) const
790 {
791  // Error code used to check the status of all PETSc function calls.
792  PetscErrorCode ierr = 0;
793 
794  // Create a sequential destination Vec with the right number of entries on each proc.
795  Vec dest;
796  ierr = VecCreateSeq(PETSC_COMM_SELF, indices.size(), &dest);
797  LIBMESH_CHKERR(ierr);
798 
799  // Create an IS using the libmesh routine. PETSc does not own the
800  // IS memory in this case, it is automatically cleaned up by the
801  // std::vector destructor.
802  IS is;
803  ierr = ISCreateLibMesh(this->comm().get(),
804  indices.size(),
805  numeric_petsc_cast(&indices[0]),
807  &is);
808  LIBMESH_CHKERR(ierr);
809 
810  // Create the VecScatter object. "NULL" means "use the identity IS".
811  VecScatter scatter;
812  ierr = VecScatterCreate(_vec,
813  /*src is=*/is,
814  /*dest vec=*/dest,
815  /*dest is=*/NULL,
816  &scatter);
817  LIBMESH_CHKERR(ierr);
818 
819  // Do the scatter
820  ierr = VecScatterBegin(scatter, _vec, dest, INSERT_VALUES, SCATTER_FORWARD);
821  LIBMESH_CHKERR(ierr);
822 
823  ierr = VecScatterEnd(scatter, _vec, dest, INSERT_VALUES, SCATTER_FORWARD);
824  LIBMESH_CHKERR(ierr);
825 
826  // Get access to the values stored in dest.
827  PetscScalar * values;
828  ierr = VecGetArray (dest, &values);
829  LIBMESH_CHKERR(ierr);
830 
831  // Store values into the provided v_local. Make sure there is enough
832  // space reserved and then clear out any existing entries before
833  // inserting.
834  v_local.reserve(indices.size());
835  v_local.clear();
836  v_local.insert(v_local.begin(), values, values+indices.size());
837 
838  // We are done using it, so restore the array.
839  ierr = VecRestoreArray (dest, &values);
840  LIBMESH_CHKERR(ierr);
841 
842  // Clean up PETSc data structures.
843  ierr = LibMeshVecScatterDestroy(&scatter);
844  LIBMESH_CHKERR(ierr);
845 
846  ierr = LibMeshISDestroy (&is);
847  LIBMESH_CHKERR(ierr);
848 
849  ierr = LibMeshVecDestroy(&dest);
850  LIBMESH_CHKERR(ierr);
851 }
852 
853 
854 
855 template <typename T>
856 void PetscVector<T>::localize (const numeric_index_type first_local_idx,
857  const numeric_index_type last_local_idx,
858  const std::vector<numeric_index_type> & send_list)
859 {
860  this->_restore_array();
861 
862  libmesh_assert_less_equal (send_list.size(), this->size());
863  libmesh_assert_less_equal (last_local_idx+1, this->size());
864 
865  const numeric_index_type my_size = this->size();
866  const numeric_index_type my_local_size = (last_local_idx - first_local_idx + 1);
867  PetscErrorCode ierr=0;
868 
869  // Don't bother for serial cases
870  // if ((first_local_idx == 0) &&
871  // (my_local_size == my_size))
872  // But we do need to stay in sync for degenerate cases
873  if (this->n_processors() == 1)
874  return;
875 
876 
877  // Build a parallel vector, initialize it with the local
878  // parts of (*this)
879  PetscVector<T> parallel_vec(this->comm(), PARALLEL);
880 
881  parallel_vec.init (my_size, my_local_size, true, PARALLEL);
882 
883 
884  // Copy part of *this into the parallel_vec
885  {
886  IS is;
887  VecScatter scatter;
888 
889  // Create idx, idx[i] = i+first_local_idx;
890  std::vector<PetscInt> idx(my_local_size);
891  Utility::iota (idx.begin(), idx.end(), first_local_idx);
892 
893  // Create the index set & scatter object
894  ierr = ISCreateLibMesh(this->comm().get(), my_local_size,
895  my_local_size ? &idx[0] : libmesh_nullptr, PETSC_USE_POINTER, &is);
896  LIBMESH_CHKERR(ierr);
897 
898  ierr = VecScatterCreate(_vec, is,
899  parallel_vec._vec, is,
900  &scatter);
901  LIBMESH_CHKERR(ierr);
902 
903  ierr = VecScatterBegin(scatter, _vec, parallel_vec._vec,
904  INSERT_VALUES, SCATTER_FORWARD);
905  LIBMESH_CHKERR(ierr);
906 
907  ierr = VecScatterEnd (scatter, _vec, parallel_vec._vec,
908  INSERT_VALUES, SCATTER_FORWARD);
909  LIBMESH_CHKERR(ierr);
910 
911  // Clean up
912  ierr = LibMeshISDestroy (&is);
913  LIBMESH_CHKERR(ierr);
914 
915  ierr = LibMeshVecScatterDestroy(&scatter);
916  LIBMESH_CHKERR(ierr);
917  }
918 
919  // localize like normal
920  parallel_vec.close();
921  parallel_vec.localize (*this, send_list);
922  this->close();
923 }
924 
925 
926 
927 template <typename T>
928 void PetscVector<T>::localize (std::vector<T> & v_local) const
929 {
930  this->_restore_array();
931 
932  // This function must be run on all processors at once
933  parallel_object_only();
934 
935  PetscErrorCode ierr=0;
936  const PetscInt n = this->size();
937  const PetscInt nl = this->local_size();
938  PetscScalar * values;
939 
940  v_local.clear();
941  v_local.resize(n, 0.);
942 
943  ierr = VecGetArray (_vec, &values);
944  LIBMESH_CHKERR(ierr);
945 
946  numeric_index_type ioff = first_local_index();
947 
948  for (PetscInt i=0; i<nl; i++)
949  v_local[i+ioff] = static_cast<T>(values[i]);
950 
951  ierr = VecRestoreArray (_vec, &values);
952  LIBMESH_CHKERR(ierr);
953 
954  this->comm().sum(v_local);
955 }
956 
957 
958 
959 // Full specialization for Real datatypes
960 #ifdef LIBMESH_USE_REAL_NUMBERS
961 
962 template <>
963 void PetscVector<Real>::localize_to_one (std::vector<Real> & v_local,
964  const processor_id_type pid) const
965 {
966  this->_restore_array();
967 
968  PetscErrorCode ierr=0;
969  const PetscInt n = size();
970  const PetscInt nl = local_size();
971  PetscScalar * values;
972 
973 
974  // only one processor
975  if (n_processors() == 1)
976  {
977  v_local.resize(n);
978 
979  ierr = VecGetArray (_vec, &values);
980  LIBMESH_CHKERR(ierr);
981 
982  for (PetscInt i=0; i<n; i++)
983  v_local[i] = static_cast<Real>(values[i]);
984 
985  ierr = VecRestoreArray (_vec, &values);
986  LIBMESH_CHKERR(ierr);
987  }
988 
989  // otherwise multiple processors
990  else
991  {
992  if (pid == 0) // optimized version for localizing to 0
993  {
994  Vec vout;
995  VecScatter ctx;
996 
997  ierr = VecScatterCreateToZero(_vec, &ctx, &vout);
998  LIBMESH_CHKERR(ierr);
999 
1000  ierr = VecScatterBegin(ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
1001  LIBMESH_CHKERR(ierr);
1002  ierr = VecScatterEnd(ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
1003  LIBMESH_CHKERR(ierr);
1004 
1005  if (processor_id() == 0)
1006  {
1007  v_local.resize(n);
1008 
1009  ierr = VecGetArray (vout, &values);
1010  LIBMESH_CHKERR(ierr);
1011 
1012  for (PetscInt i=0; i<n; i++)
1013  v_local[i] = static_cast<Real>(values[i]);
1014 
1015  ierr = VecRestoreArray (vout, &values);
1016  LIBMESH_CHKERR(ierr);
1017  }
1018 
1019  ierr = LibMeshVecScatterDestroy(&ctx);
1020  LIBMESH_CHKERR(ierr);
1021  ierr = LibMeshVecDestroy(&vout);
1022  LIBMESH_CHKERR(ierr);
1023 
1024  }
1025  else
1026  {
1027  v_local.resize(n);
1028 
1029  numeric_index_type ioff = this->first_local_index();
1030  std::vector<Real> local_values (n, 0.);
1031 
1032  {
1033  ierr = VecGetArray (_vec, &values);
1034  LIBMESH_CHKERR(ierr);
1035 
1036  for (PetscInt i=0; i<nl; i++)
1037  local_values[i+ioff] = static_cast<Real>(values[i]);
1038 
1039  ierr = VecRestoreArray (_vec, &values);
1040  LIBMESH_CHKERR(ierr);
1041  }
1042 
1043 
1044  MPI_Reduce (&local_values[0], &v_local[0], n, MPI_REAL, MPI_SUM,
1045  pid, this->comm().get());
1046  }
1047  }
1048 }
1049 
1050 #endif
1051 
1052 
1053 // Full specialization for Complex datatypes
1054 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1055 
1056 template <>
1057 void PetscVector<Complex>::localize_to_one (std::vector<Complex> & v_local,
1058  const processor_id_type pid) const
1059 {
1060  this->_restore_array();
1061 
1062  PetscErrorCode ierr=0;
1063  const PetscInt n = size();
1064  const PetscInt nl = local_size();
1065  PetscScalar * values;
1066 
1067 
1068  v_local.resize(n);
1069 
1070 
1071  for (PetscInt i=0; i<n; i++)
1072  v_local[i] = 0.;
1073 
1074  // only one processor
1075  if (n == nl)
1076  {
1077  ierr = VecGetArray (_vec, &values);
1078  LIBMESH_CHKERR(ierr);
1079 
1080  for (PetscInt i=0; i<n; i++)
1081  v_local[i] = static_cast<Complex>(values[i]);
1082 
1083  ierr = VecRestoreArray (_vec, &values);
1084  LIBMESH_CHKERR(ierr);
1085  }
1086 
1087  // otherwise multiple processors
1088  else
1089  {
1090  numeric_index_type ioff = this->first_local_index();
1091 
1092  /* in here the local values are stored, acting as send buffer for MPI
1093  * initialize to zero, since we collect using MPI_SUM
1094  */
1095  std::vector<Real> real_local_values(n, 0.);
1096  std::vector<Real> imag_local_values(n, 0.);
1097 
1098  {
1099  ierr = VecGetArray (_vec, &values);
1100  LIBMESH_CHKERR(ierr);
1101 
1102  // provide my local share to the real and imag buffers
1103  for (PetscInt i=0; i<nl; i++)
1104  {
1105  real_local_values[i+ioff] = static_cast<Complex>(values[i]).real();
1106  imag_local_values[i+ioff] = static_cast<Complex>(values[i]).imag();
1107  }
1108 
1109  ierr = VecRestoreArray (_vec, &values);
1110  LIBMESH_CHKERR(ierr);
1111  }
1112 
1113  /* have buffers of the real and imaginary part of v_local.
1114  * Once MPI_Reduce() collected all the real and imaginary
1115  * parts in these std::vector<Real>, the values can be
1116  * copied to v_local
1117  */
1118  std::vector<Real> real_v_local(n);
1119  std::vector<Real> imag_v_local(n);
1120 
1121  // collect entries from other proc's in real_v_local, imag_v_local
1122  MPI_Reduce (&real_local_values[0], &real_v_local[0], n,
1123  MPI_REAL, MPI_SUM,
1124  pid, this->comm().get());
1125 
1126  MPI_Reduce (&imag_local_values[0], &imag_v_local[0], 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  Utility::iota (idx.begin(), idx.end(), 0);
1303 
1304  // Construct index sets
1305  ierr = ISCreateLibMesh(this->comm().get(),
1306  rows.size(),
1307  numeric_petsc_cast(&rows[0]),
1309  &parent_is);
1310  LIBMESH_CHKERR(ierr);
1311 
1312  ierr = ISCreateLibMesh(this->comm().get(),
1313  rows.size(),
1314  &idx[0],
1316  &subvector_is);
1317  LIBMESH_CHKERR(ierr);
1318 
1319  // Construct the scatter object
1320  ierr = VecScatterCreate(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 = libmesh_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 = libmesh_nullptr;
1505  ierr = VecGhostRestoreLocalForm (_vec,&_local_form);
1506  LIBMESH_CHKERR(ierr);
1507  _local_form = libmesh_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: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)
Definition: petsc_vector.C:277
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:125
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:399
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:625
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
virtual T dot(const NumericVector< T > &v) const libmesh_override
Definition: petsc_vector.C:465
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:57
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 NumericVector< T > & operator/=(NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:426
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual NumericVector< T > & operator=(const T s) libmesh_override
Definition: petsc_vector.C:509
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:380
virtual numeric_index_type size() const libmesh_override
Definition: petsc_vector.h:919
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) libmesh_override
Definition: petsc_vector.C:111
T indefinite_dot(const NumericVector< T > &v) const
Definition: petsc_vector.C:486
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) libmesh_override
Definition: petsc_vector.C:249
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:439
std::complex< Real > Complex
virtual void close() libmesh_override
Definition: petsc_vector.h:807
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:928
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
Definition: petsc_vector.h:939
ParallelType type() const
SparseMatrix interface to PETSc Mat.
bool initialized()
Definition: libmesh.C:272
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