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 12 : T PetscVector<T>::sum () const
45 : {
46 12 : this->_restore_array();
47 4 : libmesh_assert(this->closed());
48 :
49 12 : PetscScalar value=0.;
50 :
51 12 : LibmeshPetscCall(VecSum (_vec, &value));
52 :
53 12 : return static_cast<T>(value);
54 : }
55 :
56 :
57 :
58 : template <typename T>
59 : template <NormType N>
60 22872 : Real PetscVector<T>::norm () const
61 : {
62 7654 : parallel_object_only();
63 :
64 22872 : this->_restore_array();
65 7654 : libmesh_assert(this->closed());
66 :
67 22872 : PetscReal value=0.;
68 :
69 22872 : LibmeshPetscCall(VecNorm (_vec, N, &value));
70 :
71 22872 : return static_cast<Real>(value);
72 : }
73 : template <typename T>
74 132 : Real PetscVector<T>::l1_norm () const
75 : {
76 132 : return PetscVector<T>::norm<NORM_1>();
77 : }
78 : template <typename T>
79 22722 : Real PetscVector<T>::l2_norm () const
80 : {
81 22722 : return PetscVector<T>::norm<NORM_2>();
82 : }
83 : template <typename T>
84 18 : Real PetscVector<T>::linfty_norm () const
85 : {
86 18 : return PetscVector<T>::norm<NORM_INFINITY>();
87 : }
88 :
89 :
90 :
91 : template <typename T>
92 : NumericVector<T> &
93 2814 : PetscVector<T>::operator += (const NumericVector<T> & v)
94 : {
95 938 : parallel_object_only();
96 :
97 2814 : this->_restore_array();
98 938 : libmesh_assert(this->closed());
99 :
100 2814 : this->add(1., v);
101 :
102 2814 : return *this;
103 : }
104 :
105 :
106 :
107 : template <typename T>
108 : NumericVector<T> &
109 6924 : PetscVector<T>::operator -= (const NumericVector<T> & v)
110 : {
111 2308 : parallel_object_only();
112 :
113 6924 : this->_restore_array();
114 2308 : libmesh_assert(this->closed());
115 :
116 6924 : this->add(-1., v);
117 :
118 6924 : return *this;
119 : }
120 :
121 :
122 :
123 : template <typename T>
124 11987710 : void PetscVector<T>::set (const numeric_index_type i, const T value)
125 : {
126 11987710 : this->_restore_array();
127 4052344 : libmesh_assert_less (i, size());
128 :
129 11987710 : PetscInt i_val = static_cast<PetscInt>(i);
130 11987710 : PetscScalar petsc_value = PS(value);
131 :
132 11987710 : std::scoped_lock lock(this->_numeric_vector_mutex);
133 11987710 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES));
134 :
135 11987710 : this->_is_closed = false;
136 11987710 : }
137 :
138 :
139 :
140 : template <typename T>
141 12 : 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 12 : LibmeshPetscCall(VecReciprocal(_vec));
148 12 : }
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 51081004 : void PetscVector<T>::add (const numeric_index_type i, const T value)
166 : {
167 51081004 : this->_restore_array();
168 16671969 : libmesh_assert_less (i, size());
169 :
170 51081004 : PetscInt i_val = static_cast<PetscInt>(i);
171 51081004 : PetscScalar petsc_value = PS(value);
172 :
173 51081004 : std::scoped_lock lock(this->_numeric_vector_mutex);
174 51081004 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES));
175 :
176 51081004 : this->_is_closed = false;
177 51081004 : }
178 :
179 :
180 :
181 : template <typename T>
182 18643977 : 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 18643977 : if (dof_indices.empty())
187 0 : return;
188 :
189 18643977 : 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 18643977 : std::scoped_lock lock(this->_numeric_vector_mutex);
195 24217539 : LibmeshPetscCall(VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
196 : i_val, petsc_value, ADD_VALUES));
197 :
198 18643977 : this->_is_closed = false;
199 : }
200 :
201 :
202 :
203 : template <typename T>
204 255414 : void PetscVector<T>::add_vector (const NumericVector<T> & v_in,
205 : const SparseMatrix<T> & A_in)
206 : {
207 85138 : parallel_object_only();
208 :
209 255414 : 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 255414 : 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 255414 : LibmeshPetscCall(MatMultAdd(const_cast<PetscMatrixBase<T> *>(A)->mat(), v->_vec, _vec, _vec));
227 255414 : }
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 26 : void PetscVector<T>::add (const T v_in)
295 : {
296 26 : this->_get_array(false);
297 :
298 2718 : for (numeric_index_type i=0; i<_local_size; i++)
299 2692 : _values[i] += PetscScalar(v_in);
300 26 : }
301 :
302 :
303 :
304 : template <typename T>
305 51864 : void PetscVector<T>::add (const NumericVector<T> & v)
306 : {
307 17288 : parallel_object_only();
308 :
309 51864 : this->add (1., v);
310 51864 : }
311 :
312 :
313 :
314 : template <typename T>
315 374928 : void PetscVector<T>::add (const T a_in, const NumericVector<T> & v_in)
316 : {
317 124982 : parallel_object_only();
318 :
319 374928 : this->_restore_array();
320 :
321 : // VecAXPY doesn't support &x==&y
322 374928 : if (this == &v_in)
323 : {
324 12 : this->scale(a_in+1);
325 12 : 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 374916 : v->_restore_array();
333 :
334 124978 : libmesh_assert_equal_to (this->size(), v->size());
335 :
336 374916 : 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 8352 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
343 :
344 374916 : this->_is_closed = true;
345 : }
346 :
347 :
348 :
349 : template <typename T>
350 1385493 : void PetscVector<T>::insert (const T * v,
351 : const std::vector<numeric_index_type> & dof_indices)
352 : {
353 1385493 : if (dof_indices.empty())
354 0 : return;
355 :
356 1385493 : this->_restore_array();
357 :
358 923670 : PetscInt * idx_values = numeric_petsc_cast(dof_indices.data());
359 1385493 : std::scoped_lock lock(this->_numeric_vector_mutex);
360 1847328 : LibmeshPetscCall(VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
361 : idx_values, pPS(v), INSERT_VALUES));
362 :
363 1385493 : this->_is_closed = false;
364 : }
365 :
366 :
367 :
368 : template <typename T>
369 54228 : void PetscVector<T>::scale (const T factor_in)
370 : {
371 18076 : parallel_object_only();
372 :
373 54228 : this->_restore_array();
374 :
375 18076 : PetscScalar factor = PS(factor_in);
376 :
377 54228 : 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 54228 : }
385 :
386 : template <typename T>
387 12 : 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 12 : LibmeshPetscCall(VecPointwiseMult(_vec, _vec, v_vec->_vec));
395 :
396 12 : return *this;
397 : }
398 :
399 : template <typename T>
400 7258 : 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 7258 : LibmeshPetscCall(VecPointwiseDivide(_vec, _vec, v_vec->_vec));
408 :
409 7258 : 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 888348 : T PetscVector<T>::dot (const NumericVector<T> & v_in) const
430 : {
431 296116 : parallel_object_only();
432 :
433 888348 : this->_restore_array();
434 :
435 : // Error flag
436 :
437 : // Return value
438 888348 : 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 888348 : LibmeshPetscCall(VecDot(this->_vec, v->_vec, &value));
445 :
446 888348 : 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 12 : PetscVector<T>::operator = (const T s_in)
474 : {
475 4 : parallel_object_only();
476 :
477 12 : this->_restore_array();
478 4 : libmesh_assert(this->closed());
479 :
480 4 : PetscScalar s = PS(s_in);
481 :
482 12 : if (this->size() != 0)
483 : {
484 12 : 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 12 : return *this;
494 : }
495 :
496 :
497 :
498 : template <typename T>
499 : NumericVector<T> &
500 333832 : 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 333832 : *this = *v;
508 :
509 333832 : return *this;
510 : }
511 :
512 :
513 :
514 : template <typename T>
515 : PetscVector<T> &
516 333838 : 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 333838 : if (this == &v)
527 480 : return *this;
528 :
529 332398 : this->_restore_array();
530 332398 : 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 332254 : 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 144 : LibmeshPetscCall(VecScatterCreateToAll(v._vec, scatter.get(), nullptr));
554 144 : 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 332254 : 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 245046 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
574 :
575 332398 : this->_is_closed = true;
576 :
577 332398 : return *this;
578 : }
579 :
580 :
581 :
582 : template <typename T>
583 : NumericVector<T> &
584 6 : PetscVector<T>::operator = (const std::vector<T> & v)
585 : {
586 2 : parallel_object_only();
587 :
588 6 : 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 6 : if (this->size() == v.size())
595 : {
596 6 : numeric_index_type first = first_local_index();
597 6 : numeric_index_type last = last_local_index();
598 25788 : for (numeric_index_type i=0; i<last-first; i++)
599 34376 : _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 6 : return *this;
617 : }
618 :
619 :
620 :
621 : template <typename T>
622 124774 : 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 124774 : v_local_in = *this;
632 124774 : }
633 :
634 :
635 :
636 : template <typename T>
637 191098 : 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 125092 : }
642 :
643 :
644 :
645 : template <typename T>
646 12 : 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 16 : 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 12 : indices.empty() ? nullptr : numeric_petsc_cast(indices.data());
662 8 : WrappedPetsc<IS> is;
663 12 : 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 12 : 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 12 : VecScatterBeginEnd(this->comm(), scatter, _vec, dest, INSERT_VALUES, SCATTER_FORWARD);
676 :
677 : // Get access to the values stored in dest.
678 : PetscScalar * values;
679 12 : 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 16 : v_local.reserve(indices.size());
685 4 : v_local.clear();
686 16 : v_local.insert(v_local.begin(), values, values+indices.size());
687 :
688 : // We are done using it, so restore the array.
689 12 : LibmeshPetscCall(VecRestoreArray (dest, &values));
690 12 : }
691 :
692 :
693 :
694 : template <typename T>
695 1812 : 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 1812 : 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 1812 : const numeric_index_type my_size = this->size();
707 1812 : 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 2484 : if (this->n_processors() == 1)
714 0 : return;
715 :
716 :
717 : // Build a parallel vector, initialize it with the local
718 : // parts of (*this)
719 3156 : PetscVector<T> parallel_vec(this->comm(), PARALLEL);
720 :
721 1812 : 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 2484 : 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 2928 : 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 1812 : LibmeshPetscCall(VecScatterCreate(_vec, is,
737 : parallel_vec._vec, is,
738 : scatter.get()));
739 :
740 : // Perform the scatter
741 1812 : VecScatterBeginEnd(this->comm(), scatter, _vec, parallel_vec._vec, INSERT_VALUES, SCATTER_FORWARD);
742 : }
743 :
744 : // localize like normal
745 1812 : parallel_vec.close();
746 1812 : parallel_vec.localize (*this, send_list);
747 1812 : this->close();
748 468 : }
749 :
750 :
751 :
752 : template <typename T>
753 1830 : void PetscVector<T>::localize (std::vector<T> & v_local) const
754 : {
755 610 : parallel_object_only();
756 :
757 1830 : this->_restore_array();
758 :
759 : // This function must be run on all processors at once
760 610 : parallel_object_only();
761 :
762 1830 : const PetscInt n = this->size();
763 1830 : const PetscInt nl = this->local_size();
764 : PetscScalar * values;
765 :
766 610 : v_local.clear();
767 1830 : v_local.resize(n, 0.);
768 :
769 1830 : LibmeshPetscCall(VecGetArray (_vec, &values));
770 :
771 1830 : numeric_index_type ioff = first_local_index();
772 :
773 3594273 : for (PetscInt i=0; i<nl; i++)
774 4789924 : v_local[i+ioff] = static_cast<T>(values[i]);
775 :
776 1830 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
777 :
778 1830 : this->comm().sum(v_local);
779 1830 : }
780 :
781 :
782 :
783 : // Full specialization for Real datatypes
784 : #ifdef LIBMESH_USE_REAL_NUMBERS
785 :
786 : template <>
787 7948 : 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 7948 : this->_restore_array();
794 :
795 7948 : const PetscInt n = size();
796 : PetscScalar * values;
797 :
798 : // only one processor
799 10642 : if (n_processors() == 1)
800 : {
801 0 : v_local.resize(n);
802 :
803 0 : LibmeshPetscCall(VecGetArray (_vec, &values));
804 :
805 0 : for (PetscInt i=0; i<n; i++)
806 0 : v_local[i] = static_cast<Real>(values[i]);
807 :
808 0 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
809 : }
810 :
811 : // otherwise multiple processors
812 : #ifdef LIBMESH_HAVE_MPI
813 : else
814 : {
815 7948 : if (pid == 0) // optimized version for localizing to 0
816 : {
817 5388 : WrappedPetsc<Vec> vout;
818 5388 : WrappedPetsc<VecScatter> ctx;
819 :
820 7948 : LibmeshPetscCall(VecScatterCreateToZero(_vec, ctx.get(), vout.get()));
821 :
822 7948 : VecScatterBeginEnd(this->comm(), ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
823 :
824 10642 : if (processor_id() == 0)
825 : {
826 3974 : v_local.resize(n);
827 :
828 3974 : LibmeshPetscCall(VecGetArray (vout, &values));
829 :
830 8235040 : for (PetscInt i=0; i<n; i++)
831 10922062 : v_local[i] = static_cast<Real>(values[i]);
832 :
833 3974 : 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 7948 : }
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 6 : void PetscVector<T>::pointwise_mult (const NumericVector<T> & vec1,
956 : const NumericVector<T> & vec2)
957 : {
958 2 : parallel_object_only();
959 :
960 6 : 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 6 : 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 6 : this->_is_closed = true;
976 6 : }
977 :
978 : template <typename T>
979 6 : void PetscVector<T>::pointwise_divide (const NumericVector<T> & vec1,
980 : const NumericVector<T> & vec2)
981 : {
982 2 : parallel_object_only();
983 :
984 6 : 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 6 : 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 6 : this->_is_closed = true;
1000 6 : }
1001 :
1002 : template <typename T>
1003 12 : void PetscVector<T>::print_matlab (const std::string & name) const
1004 : {
1005 4 : parallel_object_only();
1006 :
1007 12 : 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 12 : if (name != "")
1013 : {
1014 8 : WrappedPetsc<PetscViewer> petsc_viewer;
1015 :
1016 12 : 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 12 : LibmeshPetscCall(PetscViewerPushFormat (petsc_viewer,
1025 : PETSC_VIEWER_ASCII_MATLAB));
1026 : #endif
1027 :
1028 12 : 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 12 : }
1046 :
1047 :
1048 :
1049 :
1050 :
1051 : template <typename T>
1052 738 : 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 738 : 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 738 : 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 738 : petsc_subvector->_restore_array();
1140 : }
1141 :
1142 1230 : std::vector<PetscInt> idx(rows.size());
1143 738 : if (supplying_global_rows)
1144 4 : std::iota (idx.begin(), idx.end(), 0);
1145 : else
1146 : {
1147 : PetscInt start;
1148 726 : LibmeshPetscCall(VecGetOwnershipRange(petsc_subvector->_vec, &start, nullptr));
1149 726 : std::iota (idx.begin(), idx.end(), start);
1150 : }
1151 :
1152 : // Construct index sets
1153 492 : WrappedPetsc<IS> parent_is;
1154 984 : 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 984 : 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 738 : LibmeshPetscCall(VecScatterCreate(this->_vec,
1170 : parent_is,
1171 : petsc_subvector->_vec,
1172 : subvector_is,
1173 : scatter.get()));
1174 :
1175 : // Actually perform the scatter
1176 738 : VecScatterBeginEnd(this->comm(), scatter, this->_vec, petsc_subvector->_vec, INSERT_VALUES, SCATTER_FORWARD);
1177 :
1178 738 : petsc_subvector->_is_closed = true;
1179 738 : }
1180 :
1181 :
1182 :
1183 : template <typename T>
1184 1889369388 : void PetscVector<T>::_get_array(bool read_only) const
1185 : {
1186 589769070 : libmesh_assert (this->initialized());
1187 :
1188 589769070 : 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 1889369388 : if (initially_array_is_present && read_only && !_values_read_only)
1193 6765258 : _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 1889369388 : if (initially_array_is_present && !read_only && _values_read_only)
1200 : {
1201 18 : _restore_array();
1202 6 : initially_array_is_present = false;
1203 : }
1204 :
1205 1889369376 : if (!initially_array_is_present)
1206 : {
1207 453847 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1208 338271 : if (!_array_is_present.load(std::memory_order_relaxed))
1209 : {
1210 115099 : if (!this->is_effectively_ghosted())
1211 : {
1212 4339 : if (read_only)
1213 : {
1214 4289 : LibmeshPetscCall(VecGetArrayRead(_vec, &_read_only_values));
1215 4289 : _values_read_only = true;
1216 : }
1217 : else
1218 : {
1219 50 : LibmeshPetscCall(VecGetArray(_vec, &_values));
1220 50 : _values_read_only = false;
1221 : }
1222 4339 : _local_size = this->local_size();
1223 : }
1224 : else
1225 : {
1226 332688 : LibmeshPetscCall(VecGhostGetLocalForm (_vec,&_local_form));
1227 :
1228 332688 : if (read_only)
1229 : {
1230 332688 : LibmeshPetscCall(VecGetArrayRead(_local_form, &_read_only_values));
1231 332688 : _values_read_only = true;
1232 : }
1233 : else
1234 : {
1235 0 : LibmeshPetscCall(VecGetArray(_local_form, &_values));
1236 0 : _values_read_only = false;
1237 : }
1238 :
1239 332688 : PetscInt my_local_size = 0;
1240 332688 : LibmeshPetscCall(VecGetLocalSize(_local_form, &my_local_size));
1241 332688 : _local_size = static_cast<numeric_index_type>(my_local_size);
1242 : }
1243 :
1244 : { // cache ownership range
1245 337027 : PetscInt petsc_first=0, petsc_last=0;
1246 337027 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1247 337027 : _first = static_cast<numeric_index_type>(petsc_first);
1248 337027 : _last = static_cast<numeric_index_type>(petsc_last);
1249 : }
1250 115099 : _array_is_present.store(true, std::memory_order_release);
1251 : }
1252 : }
1253 1889369388 : }
1254 :
1255 :
1256 :
1257 : template <typename T>
1258 86991590 : void PetscVector<T>::_restore_array() const
1259 : {
1260 86991590 : 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 86991590 : if (_array_is_present.load(std::memory_order_acquire))
1265 : {
1266 452126 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1267 337027 : if (_array_is_present.load(std::memory_order_relaxed))
1268 : {
1269 115099 : if (!this->is_effectively_ghosted())
1270 : {
1271 4339 : if (_values_read_only)
1272 4289 : LibmeshPetscCall(VecRestoreArrayRead (_vec, &_read_only_values));
1273 : else
1274 50 : LibmeshPetscCall(VecRestoreArray (_vec, &_values));
1275 :
1276 4339 : _values = nullptr;
1277 : }
1278 : else
1279 : {
1280 332688 : if (_values_read_only)
1281 332688 : LibmeshPetscCall(VecRestoreArrayRead (_local_form, &_read_only_values));
1282 : else
1283 0 : LibmeshPetscCall(VecRestoreArray (_local_form, &_values));
1284 :
1285 332688 : _values = nullptr;
1286 332688 : LibmeshPetscCall(VecGhostRestoreLocalForm (_vec,&_local_form));
1287 332688 : _local_form = nullptr;
1288 332688 : _local_size = 0;
1289 : }
1290 115099 : _array_is_present.store(false, std::memory_order_release);
1291 : }
1292 : }
1293 86991590 : }
1294 :
1295 : template <typename T>
1296 : std::unique_ptr<NumericVector<T>>
1297 946 : PetscVector<T>::get_subvector(const std::vector<numeric_index_type> & rows)
1298 : {
1299 : // Construct index set
1300 242 : WrappedPetsc<IS> parent_is;
1301 1188 : 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 946 : LibmeshPetscCall(VecGetSubVector(_vec, parent_is, &subvec));
1309 :
1310 946 : this->_is_closed = false;
1311 :
1312 1188 : return std::make_unique<PetscVector<T>>(subvec, this->comm());
1313 : }
1314 :
1315 : template <typename T>
1316 : void
1317 946 : 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 1188 : 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 946 : Vec subvec = petsc_subvector->vec();
1331 946 : 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 946 : this->_is_closed = true;
1337 946 : }
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
|
