Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : #include "libmesh/petsc_vector.h"
19 :
20 : #ifdef LIBMESH_HAVE_PETSC
21 :
22 : // libMesh includes
23 : #include "libmesh/petsc_matrix_base.h"
24 : #include "libmesh/dense_subvector.h"
25 : #include "libmesh/dense_vector.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/petsc_macro.h"
28 : #include "libmesh/wrapped_petsc.h"
29 :
30 : // TIMPI includes
31 : #include "timpi/op_function.h"
32 : #include "timpi/parallel_implementation.h"
33 : #include "timpi/standard_type.h"
34 :
35 : // C++ includes
36 : #include <numeric> // std::iota
37 :
38 : namespace libMesh
39 : {
40 : //-----------------------------------------------------------------------
41 : // PetscVector members
42 :
43 : template <typename T>
44 140 : T PetscVector<T>::sum () const
45 : {
46 140 : this->_restore_array();
47 4 : libmesh_assert(this->closed());
48 :
49 140 : PetscScalar value=0.;
50 :
51 140 : LibmeshPetscCall(VecSum (_vec, &value));
52 :
53 140 : return static_cast<T>(value);
54 : }
55 :
56 :
57 :
58 : template <typename T>
59 : template <NormType N>
60 263751 : Real PetscVector<T>::norm () const
61 : {
62 7654 : parallel_object_only();
63 :
64 263751 : this->_restore_array();
65 7654 : libmesh_assert(this->closed());
66 :
67 263751 : PetscReal value=0.;
68 :
69 263751 : LibmeshPetscCall(VecNorm (_vec, N, &value));
70 :
71 263751 : return static_cast<Real>(value);
72 : }
73 : template <typename T>
74 1540 : Real PetscVector<T>::l1_norm () const
75 : {
76 1540 : return PetscVector<T>::norm<NORM_1>();
77 : }
78 : template <typename T>
79 262001 : Real PetscVector<T>::l2_norm () const
80 : {
81 262001 : return PetscVector<T>::norm<NORM_2>();
82 : }
83 : template <typename T>
84 210 : Real PetscVector<T>::linfty_norm () const
85 : {
86 210 : return PetscVector<T>::norm<NORM_INFINITY>();
87 : }
88 :
89 :
90 :
91 : template <typename T>
92 : NumericVector<T> &
93 32830 : PetscVector<T>::operator += (const NumericVector<T> & v)
94 : {
95 938 : parallel_object_only();
96 :
97 32830 : this->_restore_array();
98 938 : libmesh_assert(this->closed());
99 :
100 32830 : this->add(1., v);
101 :
102 32830 : return *this;
103 : }
104 :
105 :
106 :
107 : template <typename T>
108 : NumericVector<T> &
109 80780 : PetscVector<T>::operator -= (const NumericVector<T> & v)
110 : {
111 2308 : parallel_object_only();
112 :
113 80780 : this->_restore_array();
114 2308 : libmesh_assert(this->closed());
115 :
116 80780 : this->add(-1., v);
117 :
118 80780 : return *this;
119 : }
120 :
121 :
122 :
123 : template <typename T>
124 44138174 : void PetscVector<T>::set (const numeric_index_type i, const T value)
125 : {
126 44138174 : this->_restore_array();
127 4052344 : libmesh_assert_less (i, size());
128 :
129 44138174 : PetscInt i_val = static_cast<PetscInt>(i);
130 44138174 : PetscScalar petsc_value = PS(value);
131 :
132 44138174 : std::scoped_lock lock(this->_numeric_vector_mutex);
133 44138174 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES));
134 :
135 44138174 : this->_is_closed = false;
136 44138174 : }
137 :
138 :
139 :
140 : template <typename T>
141 140 : void PetscVector<T>::reciprocal()
142 : {
143 4 : parallel_object_only();
144 :
145 :
146 : // VecReciprocal has been in PETSc since at least 2.3.3 days
147 140 : LibmeshPetscCall(VecReciprocal(_vec));
148 140 : }
149 :
150 :
151 :
152 : template <typename T>
153 0 : void PetscVector<T>::conjugate()
154 : {
155 0 : parallel_object_only();
156 :
157 :
158 : // We just call the PETSc VecConjugate
159 0 : LibmeshPetscCall(VecConjugate(_vec));
160 0 : }
161 :
162 :
163 :
164 : template <typename T>
165 190413154 : void PetscVector<T>::add (const numeric_index_type i, const T value)
166 : {
167 190413154 : this->_restore_array();
168 16671969 : libmesh_assert_less (i, size());
169 :
170 190413154 : PetscInt i_val = static_cast<PetscInt>(i);
171 190413154 : PetscScalar petsc_value = PS(value);
172 :
173 190413154 : std::scoped_lock lock(this->_numeric_vector_mutex);
174 190413154 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES));
175 :
176 190413154 : this->_is_closed = false;
177 190413154 : }
178 :
179 :
180 :
181 : template <typename T>
182 74796616 : void PetscVector<T>::add_vector (const T * v,
183 : const std::vector<numeric_index_type> & dof_indices)
184 : {
185 : // If we aren't adding anything just return
186 74796616 : if (dof_indices.empty())
187 0 : return;
188 :
189 74796616 : this->_restore_array();
190 :
191 11356020 : const PetscInt * i_val = reinterpret_cast<const PetscInt *>(dof_indices.data());
192 5782458 : const PetscScalar * petsc_value = pPS(v);
193 :
194 74796616 : std::scoped_lock lock(this->_numeric_vector_mutex);
195 80370178 : LibmeshPetscCall(VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
196 : i_val, petsc_value, ADD_VALUES));
197 :
198 74796616 : this->_is_closed = false;
199 : }
200 :
201 :
202 :
203 : template <typename T>
204 2978796 : void PetscVector<T>::add_vector (const NumericVector<T> & v_in,
205 : const SparseMatrix<T> & A_in)
206 : {
207 85138 : parallel_object_only();
208 :
209 2978796 : this->_restore_array();
210 : // Make sure the data passed in are really of Petsc types
211 85138 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
212 85138 : const PetscMatrixBase<T> * A = cast_ptr<const PetscMatrixBase<T> *>(&A_in);
213 :
214 :
215 : // We shouldn't close() the matrix for you, as that would potentially modify the state of a const object.
216 2978796 : if (!A->closed())
217 : {
218 : libmesh_warning("Matrix A must be assembled before calling PetscVector::add_vector(v, A).\n"
219 : "Please update your code, as this warning will become an error in a future release.");
220 : libmesh_deprecated();
221 0 : const_cast<PetscMatrixBase<T> *>(A)->close();
222 : }
223 :
224 : // The const_cast<> is not elegant, but it is required since PETSc
225 : // expects a non-const Mat.
226 2978796 : LibmeshPetscCall(MatMultAdd(const_cast<PetscMatrixBase<T> *>(A)->mat(), v->_vec, _vec, _vec));
227 2978796 : }
228 :
229 :
230 :
231 : template <typename T>
232 0 : void PetscVector<T>::add_vector_transpose (const NumericVector<T> & v_in,
233 : const SparseMatrix<T> & A_in)
234 : {
235 0 : parallel_object_only();
236 :
237 0 : this->_restore_array();
238 : // Make sure the data passed in are really of Petsc types
239 0 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
240 0 : const PetscMatrixBase<T> * A = cast_ptr<const PetscMatrixBase<T> *>(&A_in);
241 :
242 :
243 : // We shouldn't close() the matrix for you, as that would potentially modify the state of a const object.
244 0 : if (!A->closed())
245 : {
246 : libmesh_warning("Matrix A must be assembled before calling PetscVector::add_vector_transpose(v, A).\n"
247 : "Please update your code, as this warning will become an error in a future release.");
248 : libmesh_deprecated();
249 0 : const_cast<PetscMatrixBase<T> *>(A)->close();
250 : }
251 :
252 : // The const_cast<> is not elegant, but it is required since PETSc
253 : // expects a non-const Mat.
254 0 : LibmeshPetscCall(MatMultTransposeAdd(const_cast<PetscMatrixBase<T> *>(A)->mat(), v->_vec, _vec, _vec));
255 0 : }
256 :
257 :
258 :
259 : template <typename T>
260 0 : void PetscVector<T>::add_vector_conjugate_transpose (const NumericVector<T> & v_in,
261 : const SparseMatrix<T> & A_in)
262 : {
263 0 : parallel_object_only();
264 :
265 0 : this->_restore_array();
266 : // Make sure the data passed in are really of Petsc types
267 0 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
268 0 : const PetscMatrixBase<T> * A = cast_ptr<const PetscMatrixBase<T> *>(&A_in);
269 :
270 : // We shouldn't close() the matrix for you, as that would potentially modify the state of a const object.
271 0 : if (!A->closed())
272 : {
273 : libmesh_warning("Matrix A must be assembled before calling PetscVector::add_vector_conjugate_transpose(v, A).\n"
274 : "Please update your code, as this warning will become an error in a future release.");
275 : libmesh_deprecated();
276 0 : const_cast<PetscMatrixBase<T> *>(A)->close();
277 : }
278 :
279 : // Store a temporary copy since MatMultHermitianTransposeAdd doesn't seem to work
280 : // TODO: Find out why MatMultHermitianTransposeAdd doesn't work, might be a PETSc bug?
281 0 : std::unique_ptr<NumericVector<Number>> this_clone = this->clone();
282 :
283 : // The const_cast<> is not elegant, but it is required since PETSc
284 : // expects a non-const Mat.
285 0 : LibmeshPetscCall(MatMultHermitianTranspose(const_cast<PetscMatrixBase<T> *>(A)->mat(), v->_vec, _vec));
286 :
287 : // Add the temporary copy to the matvec result
288 0 : this->add(1., *this_clone);
289 0 : }
290 :
291 :
292 :
293 : template <typename T>
294 346 : void PetscVector<T>::add (const T v_in)
295 : {
296 346 : this->_get_array(false);
297 :
298 17382 : for (numeric_index_type i=0; i<_local_size; i++)
299 17036 : _values[i] += PetscScalar(v_in);
300 346 : }
301 :
302 :
303 :
304 : template <typename T>
305 603238 : void PetscVector<T>::add (const NumericVector<T> & v)
306 : {
307 17288 : parallel_object_only();
308 :
309 603238 : this->add (1., v);
310 603238 : }
311 :
312 :
313 :
314 : template <typename T>
315 4369171 : void PetscVector<T>::add (const T a_in, const NumericVector<T> & v_in)
316 : {
317 124982 : parallel_object_only();
318 :
319 4369171 : this->_restore_array();
320 :
321 : // VecAXPY doesn't support &x==&y
322 4369171 : if (this == &v_in)
323 : {
324 140 : this->scale(a_in+1);
325 140 : return;
326 : }
327 :
328 124978 : PetscScalar a = PS(a_in);
329 :
330 : // Make sure the NumericVector passed in is really a PetscVector
331 124978 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
332 4369031 : v->_restore_array();
333 :
334 124978 : libmesh_assert_equal_to (this->size(), v->size());
335 :
336 4369031 : LibmeshPetscCall(VecAXPY(_vec, a, v->vec()));
337 :
338 124978 : libmesh_assert(this->comm().verify(
339 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
340 :
341 124978 : if (this->is_effectively_ghosted())
342 96048 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
343 :
344 4369031 : this->_is_closed = true;
345 : }
346 :
347 :
348 :
349 : template <typename T>
350 5109831 : void PetscVector<T>::insert (const T * v,
351 : const std::vector<numeric_index_type> & dof_indices)
352 : {
353 5109831 : if (dof_indices.empty())
354 0 : return;
355 :
356 5109831 : this->_restore_array();
357 :
358 923670 : PetscInt * idx_values = numeric_petsc_cast(dof_indices.data());
359 5109831 : std::scoped_lock lock(this->_numeric_vector_mutex);
360 5571666 : LibmeshPetscCall(VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
361 : idx_values, pPS(v), INSERT_VALUES));
362 :
363 5109831 : this->_is_closed = false;
364 : }
365 :
366 :
367 :
368 : template <typename T>
369 632902 : void PetscVector<T>::scale (const T factor_in)
370 : {
371 18076 : parallel_object_only();
372 :
373 632902 : this->_restore_array();
374 :
375 18076 : PetscScalar factor = PS(factor_in);
376 :
377 632902 : LibmeshPetscCall(VecScale(_vec, factor));
378 :
379 18076 : libmesh_assert(this->comm().verify(
380 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
381 :
382 18076 : if (this->is_effectively_ghosted())
383 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
384 632902 : }
385 :
386 : template <typename T>
387 140 : NumericVector<T> & PetscVector<T>::operator *= (const NumericVector<T> & v)
388 : {
389 4 : parallel_object_only();
390 :
391 :
392 4 : const PetscVector<T> * v_vec = cast_ptr<const PetscVector<T> *>(&v);
393 :
394 140 : LibmeshPetscCall(VecPointwiseMult(_vec, _vec, v_vec->_vec));
395 :
396 140 : return *this;
397 : }
398 :
399 : template <typename T>
400 79452 : NumericVector<T> & PetscVector<T>::operator /= (const NumericVector<T> & v)
401 : {
402 2448 : parallel_object_only();
403 :
404 :
405 2448 : const PetscVector<T> * v_vec = cast_ptr<const PetscVector<T> *>(&v);
406 :
407 79452 : LibmeshPetscCall(VecPointwiseDivide(_vec, _vec, v_vec->_vec));
408 :
409 79452 : return *this;
410 : }
411 :
412 : template <typename T>
413 0 : void PetscVector<T>::abs()
414 : {
415 0 : parallel_object_only();
416 :
417 0 : this->_restore_array();
418 :
419 0 : LibmeshPetscCall(VecAbs(_vec));
420 :
421 0 : libmesh_assert(this->comm().verify(
422 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
423 :
424 0 : if (this->is_effectively_ghosted())
425 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
426 0 : }
427 :
428 : template <typename T>
429 10363074 : T PetscVector<T>::dot (const NumericVector<T> & v_in) const
430 : {
431 296116 : parallel_object_only();
432 :
433 10363074 : this->_restore_array();
434 :
435 : // Error flag
436 :
437 : // Return value
438 10363074 : PetscScalar value=0.;
439 :
440 : // Make sure the NumericVector passed in is really a PetscVector
441 296116 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
442 :
443 : // 2.3.x (at least) style. Untested for previous versions.
444 10363074 : LibmeshPetscCall(VecDot(this->_vec, v->_vec, &value));
445 :
446 10363074 : return static_cast<T>(value);
447 : }
448 :
449 : template <typename T>
450 0 : T PetscVector<T>::indefinite_dot (const NumericVector<T> & v_in) const
451 : {
452 0 : parallel_object_only();
453 :
454 0 : this->_restore_array();
455 :
456 : // Error flag
457 :
458 : // Return value
459 0 : PetscScalar value=0.;
460 :
461 : // Make sure the NumericVector passed in is really a PetscVector
462 0 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
463 :
464 : // 2.3.x (at least) style. Untested for previous versions.
465 0 : LibmeshPetscCall(VecTDot(this->_vec, v->_vec, &value));
466 :
467 0 : return static_cast<T>(value);
468 : }
469 :
470 :
471 : template <typename T>
472 : NumericVector<T> &
473 140 : PetscVector<T>::operator = (const T s_in)
474 : {
475 4 : parallel_object_only();
476 :
477 140 : this->_restore_array();
478 4 : libmesh_assert(this->closed());
479 :
480 4 : PetscScalar s = PS(s_in);
481 :
482 140 : if (this->size() != 0)
483 : {
484 140 : LibmeshPetscCall(VecSet(_vec, s));
485 :
486 4 : libmesh_assert(this->comm().verify(
487 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
488 :
489 4 : if (this->is_effectively_ghosted())
490 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
491 : }
492 :
493 140 : return *this;
494 : }
495 :
496 :
497 :
498 : template <typename T>
499 : NumericVector<T> &
500 3940565 : PetscVector<T>::operator = (const NumericVector<T> & v_in)
501 : {
502 110384 : parallel_object_only();
503 :
504 : // Make sure the NumericVector passed in is really a PetscVector
505 110384 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
506 :
507 3940565 : *this = *v;
508 :
509 3940565 : return *this;
510 : }
511 :
512 :
513 :
514 : template <typename T>
515 : PetscVector<T> &
516 3940635 : PetscVector<T>::operator = (const PetscVector<T> & v)
517 : {
518 : enum AssignmentType { ParallelToSerial, SerialToParallel, SameToSame, Error };
519 :
520 110386 : parallel_object_only();
521 110386 : libmesh_assert(this->comm().verify(
522 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
523 110386 : libmesh_assert(this->comm().verify(
524 : static_cast<typename std::underlying_type<ParallelType>::type>(v.type())));
525 :
526 3940635 : if (this == &v)
527 480 : return *this;
528 :
529 3923835 : this->_restore_array();
530 3923835 : v._restore_array();
531 :
532 109906 : libmesh_assert_equal_to (this->size(), v.size());
533 109906 : libmesh_assert (v.closed());
534 :
535 109906 : AssignmentType assign_type = Error;
536 109906 : if (this->is_effectively_serial() && !v.is_effectively_serial())
537 48 : assign_type = ParallelToSerial;
538 109858 : else if (!this->is_effectively_serial() && v.is_effectively_serial())
539 0 : assign_type = SerialToParallel;
540 3922179 : else if (this->local_size() == v.local_size())
541 109858 : assign_type = SameToSame;
542 :
543 109906 : libmesh_assert(this->comm().verify(
544 : static_cast<typename std::underlying_type<AssignmentType>::type>(assign_type)));
545 :
546 109906 : switch (assign_type)
547 : {
548 48 : case ParallelToSerial:
549 : {
550 : // scatter from parallel to serial
551 48 : libmesh_assert(v.comm().size() > 1);
552 96 : WrappedPetsc<VecScatter> scatter;
553 1656 : LibmeshPetscCall(VecScatterCreateToAll(v._vec, scatter.get(), nullptr));
554 1656 : VecScatterBeginEnd(v.comm(), scatter, v._vec, _vec, INSERT_VALUES, SCATTER_FORWARD);
555 48 : break;
556 : }
557 109858 : case SameToSame:
558 : // serial to serial or parallel to parallel
559 3922179 : LibmeshPetscCall(VecCopy(v._vec, _vec));
560 109858 : break;
561 0 : case SerialToParallel:
562 0 : libmesh_not_implemented_msg(
563 : "Scattering from a serial vector on every rank to a parallel vector is not behavior we "
564 : "define because we do not verify the serial vector is the same on each rank");
565 0 : default:
566 0 : libmesh_error_msg("Unhandled vector combination");
567 : }
568 :
569 109906 : libmesh_assert(this->comm().verify(
570 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
571 :
572 109906 : if (this->is_effectively_ghosted())
573 2867104 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
574 :
575 3923835 : this->_is_closed = true;
576 :
577 3923835 : return *this;
578 : }
579 :
580 :
581 :
582 : template <typename T>
583 : NumericVector<T> &
584 70 : PetscVector<T>::operator = (const std::vector<T> & v)
585 : {
586 2 : parallel_object_only();
587 :
588 70 : this->_get_array(false);
589 :
590 : /**
591 : * Case 1: The vector is the same size of
592 : * The global vector. Only add the local components.
593 : */
594 70 : if (this->size() == v.size())
595 : {
596 70 : numeric_index_type first = first_local_index();
597 70 : numeric_index_type last = last_local_index();
598 300860 : for (numeric_index_type i=0; i<last-first; i++)
599 309384 : _values[i] = PS(v[first + i]);
600 : }
601 :
602 : /**
603 : * Case 2: The vector is the same size as our local
604 : * piece. Insert directly to the local piece.
605 : */
606 : else
607 : {
608 0 : for (numeric_index_type i=0; i<_local_size; i++)
609 0 : _values[i] = PS(v[i]);
610 : }
611 :
612 : // Make sure ghost dofs are up to date
613 2 : if (this->is_effectively_ghosted())
614 0 : this->close();
615 :
616 70 : return *this;
617 : }
618 :
619 :
620 :
621 : template <typename T>
622 126761 : void PetscVector<T>::localize (NumericVector<T> & v_local_in) const
623 : {
624 62864 : parallel_object_only();
625 :
626 62864 : libmesh_assert(this->comm().verify(
627 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
628 62864 : libmesh_assert(this->comm().verify(
629 : static_cast<typename std::underlying_type<ParallelType>::type>(v_local_in.type())));
630 :
631 126761 : v_local_in = *this;
632 126761 : }
633 :
634 :
635 :
636 : template <typename T>
637 2274133 : void PetscVector<T>::localize (NumericVector<T> & v_local_in,
638 : const std::vector<numeric_index_type> & /*send_list*/) const
639 : {
640 124624 : this->localize(v_local_in);
641 140464 : }
642 :
643 :
644 :
645 : template <typename T>
646 140 : void PetscVector<T>::localize (std::vector<T> & v_local,
647 : const std::vector<numeric_index_type> & indices) const
648 : {
649 4 : parallel_object_only();
650 :
651 : // Error code used to check the status of all PETSc function calls.
652 :
653 : // Create a sequential destination Vec with the right number of entries on each proc.
654 8 : WrappedPetsc<Vec> dest;
655 144 : LibmeshPetscCall(VecCreateSeq(PETSC_COMM_SELF, cast_int<PetscInt>(indices.size()), dest.get()));
656 :
657 : // Create an IS using the libmesh routine. PETSc does not own the
658 : // IS memory in this case, it is automatically cleaned up by the
659 : // std::vector destructor.
660 8 : PetscInt * idxptr =
661 140 : indices.empty() ? nullptr : numeric_petsc_cast(indices.data());
662 8 : WrappedPetsc<IS> is;
663 140 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(), cast_int<PetscInt>(indices.size()), idxptr,
664 : PETSC_USE_POINTER, is.get()));
665 :
666 : // Create the VecScatter object. "PETSC_NULL" means "use the identity IS".
667 8 : WrappedPetsc<VecScatter> scatter;
668 140 : LibmeshPetscCall(VecScatterCreate(_vec,
669 : /*src is=*/is,
670 : /*dest vec=*/dest,
671 : /*dest is=*/LIBMESH_PETSC_NULLPTR,
672 : scatter.get()));
673 :
674 : // Do the scatter
675 140 : VecScatterBeginEnd(this->comm(), scatter, _vec, dest, INSERT_VALUES, SCATTER_FORWARD);
676 :
677 : // Get access to the values stored in dest.
678 : PetscScalar * values;
679 140 : LibmeshPetscCall(VecGetArray (dest, &values));
680 :
681 : // Store values into the provided v_local. Make sure there is enough
682 : // space reserved and then clear out any existing entries before
683 : // inserting.
684 144 : v_local.reserve(indices.size());
685 4 : v_local.clear();
686 144 : v_local.insert(v_local.begin(), values, values+indices.size());
687 :
688 : // We are done using it, so restore the array.
689 140 : LibmeshPetscCall(VecRestoreArray (dest, &values));
690 140 : }
691 :
692 :
693 :
694 : template <typename T>
695 17428 : void PetscVector<T>::localize (const numeric_index_type first_local_idx,
696 : const numeric_index_type last_local_idx,
697 : const std::vector<numeric_index_type> & send_list)
698 : {
699 672 : parallel_object_only();
700 :
701 17428 : this->_restore_array();
702 :
703 672 : libmesh_assert_less_equal (send_list.size(), this->size());
704 672 : libmesh_assert_less_equal (last_local_idx+1, this->size());
705 :
706 17428 : const numeric_index_type my_size = this->size();
707 17428 : const numeric_index_type my_local_size = (last_local_idx + 1 - first_local_idx);
708 :
709 : // Don't bother for serial cases
710 : // if ((first_local_idx == 0) &&
711 : // (my_local_size == my_size))
712 : // But we do need to stay in sync for degenerate cases
713 18100 : if (this->n_processors() == 1)
714 244 : return;
715 :
716 :
717 : // Build a parallel vector, initialize it with the local
718 : // parts of (*this)
719 18528 : PetscVector<T> parallel_vec(this->comm(), PARALLEL);
720 :
721 17184 : parallel_vec.init (my_size, my_local_size, true, PARALLEL);
722 :
723 :
724 : // Copy part of *this into the parallel_vec
725 : {
726 : // Create idx, idx[i] = i+first_local_idx;
727 17856 : std::vector<PetscInt> idx(my_local_size);
728 672 : std::iota (idx.begin(), idx.end(), first_local_idx);
729 :
730 : // Create the index set & scatter objects
731 1344 : WrappedPetsc<IS> is;
732 32470 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(), my_local_size,
733 : my_local_size ? idx.data() : nullptr, PETSC_USE_POINTER, is.get()));
734 :
735 672 : WrappedPetsc<VecScatter> scatter;
736 17184 : LibmeshPetscCall(VecScatterCreate(_vec, is,
737 : parallel_vec._vec, is,
738 : scatter.get()));
739 :
740 : // Perform the scatter
741 17184 : VecScatterBeginEnd(this->comm(), scatter, _vec, parallel_vec._vec, INSERT_VALUES, SCATTER_FORWARD);
742 : }
743 :
744 : // localize like normal
745 17184 : parallel_vec.close();
746 17184 : parallel_vec.localize (*this, send_list);
747 17184 : this->close();
748 15840 : }
749 :
750 :
751 :
752 : template <typename T>
753 21350 : void PetscVector<T>::localize (std::vector<T> & v_local) const
754 : {
755 610 : parallel_object_only();
756 :
757 21350 : this->_restore_array();
758 :
759 : // This function must be run on all processors at once
760 610 : parallel_object_only();
761 :
762 21350 : const PetscInt n = this->size();
763 21350 : const PetscInt nl = this->local_size();
764 : PetscScalar * values;
765 :
766 610 : v_local.clear();
767 21350 : v_local.resize(n, 0.);
768 :
769 21350 : LibmeshPetscCall(VecGetArray (_vec, &values));
770 :
771 21350 : numeric_index_type ioff = first_local_index();
772 :
773 13196761 : for (PetscInt i=0; i<nl; i++)
774 14372892 : v_local[i+ioff] = static_cast<T>(values[i]);
775 :
776 21350 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
777 :
778 21350 : this->comm().sum(v_local);
779 21350 : }
780 :
781 :
782 :
783 : // Full specialization for Real datatypes
784 : #ifdef LIBMESH_USE_REAL_NUMBERS
785 :
786 : template <>
787 85454 : void PetscVector<Real>::localize_to_one (std::vector<Real> & v_local,
788 : const processor_id_type
789 : timpi_mpi_var(pid)) const
790 : {
791 2694 : parallel_object_only();
792 :
793 85454 : this->_restore_array();
794 :
795 85454 : const PetscInt n = size();
796 : PetscScalar * values;
797 :
798 : // only one processor
799 88148 : if (n_processors() == 1)
800 : {
801 1213 : v_local.resize(n);
802 :
803 1213 : LibmeshPetscCall(VecGetArray (_vec, &values));
804 :
805 2717154 : for (PetscInt i=0; i<n; i++)
806 2715941 : v_local[i] = static_cast<Real>(values[i]);
807 :
808 1213 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
809 : }
810 :
811 : // otherwise multiple processors
812 : #ifdef LIBMESH_HAVE_MPI
813 : else
814 : {
815 84241 : if (pid == 0) // optimized version for localizing to 0
816 : {
817 5388 : WrappedPetsc<Vec> vout;
818 5388 : WrappedPetsc<VecScatter> ctx;
819 :
820 84241 : LibmeshPetscCall(VecScatterCreateToZero(_vec, ctx.get(), vout.get()));
821 :
822 84241 : VecScatterBeginEnd(this->comm(), ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
823 :
824 86935 : if (processor_id() == 0)
825 : {
826 12451 : v_local.resize(n);
827 :
828 12451 : LibmeshPetscCall(VecGetArray (vout, &values));
829 :
830 27244108 : for (PetscInt i=0; i<n; i++)
831 29922653 : v_local[i] = static_cast<Real>(values[i]);
832 :
833 12451 : LibmeshPetscCall(VecRestoreArray (vout, &values));
834 : }
835 : }
836 : else
837 : {
838 0 : v_local.resize(n);
839 :
840 0 : numeric_index_type ioff = this->first_local_index();
841 0 : std::vector<Real> local_values (n, 0.);
842 :
843 : {
844 0 : LibmeshPetscCall(VecGetArray (_vec, &values));
845 :
846 0 : const PetscInt nl = local_size();
847 0 : for (PetscInt i=0; i<nl; i++)
848 0 : local_values[i+ioff] = static_cast<Real>(values[i]);
849 :
850 0 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
851 : }
852 :
853 :
854 0 : MPI_Reduce (local_values.data(), v_local.data(), n,
855 0 : Parallel::StandardType<Real>(),
856 : Parallel::OpFunction<Real>::sum(), pid,
857 0 : this->comm().get());
858 : }
859 : }
860 : #endif // LIBMESH_HAVE_MPI
861 85454 : }
862 :
863 : #endif
864 :
865 :
866 : // Full specialization for Complex datatypes
867 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
868 :
869 : template <>
870 : void PetscVector<Complex>::localize_to_one (std::vector<Complex> & v_local,
871 : const processor_id_type pid) const
872 : {
873 : parallel_object_only();
874 :
875 : this->_restore_array();
876 :
877 : const PetscInt n = size();
878 : const PetscInt nl = local_size();
879 : PetscScalar * values;
880 :
881 :
882 : v_local.resize(n);
883 :
884 :
885 : for (PetscInt i=0; i<n; i++)
886 : v_local[i] = 0.;
887 :
888 : // only one processor
889 : if (n_processors() == 1)
890 : {
891 : LibmeshPetscCall(VecGetArray (_vec, &values));
892 :
893 : for (PetscInt i=0; i<n; i++)
894 : v_local[i] = static_cast<Complex>(values[i]);
895 :
896 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
897 : }
898 :
899 : // otherwise multiple processors
900 : else
901 : {
902 : numeric_index_type ioff = this->first_local_index();
903 :
904 : /* in here the local values are stored, acting as send buffer for MPI
905 : * initialize to zero, since we collect using MPI_SUM
906 : */
907 : std::vector<Real> real_local_values(n, 0.);
908 : std::vector<Real> imag_local_values(n, 0.);
909 :
910 : {
911 : LibmeshPetscCall(VecGetArray (_vec, &values));
912 :
913 : // provide my local share to the real and imag buffers
914 : for (PetscInt i=0; i<nl; i++)
915 : {
916 : real_local_values[i+ioff] = static_cast<Complex>(values[i]).real();
917 : imag_local_values[i+ioff] = static_cast<Complex>(values[i]).imag();
918 : }
919 :
920 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
921 : }
922 :
923 : /* have buffers of the real and imaginary part of v_local.
924 : * Once MPI_Reduce() collected all the real and imaginary
925 : * parts in these std::vector<Real>, the values can be
926 : * copied to v_local
927 : */
928 : std::vector<Real> real_v_local(n);
929 : std::vector<Real> imag_v_local(n);
930 :
931 : // collect entries from other proc's in real_v_local, imag_v_local
932 : MPI_Reduce (real_local_values.data(),
933 : real_v_local.data(), n,
934 : Parallel::StandardType<Real>(),
935 : Parallel::OpFunction<Real>::sum(),
936 : pid, this->comm().get());
937 :
938 : MPI_Reduce (imag_local_values.data(),
939 : imag_v_local.data(), n,
940 : Parallel::StandardType<Real>(),
941 : Parallel::OpFunction<Real>::sum(),
942 : pid, this->comm().get());
943 :
944 : // copy real_v_local and imag_v_local to v_local
945 : for (PetscInt i=0; i<n; i++)
946 : v_local[i] = Complex(real_v_local[i], imag_v_local[i]);
947 : }
948 : }
949 :
950 : #endif
951 :
952 :
953 :
954 : template <typename T>
955 70 : void PetscVector<T>::pointwise_mult (const NumericVector<T> & vec1,
956 : const NumericVector<T> & vec2)
957 : {
958 2 : parallel_object_only();
959 :
960 70 : this->_restore_array();
961 :
962 : // Convert arguments to PetscVector*.
963 2 : const PetscVector<T> * vec1_petsc = cast_ptr<const PetscVector<T> *>(&vec1);
964 2 : const PetscVector<T> * vec2_petsc = cast_ptr<const PetscVector<T> *>(&vec2);
965 :
966 : // Call PETSc function.
967 70 : LibmeshPetscCall(VecPointwiseMult(_vec, vec1_petsc->vec(), vec2_petsc->vec()));
968 :
969 2 : libmesh_assert(this->comm().verify(
970 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
971 :
972 2 : if (this->is_effectively_ghosted())
973 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
974 :
975 70 : this->_is_closed = true;
976 70 : }
977 :
978 : template <typename T>
979 70 : void PetscVector<T>::pointwise_divide (const NumericVector<T> & vec1,
980 : const NumericVector<T> & vec2)
981 : {
982 2 : parallel_object_only();
983 :
984 70 : this->_restore_array();
985 :
986 : // Convert arguments to PetscVector*.
987 2 : const PetscVector<T> * const vec1_petsc = cast_ptr<const PetscVector<T> *>(&vec1);
988 2 : const PetscVector<T> * const vec2_petsc = cast_ptr<const PetscVector<T> *>(&vec2);
989 :
990 : // Call PETSc function.
991 70 : LibmeshPetscCall(VecPointwiseDivide(_vec, vec1_petsc->vec(), vec2_petsc->vec()));
992 :
993 2 : libmesh_assert(this->comm().verify(
994 : static_cast<typename std::underlying_type<ParallelType>::type>(this->type())));
995 :
996 2 : if (this->is_effectively_ghosted())
997 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
998 :
999 70 : this->_is_closed = true;
1000 70 : }
1001 :
1002 : template <typename T>
1003 140 : void PetscVector<T>::print_matlab (const std::string & name) const
1004 : {
1005 4 : parallel_object_only();
1006 :
1007 140 : this->_restore_array();
1008 4 : libmesh_assert (this->closed());
1009 :
1010 : // Create an ASCII file containing the matrix
1011 : // if a filename was provided.
1012 140 : if (name != "")
1013 : {
1014 8 : WrappedPetsc<PetscViewer> petsc_viewer;
1015 :
1016 140 : LibmeshPetscCall(PetscViewerASCIIOpen(this->comm().get(),
1017 : name.c_str(),
1018 : petsc_viewer.get()));
1019 :
1020 : #if PETSC_VERSION_LESS_THAN(3,7,0)
1021 : LibmeshPetscCall(PetscViewerSetFormat (petsc_viewer,
1022 : PETSC_VIEWER_ASCII_MATLAB));
1023 : #else
1024 140 : LibmeshPetscCall(PetscViewerPushFormat (petsc_viewer,
1025 : PETSC_VIEWER_ASCII_MATLAB));
1026 : #endif
1027 :
1028 140 : LibmeshPetscCall(VecView (_vec, petsc_viewer));
1029 : }
1030 :
1031 : // Otherwise the matrix will be dumped to the screen.
1032 : else
1033 : {
1034 :
1035 : #if PETSC_VERSION_LESS_THAN(3,7,0)
1036 : LibmeshPetscCall(PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
1037 : PETSC_VIEWER_ASCII_MATLAB));
1038 : #else
1039 0 : LibmeshPetscCall(PetscViewerPushFormat (PETSC_VIEWER_STDOUT_WORLD,
1040 : PETSC_VIEWER_ASCII_MATLAB));
1041 : #endif
1042 :
1043 0 : LibmeshPetscCall(VecView (_vec, PETSC_VIEWER_STDOUT_WORLD));
1044 : }
1045 140 : }
1046 :
1047 :
1048 :
1049 :
1050 :
1051 : template <typename T>
1052 8610 : void PetscVector<T>::create_subvector(NumericVector<T> & subvector,
1053 : const std::vector<numeric_index_type> & rows,
1054 : const bool supplying_global_rows) const
1055 : {
1056 246 : parallel_object_only();
1057 :
1058 246 : libmesh_error_msg_if(
1059 : subvector.is_effectively_ghosted(),
1060 : "We do not support scattering parallel information to ghosts for subvectors");
1061 :
1062 8610 : this->_restore_array();
1063 :
1064 : // PETSc data structures
1065 :
1066 : // Make sure the passed in subvector is really a PetscVector
1067 246 : PetscVector<T> * petsc_subvector = cast_ptr<PetscVector<T> *>(&subvector);
1068 :
1069 : // If the petsc_subvector is already initialized, we assume that the
1070 : // user has already allocated the *correct* amount of space for it.
1071 : // If not, we use the appropriate PETSc routines to initialize it.
1072 8610 : if (!petsc_subvector->initialized())
1073 : {
1074 0 : if (this->is_effectively_serial())
1075 : {
1076 0 : libmesh_assert(this->comm().verify(rows.size()));
1077 0 : LibmeshPetscCall(VecCreateSeq(this->comm().get(), rows.size(), &(petsc_subvector->_vec)));
1078 :
1079 : #ifdef LIBMESH_ENABLE_DEPRECATED
1080 : // In the past we always created PARALLEL subvectors from
1081 : // GHOSTED vectors, but this behavior must now be specified
1082 : // downstream when it's desired, by setting the subvector
1083 : // type in advance.
1084 0 : petsc_subvector->_type = (this->type() == SERIAL) ? SERIAL : PARALLEL;
1085 :
1086 0 : if (petsc_subvector->type() == AUTOMATIC &&
1087 0 : this->type() != petsc_subvector->type())
1088 : libmesh_deprecated();
1089 : #else
1090 : // If the subvector already had a type, let's not change it.
1091 : // If it was just at the default AUTOMATIC, let's pick the
1092 : // true type from the supervector.
1093 : if (petsc_subvector->type() == AUTOMATIC)
1094 : petsc_subvector->_type = this->type();
1095 : #endif
1096 : }
1097 : else
1098 : {
1099 : #ifndef LIBMESH_ENABLE_DEPRECATED
1100 : // We're only supporting PARALLEL subvectors in parallel
1101 : // so far, but we want to leave other possible API inputs
1102 : // marked as invalid so we can add support for other
1103 : // subvector types later.
1104 : if (petsc_subvector->type() != PARALLEL &&
1105 : (petsc_subvector->type() != AUTOMATIC ||
1106 : this->type() != PARALLEL))
1107 : libmesh_not_implemented_msg("Cannot yet create non-PARALLEL subvectors in parallel");
1108 : #endif
1109 :
1110 0 : if (supplying_global_rows)
1111 : // Initialize the petsc_subvector to have enough space to hold
1112 : // the entries which will be scattered into it. Note: such an
1113 : // init() function (where we let PETSc decide the number of local
1114 : // entries) is not currently offered by the PetscVector
1115 : // class. Should we differentiate here between sequential and
1116 : // parallel vector creation based on this->n_processors() ?
1117 0 : LibmeshPetscCall(VecCreateMPI(this->comm().get(),
1118 : PETSC_DECIDE, // n_local
1119 : cast_int<PetscInt>(rows.size()), // n_global
1120 : &(petsc_subvector->_vec)));
1121 : else
1122 0 : LibmeshPetscCall(VecCreateMPI(this->comm().get(),
1123 : cast_int<PetscInt>(rows.size()),
1124 : PETSC_DETERMINE,
1125 : &(petsc_subvector->_vec)));
1126 :
1127 : // We created a parallel vector
1128 0 : petsc_subvector->_type = PARALLEL;
1129 : }
1130 :
1131 0 : LibmeshPetscCall(VecSetFromOptions (petsc_subvector->_vec));
1132 :
1133 :
1134 : // Mark the subvector as initialized
1135 0 : petsc_subvector->_is_initialized = true;
1136 : }
1137 : else
1138 : {
1139 8610 : petsc_subvector->_restore_array();
1140 : }
1141 :
1142 9102 : std::vector<PetscInt> idx(rows.size());
1143 8610 : if (supplying_global_rows)
1144 4 : std::iota (idx.begin(), idx.end(), 0);
1145 : else
1146 : {
1147 : PetscInt start;
1148 8470 : LibmeshPetscCall(VecGetOwnershipRange(petsc_subvector->_vec, &start, nullptr));
1149 8470 : std::iota (idx.begin(), idx.end(), start);
1150 : }
1151 :
1152 : // Construct index sets
1153 492 : WrappedPetsc<IS> parent_is;
1154 8856 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1155 : cast_int<PetscInt>(rows.size()),
1156 : numeric_petsc_cast(rows.data()),
1157 : PETSC_USE_POINTER,
1158 : parent_is.get()));
1159 :
1160 492 : WrappedPetsc<IS> subvector_is;
1161 8856 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1162 : cast_int<PetscInt>(rows.size()),
1163 : idx.data(),
1164 : PETSC_USE_POINTER,
1165 : subvector_is.get()));
1166 :
1167 : // Construct the scatter object
1168 246 : WrappedPetsc<VecScatter> scatter;
1169 8610 : LibmeshPetscCall(VecScatterCreate(this->_vec,
1170 : parent_is,
1171 : petsc_subvector->_vec,
1172 : subvector_is,
1173 : scatter.get()));
1174 :
1175 : // Actually perform the scatter
1176 8610 : VecScatterBeginEnd(this->comm(), scatter, this->_vec, petsc_subvector->_vec, INSERT_VALUES, SCATTER_FORWARD);
1177 :
1178 8610 : petsc_subvector->_is_closed = true;
1179 8610 : }
1180 :
1181 :
1182 :
1183 : template <typename T>
1184 7454797555 : void PetscVector<T>::_get_array(bool read_only) const
1185 : {
1186 589768282 : libmesh_assert (this->initialized());
1187 :
1188 589768282 : bool initially_array_is_present = _array_is_present.load(std::memory_order_acquire);
1189 :
1190 : // If we already have a read/write array - and we're trying
1191 : // to get a read only array - let's just use the read write
1192 7454797555 : if (initially_array_is_present && read_only && !_values_read_only)
1193 24807602 : _read_only_values = _values;
1194 :
1195 : // If the values have already been retrieved and we're currently
1196 : // trying to get a non-read only view (ie read/write) and the
1197 : // values are currently read only... then we need to restore
1198 : // the array first... and then retrieve it again.
1199 7454797555 : if (initially_array_is_present && !read_only && _values_read_only)
1200 : {
1201 210 : _restore_array();
1202 6 : initially_array_is_present = false;
1203 : }
1204 :
1205 7454797351 : if (!initially_array_is_present)
1206 : {
1207 1956665 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1208 1841419 : if (!_array_is_present.load(std::memory_order_relaxed))
1209 : {
1210 115099 : if (!this->is_effectively_ghosted())
1211 : {
1212 124789 : if (read_only)
1213 : {
1214 124163 : LibmeshPetscCall(VecGetArrayRead(_vec, &_read_only_values));
1215 124163 : _values_read_only = true;
1216 : }
1217 : else
1218 : {
1219 626 : LibmeshPetscCall(VecGetArray(_vec, &_values));
1220 626 : _values_read_only = false;
1221 : }
1222 124789 : _local_size = this->local_size();
1223 : }
1224 : else
1225 : {
1226 1716247 : LibmeshPetscCall(VecGhostGetLocalForm (_vec,&_local_form));
1227 :
1228 1716247 : if (read_only)
1229 : {
1230 1716247 : LibmeshPetscCall(VecGetArrayRead(_local_form, &_read_only_values));
1231 1716247 : _values_read_only = true;
1232 : }
1233 : else
1234 : {
1235 0 : LibmeshPetscCall(VecGetArray(_local_form, &_values));
1236 0 : _values_read_only = false;
1237 : }
1238 :
1239 1716247 : PetscInt my_local_size = 0;
1240 1716247 : LibmeshPetscCall(VecGetLocalSize(_local_form, &my_local_size));
1241 1716247 : _local_size = static_cast<numeric_index_type>(my_local_size);
1242 : }
1243 :
1244 : { // cache ownership range
1245 1841036 : PetscInt petsc_first=0, petsc_last=0;
1246 1841036 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1247 1841036 : _first = static_cast<numeric_index_type>(petsc_first);
1248 1841036 : _last = static_cast<numeric_index_type>(petsc_last);
1249 : }
1250 115099 : _array_is_present.store(true, std::memory_order_release);
1251 : }
1252 : }
1253 7454797555 : }
1254 :
1255 :
1256 :
1257 : template <typename T>
1258 360100611 : void PetscVector<T>::_restore_array() const
1259 : {
1260 360100611 : libmesh_error_msg_if(_values_manually_retrieved,
1261 : "PetscVector values were manually retrieved but are being automatically restored!");
1262 :
1263 28260730 : libmesh_assert (this->initialized());
1264 360100611 : if (_array_is_present.load(std::memory_order_acquire))
1265 : {
1266 1956135 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1267 1841036 : if (_array_is_present.load(std::memory_order_relaxed))
1268 : {
1269 115099 : if (!this->is_effectively_ghosted())
1270 : {
1271 124789 : if (_values_read_only)
1272 124163 : LibmeshPetscCall(VecRestoreArrayRead (_vec, &_read_only_values));
1273 : else
1274 626 : LibmeshPetscCall(VecRestoreArray (_vec, &_values));
1275 :
1276 124789 : _values = nullptr;
1277 : }
1278 : else
1279 : {
1280 1716247 : if (_values_read_only)
1281 1716247 : LibmeshPetscCall(VecRestoreArrayRead (_local_form, &_read_only_values));
1282 : else
1283 0 : LibmeshPetscCall(VecRestoreArray (_local_form, &_values));
1284 :
1285 1716247 : _values = nullptr;
1286 1716247 : LibmeshPetscCall(VecGhostRestoreLocalForm (_vec,&_local_form));
1287 1716247 : _local_form = nullptr;
1288 1716247 : _local_size = 0;
1289 : }
1290 115099 : _array_is_present.store(false, std::memory_order_release);
1291 : }
1292 : }
1293 360100611 : }
1294 :
1295 : template <typename T>
1296 : std::unique_ptr<NumericVector<T>>
1297 15730 : PetscVector<T>::get_subvector(const std::vector<numeric_index_type> & rows)
1298 : {
1299 : // Construct index set
1300 242 : WrappedPetsc<IS> parent_is;
1301 15972 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1302 : cast_int<PetscInt>(rows.size()),
1303 : numeric_petsc_cast(rows.data()),
1304 : PETSC_USE_POINTER,
1305 : parent_is.get()));
1306 :
1307 : Vec subvec;
1308 15730 : LibmeshPetscCall(VecGetSubVector(_vec, parent_is, &subvec));
1309 :
1310 15730 : this->_is_closed = false;
1311 :
1312 15972 : return std::make_unique<PetscVector<T>>(subvec, this->comm());
1313 : }
1314 :
1315 : template <typename T>
1316 : void
1317 15730 : PetscVector<T>::restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
1318 : const std::vector<numeric_index_type> & rows)
1319 : {
1320 242 : auto * const petsc_subvector = cast_ptr<PetscVector<T> *>(subvector.get());
1321 :
1322 : // Construct index set
1323 242 : WrappedPetsc<IS> parent_is;
1324 15972 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1325 : cast_int<PetscInt>(rows.size()),
1326 : numeric_petsc_cast(rows.data()),
1327 : PETSC_USE_POINTER,
1328 : parent_is.get()));
1329 :
1330 15730 : Vec subvec = petsc_subvector->vec();
1331 15730 : LibmeshPetscCall(VecRestoreSubVector(_vec, parent_is, &subvec));
1332 :
1333 242 : if (this->is_effectively_ghosted())
1334 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
1335 :
1336 15730 : this->_is_closed = true;
1337 15730 : }
1338 :
1339 : //------------------------------------------------------------------
1340 : // Explicit instantiations
1341 : template class LIBMESH_EXPORT PetscVector<Number>;
1342 :
1343 : } // namespace libMesh
1344 :
1345 :
1346 :
1347 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
