Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 276316 : Real PetscVector<T>::norm () const
61 : {
62 8070 : parallel_object_only();
63 :
64 276316 : this->_restore_array();
65 8070 : libmesh_assert(this->closed());
66 :
67 276316 : PetscReal value=0.;
68 :
69 276316 : LibmeshPetscCall(VecNorm (_vec, N, &value));
70 :
71 276316 : 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 274566 : Real PetscVector<T>::l2_norm () const
80 : {
81 274566 : 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 38902412 : void PetscVector<T>::set (const numeric_index_type i, const T value)
125 : {
126 38902412 : this->_restore_array();
127 3655214 : libmesh_assert_less (i, size());
128 :
129 38902412 : PetscInt i_val = static_cast<PetscInt>(i);
130 38902412 : PetscScalar petsc_value = PS(value);
131 :
132 38902412 : std::scoped_lock lock(this->_numeric_vector_mutex);
133 38902412 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES));
134 :
135 38902412 : this->_is_closed = false;
136 38902412 : }
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 190919550 : void PetscVector<T>::add (const numeric_index_type i, const T value)
166 : {
167 190919550 : this->_restore_array();
168 17728615 : libmesh_assert_less (i, size());
169 :
170 190919550 : PetscInt i_val = static_cast<PetscInt>(i);
171 190919550 : PetscScalar petsc_value = PS(value);
172 :
173 190919550 : std::scoped_lock lock(this->_numeric_vector_mutex);
174 190919550 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES));
175 :
176 190919550 : this->_is_closed = false;
177 190919550 : }
178 :
179 :
180 :
181 : template <typename T>
182 74028276 : 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 74028276 : if (dof_indices.empty())
187 0 : return;
188 :
189 74028276 : this->_restore_array();
190 :
191 11253401 : const PetscInt * i_val = reinterpret_cast<const PetscInt *>(dof_indices.data());
192 5731173 : const PetscScalar * petsc_value = pPS(v);
193 :
194 74028276 : std::scoped_lock lock(this->_numeric_vector_mutex);
195 79550504 : LibmeshPetscCall(VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
196 : i_val, petsc_value, ADD_VALUES));
197 :
198 74028276 : 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 4370942 : void PetscVector<T>::add (const T a_in, const NumericVector<T> & v_in)
316 : {
317 125044 : parallel_object_only();
318 :
319 4370942 : this->_restore_array();
320 :
321 : // VecAXPY doesn't support &x==&y
322 4370942 : if (this == &v_in)
323 : {
324 140 : this->scale(a_in+1);
325 140 : return;
326 : }
327 :
328 125040 : PetscScalar a = PS(a_in);
329 :
330 : // Make sure the NumericVector passed in is really a PetscVector
331 125040 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
332 4370802 : v->_restore_array();
333 :
334 125040 : libmesh_assert_equal_to (this->size(), v->size());
335 :
336 4370802 : LibmeshPetscCall(VecAXPY(_vec, a, v->vec()));
337 :
338 125040 : libmesh_assert(this->comm().verify(int(this->type())));
339 :
340 4370802 : if (this->type() == GHOSTED)
341 97440 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
342 :
343 4370802 : this->_is_closed = true;
344 : }
345 :
346 :
347 :
348 : template <typename T>
349 4764364 : void PetscVector<T>::insert (const T * v,
350 : const std::vector<numeric_index_type> & dof_indices)
351 : {
352 4764364 : if (dof_indices.empty())
353 0 : return;
354 :
355 4764364 : this->_restore_array();
356 :
357 865384 : PetscInt * idx_values = numeric_petsc_cast(dof_indices.data());
358 4764364 : std::scoped_lock lock(this->_numeric_vector_mutex);
359 5197056 : LibmeshPetscCall(VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
360 : idx_values, pPS(v), INSERT_VALUES));
361 :
362 4764364 : this->_is_closed = false;
363 : }
364 :
365 :
366 :
367 : template <typename T>
368 632896 : void PetscVector<T>::scale (const T factor_in)
369 : {
370 18076 : parallel_object_only();
371 :
372 632896 : this->_restore_array();
373 :
374 18076 : PetscScalar factor = PS(factor_in);
375 :
376 632896 : LibmeshPetscCall(VecScale(_vec, factor));
377 :
378 18076 : libmesh_assert(this->comm().verify(int(this->type())));
379 :
380 632896 : if (this->type() == GHOSTED)
381 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
382 632896 : }
383 :
384 : template <typename T>
385 140 : NumericVector<T> & PetscVector<T>::operator *= (const NumericVector<T> & v)
386 : {
387 4 : parallel_object_only();
388 :
389 :
390 4 : const PetscVector<T> * v_vec = cast_ptr<const PetscVector<T> *>(&v);
391 :
392 140 : LibmeshPetscCall(VecPointwiseMult(_vec, _vec, v_vec->_vec));
393 :
394 140 : return *this;
395 : }
396 :
397 : template <typename T>
398 78870 : NumericVector<T> & PetscVector<T>::operator /= (const NumericVector<T> & v)
399 : {
400 2454 : parallel_object_only();
401 :
402 :
403 2454 : const PetscVector<T> * v_vec = cast_ptr<const PetscVector<T> *>(&v);
404 :
405 78870 : LibmeshPetscCall(VecPointwiseDivide(_vec, _vec, v_vec->_vec));
406 :
407 78870 : return *this;
408 : }
409 :
410 : template <typename T>
411 0 : void PetscVector<T>::abs()
412 : {
413 0 : parallel_object_only();
414 :
415 0 : this->_restore_array();
416 :
417 0 : LibmeshPetscCall(VecAbs(_vec));
418 :
419 0 : libmesh_assert(this->comm().verify(int(this->type())));
420 :
421 0 : if (this->type() == GHOSTED)
422 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
423 0 : }
424 :
425 : template <typename T>
426 10363074 : T PetscVector<T>::dot (const NumericVector<T> & v_in) const
427 : {
428 296116 : parallel_object_only();
429 :
430 10363074 : this->_restore_array();
431 :
432 : // Error flag
433 :
434 : // Return value
435 10363074 : PetscScalar value=0.;
436 :
437 : // Make sure the NumericVector passed in is really a PetscVector
438 296116 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
439 :
440 : // 2.3.x (at least) style. Untested for previous versions.
441 10363074 : LibmeshPetscCall(VecDot(this->_vec, v->_vec, &value));
442 :
443 10363074 : return static_cast<T>(value);
444 : }
445 :
446 : template <typename T>
447 0 : T PetscVector<T>::indefinite_dot (const NumericVector<T> & v_in) const
448 : {
449 0 : parallel_object_only();
450 :
451 0 : this->_restore_array();
452 :
453 : // Error flag
454 :
455 : // Return value
456 0 : PetscScalar value=0.;
457 :
458 : // Make sure the NumericVector passed in is really a PetscVector
459 0 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
460 :
461 : // 2.3.x (at least) style. Untested for previous versions.
462 0 : LibmeshPetscCall(VecTDot(this->_vec, v->_vec, &value));
463 :
464 0 : return static_cast<T>(value);
465 : }
466 :
467 :
468 : template <typename T>
469 : NumericVector<T> &
470 140 : PetscVector<T>::operator = (const T s_in)
471 : {
472 4 : parallel_object_only();
473 :
474 140 : this->_restore_array();
475 4 : libmesh_assert(this->closed());
476 :
477 4 : PetscScalar s = PS(s_in);
478 :
479 140 : if (this->size() != 0)
480 : {
481 140 : LibmeshPetscCall(VecSet(_vec, s));
482 :
483 4 : libmesh_assert(this->comm().verify(int(this->type())));
484 :
485 140 : if (this->type() == GHOSTED)
486 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
487 : }
488 :
489 140 : return *this;
490 : }
491 :
492 :
493 :
494 : template <typename T>
495 : NumericVector<T> &
496 3726025 : PetscVector<T>::operator = (const NumericVector<T> & v_in)
497 : {
498 104240 : parallel_object_only();
499 :
500 : // Make sure the NumericVector passed in is really a PetscVector
501 104240 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
502 :
503 3726025 : *this = *v;
504 :
505 3726025 : return *this;
506 : }
507 :
508 :
509 :
510 : template <typename T>
511 : PetscVector<T> &
512 3726095 : PetscVector<T>::operator = (const PetscVector<T> & v)
513 : {
514 104242 : parallel_object_only();
515 :
516 3726095 : this->_restore_array();
517 3726095 : v._restore_array();
518 :
519 104242 : libmesh_assert_equal_to (this->size(), v.size());
520 104242 : libmesh_assert_equal_to (this->local_size(), v.local_size());
521 104242 : libmesh_assert (v.closed());
522 :
523 :
524 3726095 : LibmeshPetscCall(VecCopy (v._vec, this->_vec));
525 :
526 104242 : libmesh_assert(this->comm().verify(int(this->type())));
527 :
528 3726095 : if (this->type() == GHOSTED)
529 2721427 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
530 :
531 3726095 : this->_is_closed = true;
532 :
533 3726095 : return *this;
534 : }
535 :
536 :
537 :
538 : template <typename T>
539 : NumericVector<T> &
540 70 : PetscVector<T>::operator = (const std::vector<T> & v)
541 : {
542 2 : parallel_object_only();
543 :
544 70 : this->_get_array(false);
545 :
546 : /**
547 : * Case 1: The vector is the same size of
548 : * The global vector. Only add the local components.
549 : */
550 70 : if (this->size() == v.size())
551 : {
552 70 : numeric_index_type first = first_local_index();
553 70 : numeric_index_type last = last_local_index();
554 300860 : for (numeric_index_type i=0; i<last-first; i++)
555 309384 : _values[i] = PS(v[first + i]);
556 : }
557 :
558 : /**
559 : * Case 2: The vector is the same size as our local
560 : * piece. Insert directly to the local piece.
561 : */
562 : else
563 : {
564 0 : for (numeric_index_type i=0; i<_local_size; i++)
565 0 : _values[i] = PS(v[i]);
566 : }
567 :
568 : // Make sure ghost dofs are up to date
569 70 : if (this->type() == GHOSTED)
570 0 : this->close();
571 :
572 70 : return *this;
573 : }
574 :
575 :
576 :
577 : template <typename T>
578 2150 : void PetscVector<T>::localize (NumericVector<T> & v_local_in) const
579 : {
580 50 : parallel_object_only();
581 :
582 2150 : this->_restore_array();
583 :
584 : // Make sure the NumericVector passed in is really a PetscVector
585 50 : PetscVector<T> * v_local = cast_ptr<PetscVector<T> *>(&v_local_in);
586 :
587 50 : libmesh_assert(v_local);
588 : // v_local_in should be closed
589 50 : libmesh_assert(v_local->closed());
590 50 : libmesh_assert_equal_to (v_local->size(), this->size());
591 : // 1) v_local_in is a large vector to hold the whole world
592 : // 2) v_local_in is a ghosted vector
593 : // 3) v_local_in is a parallel vector
594 : // Cases 2) and 3) should be scalable
595 50 : libmesh_assert(this->size()==v_local->local_size() || this->local_size()==v_local->local_size());
596 :
597 :
598 2150 : if (v_local->type() == SERIAL && this->size() == v_local->local_size())
599 : {
600 48 : WrappedPetsc<VecScatter> scatter;
601 2080 : LibmeshPetscCall(VecScatterCreateToAll(_vec, scatter.get(), nullptr));
602 2080 : VecScatterBeginEnd(this->comm(), scatter, _vec, v_local->_vec, INSERT_VALUES, SCATTER_FORWARD);
603 : }
604 : // Two vectors have the same size, and we should just do a simple copy.
605 : // v_local could be either a parallel or ghosted vector
606 70 : else if (this->local_size() == v_local->local_size())
607 70 : LibmeshPetscCall(VecCopy(_vec,v_local->_vec));
608 : else
609 0 : libmesh_error_msg("Vectors are inconsistent");
610 :
611 : // Make sure ghost dofs are up to date
612 : // We do not call "close" here to save a global reduction
613 2150 : if (v_local->type() == GHOSTED)
614 70 : VecGhostUpdateBeginEnd(this->comm(), v_local->_vec, INSERT_VALUES, SCATTER_FORWARD);
615 2150 : }
616 :
617 :
618 :
619 : template <typename T>
620 2244046 : void PetscVector<T>::localize (NumericVector<T> & v_local_in,
621 : const std::vector<numeric_index_type> & send_list) const
622 : {
623 61584 : parallel_object_only();
624 :
625 61584 : libmesh_assert(this->comm().verify(int(this->type())));
626 61584 : libmesh_assert(this->comm().verify(int(v_local_in.type())));
627 :
628 : // FIXME: Workaround for a strange bug at large-scale.
629 : // If we have ghosting, PETSc lets us just copy the solution, and
630 : // doing so avoids a segfault?
631 2305630 : if (v_local_in.type() == GHOSTED &&
632 122164 : this->type() == PARALLEL)
633 : {
634 2112710 : v_local_in = *this;
635 2112710 : return;
636 : }
637 :
638 : // Normal code path begins here
639 :
640 131336 : this->_restore_array();
641 :
642 : // Make sure the NumericVector passed in is really a PetscVector
643 3860 : PetscVector<T> * v_local = cast_ptr<PetscVector<T> *>(&v_local_in);
644 :
645 3860 : libmesh_assert(v_local);
646 3860 : libmesh_assert_equal_to (v_local->size(), this->size());
647 3860 : libmesh_assert_less_equal (send_list.size(), v_local->size());
648 :
649 : const numeric_index_type n_sl =
650 7532 : cast_int<numeric_index_type>(send_list.size());
651 :
652 135196 : std::vector<PetscInt> idx(n_sl + this->local_size());
653 2982828 : for (numeric_index_type i=0; i<n_sl; i++)
654 3176636 : idx[i] = static_cast<PetscInt>(send_list[i]);
655 26086052 : for (auto i : make_range(this->local_size()))
656 25954716 : idx[n_sl+i] = i + this->first_local_index();
657 :
658 : // Create the index set & scatter objects
659 7720 : WrappedPetsc<IS> is;
660 131336 : PetscInt * idxptr = idx.empty() ? nullptr : idx.data();
661 131336 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(), n_sl+this->local_size(),
662 : idxptr, PETSC_USE_POINTER, is.get()));
663 :
664 7720 : WrappedPetsc<VecScatter> scatter;
665 131336 : LibmeshPetscCall(VecScatterCreate(_vec, is,
666 : v_local->_vec, is,
667 : scatter.get()));
668 :
669 : // Perform the scatter
670 131336 : VecScatterBeginEnd(this->comm(), scatter, _vec, v_local->_vec, INSERT_VALUES, SCATTER_FORWARD);
671 :
672 : // Make sure ghost dofs are up to date
673 131336 : if (v_local->type() == GHOSTED)
674 131336 : v_local->close();
675 : }
676 :
677 :
678 :
679 : template <typename T>
680 140 : void PetscVector<T>::localize (std::vector<T> & v_local,
681 : const std::vector<numeric_index_type> & indices) const
682 : {
683 4 : parallel_object_only();
684 :
685 : // Error code used to check the status of all PETSc function calls.
686 :
687 : // Create a sequential destination Vec with the right number of entries on each proc.
688 8 : WrappedPetsc<Vec> dest;
689 144 : LibmeshPetscCall(VecCreateSeq(PETSC_COMM_SELF, cast_int<PetscInt>(indices.size()), dest.get()));
690 :
691 : // Create an IS using the libmesh routine. PETSc does not own the
692 : // IS memory in this case, it is automatically cleaned up by the
693 : // std::vector destructor.
694 8 : PetscInt * idxptr =
695 140 : indices.empty() ? nullptr : numeric_petsc_cast(indices.data());
696 8 : WrappedPetsc<IS> is;
697 140 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(), cast_int<PetscInt>(indices.size()), idxptr,
698 : PETSC_USE_POINTER, is.get()));
699 :
700 : // Create the VecScatter object. "PETSC_NULL" means "use the identity IS".
701 8 : WrappedPetsc<VecScatter> scatter;
702 140 : LibmeshPetscCall(VecScatterCreate(_vec,
703 : /*src is=*/is,
704 : /*dest vec=*/dest,
705 : /*dest is=*/LIBMESH_PETSC_NULLPTR,
706 : scatter.get()));
707 :
708 : // Do the scatter
709 140 : VecScatterBeginEnd(this->comm(), scatter, _vec, dest, INSERT_VALUES, SCATTER_FORWARD);
710 :
711 : // Get access to the values stored in dest.
712 : PetscScalar * values;
713 140 : LibmeshPetscCall(VecGetArray (dest, &values));
714 :
715 : // Store values into the provided v_local. Make sure there is enough
716 : // space reserved and then clear out any existing entries before
717 : // inserting.
718 144 : v_local.reserve(indices.size());
719 4 : v_local.clear();
720 144 : v_local.insert(v_local.begin(), values, values+indices.size());
721 :
722 : // We are done using it, so restore the array.
723 140 : LibmeshPetscCall(VecRestoreArray (dest, &values));
724 140 : }
725 :
726 :
727 :
728 : template <typename T>
729 10428 : void PetscVector<T>::localize (const numeric_index_type first_local_idx,
730 : const numeric_index_type last_local_idx,
731 : const std::vector<numeric_index_type> & send_list)
732 : {
733 472 : parallel_object_only();
734 :
735 10428 : this->_restore_array();
736 :
737 472 : libmesh_assert_less_equal (send_list.size(), this->size());
738 472 : libmesh_assert_less_equal (last_local_idx+1, this->size());
739 :
740 10428 : const numeric_index_type my_size = this->size();
741 10428 : const numeric_index_type my_local_size = (last_local_idx + 1 - first_local_idx);
742 :
743 : // Don't bother for serial cases
744 : // if ((first_local_idx == 0) &&
745 : // (my_local_size == my_size))
746 : // But we do need to stay in sync for degenerate cases
747 10900 : if (this->n_processors() == 1)
748 144 : return;
749 :
750 :
751 : // Build a parallel vector, initialize it with the local
752 : // parts of (*this)
753 11228 : PetscVector<T> parallel_vec(this->comm(), PARALLEL);
754 :
755 10284 : parallel_vec.init (my_size, my_local_size, true, PARALLEL);
756 :
757 :
758 : // Copy part of *this into the parallel_vec
759 : {
760 : // Create idx, idx[i] = i+first_local_idx;
761 10756 : std::vector<PetscInt> idx(my_local_size);
762 472 : std::iota (idx.begin(), idx.end(), first_local_idx);
763 :
764 : // Create the index set & scatter objects
765 944 : WrappedPetsc<IS> is;
766 18870 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(), my_local_size,
767 : my_local_size ? idx.data() : nullptr, PETSC_USE_POINTER, is.get()));
768 :
769 472 : WrappedPetsc<VecScatter> scatter;
770 10284 : LibmeshPetscCall(VecScatterCreate(_vec, is,
771 : parallel_vec._vec, is,
772 : scatter.get()));
773 :
774 : // Perform the scatter
775 10284 : VecScatterBeginEnd(this->comm(), scatter, _vec, parallel_vec._vec, INSERT_VALUES, SCATTER_FORWARD);
776 : }
777 :
778 : // localize like normal
779 10284 : parallel_vec.close();
780 10284 : parallel_vec.localize (*this, send_list);
781 10284 : this->close();
782 9340 : }
783 :
784 :
785 :
786 : template <typename T>
787 21350 : void PetscVector<T>::localize (std::vector<T> & v_local) const
788 : {
789 610 : parallel_object_only();
790 :
791 21350 : this->_restore_array();
792 :
793 : // This function must be run on all processors at once
794 610 : parallel_object_only();
795 :
796 21350 : const PetscInt n = this->size();
797 21350 : const PetscInt nl = this->local_size();
798 : PetscScalar * values;
799 :
800 610 : v_local.clear();
801 21350 : v_local.resize(n, 0.);
802 :
803 21350 : LibmeshPetscCall(VecGetArray (_vec, &values));
804 :
805 21350 : numeric_index_type ioff = first_local_index();
806 :
807 13196761 : for (PetscInt i=0; i<nl; i++)
808 14372892 : v_local[i+ioff] = static_cast<T>(values[i]);
809 :
810 21350 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
811 :
812 21350 : this->comm().sum(v_local);
813 21350 : }
814 :
815 :
816 :
817 : // Full specialization for Real datatypes
818 : #ifdef LIBMESH_USE_REAL_NUMBERS
819 :
820 : template <>
821 84872 : void PetscVector<Real>::localize_to_one (std::vector<Real> & v_local,
822 : const processor_id_type
823 : timpi_mpi_var(pid)) const
824 : {
825 2700 : parallel_object_only();
826 :
827 84872 : this->_restore_array();
828 :
829 84872 : const PetscInt n = size();
830 : PetscScalar * values;
831 :
832 : // only one processor
833 87572 : if (n_processors() == 1)
834 : {
835 1204 : v_local.resize(n);
836 :
837 1204 : LibmeshPetscCall(VecGetArray (_vec, &values));
838 :
839 2689350 : for (PetscInt i=0; i<n; i++)
840 2688146 : v_local[i] = static_cast<Real>(values[i]);
841 :
842 1204 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
843 : }
844 :
845 : // otherwise multiple processors
846 : #ifdef LIBMESH_HAVE_MPI
847 : else
848 : {
849 83668 : if (pid == 0) // optimized version for localizing to 0
850 : {
851 5400 : WrappedPetsc<Vec> vout;
852 5400 : WrappedPetsc<VecScatter> ctx;
853 :
854 83668 : LibmeshPetscCall(VecScatterCreateToZero(_vec, ctx.get(), vout.get()));
855 :
856 83668 : VecScatterBeginEnd(this->comm(), ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
857 :
858 86368 : if (processor_id() == 0)
859 : {
860 12385 : v_local.resize(n);
861 :
862 12385 : LibmeshPetscCall(VecGetArray (vout, &values));
863 :
864 27338944 : for (PetscInt i=0; i<n; i++)
865 30179034 : v_local[i] = static_cast<Real>(values[i]);
866 :
867 12385 : LibmeshPetscCall(VecRestoreArray (vout, &values));
868 : }
869 : }
870 : else
871 : {
872 0 : v_local.resize(n);
873 :
874 0 : numeric_index_type ioff = this->first_local_index();
875 0 : std::vector<Real> local_values (n, 0.);
876 :
877 : {
878 0 : LibmeshPetscCall(VecGetArray (_vec, &values));
879 :
880 0 : const PetscInt nl = local_size();
881 0 : for (PetscInt i=0; i<nl; i++)
882 0 : local_values[i+ioff] = static_cast<Real>(values[i]);
883 :
884 0 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
885 : }
886 :
887 :
888 0 : MPI_Reduce (local_values.data(), v_local.data(), n,
889 0 : Parallel::StandardType<Real>(),
890 : Parallel::OpFunction<Real>::sum(), pid,
891 0 : this->comm().get());
892 : }
893 : }
894 : #endif // LIBMESH_HAVE_MPI
895 84872 : }
896 :
897 : #endif
898 :
899 :
900 : // Full specialization for Complex datatypes
901 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
902 :
903 : template <>
904 : void PetscVector<Complex>::localize_to_one (std::vector<Complex> & v_local,
905 : const processor_id_type pid) const
906 : {
907 : parallel_object_only();
908 :
909 : this->_restore_array();
910 :
911 : const PetscInt n = size();
912 : const PetscInt nl = local_size();
913 : PetscScalar * values;
914 :
915 :
916 : v_local.resize(n);
917 :
918 :
919 : for (PetscInt i=0; i<n; i++)
920 : v_local[i] = 0.;
921 :
922 : // only one processor
923 : if (n_processors() == 1)
924 : {
925 : LibmeshPetscCall(VecGetArray (_vec, &values));
926 :
927 : for (PetscInt i=0; i<n; i++)
928 : v_local[i] = static_cast<Complex>(values[i]);
929 :
930 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
931 : }
932 :
933 : // otherwise multiple processors
934 : else
935 : {
936 : numeric_index_type ioff = this->first_local_index();
937 :
938 : /* in here the local values are stored, acting as send buffer for MPI
939 : * initialize to zero, since we collect using MPI_SUM
940 : */
941 : std::vector<Real> real_local_values(n, 0.);
942 : std::vector<Real> imag_local_values(n, 0.);
943 :
944 : {
945 : LibmeshPetscCall(VecGetArray (_vec, &values));
946 :
947 : // provide my local share to the real and imag buffers
948 : for (PetscInt i=0; i<nl; i++)
949 : {
950 : real_local_values[i+ioff] = static_cast<Complex>(values[i]).real();
951 : imag_local_values[i+ioff] = static_cast<Complex>(values[i]).imag();
952 : }
953 :
954 : LibmeshPetscCall(VecRestoreArray (_vec, &values));
955 : }
956 :
957 : /* have buffers of the real and imaginary part of v_local.
958 : * Once MPI_Reduce() collected all the real and imaginary
959 : * parts in these std::vector<Real>, the values can be
960 : * copied to v_local
961 : */
962 : std::vector<Real> real_v_local(n);
963 : std::vector<Real> imag_v_local(n);
964 :
965 : // collect entries from other proc's in real_v_local, imag_v_local
966 : MPI_Reduce (real_local_values.data(),
967 : real_v_local.data(), n,
968 : Parallel::StandardType<Real>(),
969 : Parallel::OpFunction<Real>::sum(),
970 : pid, this->comm().get());
971 :
972 : MPI_Reduce (imag_local_values.data(),
973 : imag_v_local.data(), n,
974 : Parallel::StandardType<Real>(),
975 : Parallel::OpFunction<Real>::sum(),
976 : pid, this->comm().get());
977 :
978 : // copy real_v_local and imag_v_local to v_local
979 : for (PetscInt i=0; i<n; i++)
980 : v_local[i] = Complex(real_v_local[i], imag_v_local[i]);
981 : }
982 : }
983 :
984 : #endif
985 :
986 :
987 :
988 : template <typename T>
989 70 : void PetscVector<T>::pointwise_mult (const NumericVector<T> & vec1,
990 : const NumericVector<T> & vec2)
991 : {
992 2 : parallel_object_only();
993 :
994 70 : this->_restore_array();
995 :
996 : // Convert arguments to PetscVector*.
997 2 : const PetscVector<T> * vec1_petsc = cast_ptr<const PetscVector<T> *>(&vec1);
998 2 : const PetscVector<T> * vec2_petsc = cast_ptr<const PetscVector<T> *>(&vec2);
999 :
1000 : // Call PETSc function.
1001 70 : LibmeshPetscCall(VecPointwiseMult(_vec, vec1_petsc->vec(), vec2_petsc->vec()));
1002 :
1003 2 : libmesh_assert(this->comm().verify(int(this->type())));
1004 :
1005 70 : if (this->type() == GHOSTED)
1006 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
1007 :
1008 70 : this->_is_closed = true;
1009 70 : }
1010 :
1011 : template <typename T>
1012 70 : void PetscVector<T>::pointwise_divide (const NumericVector<T> & vec1,
1013 : const NumericVector<T> & vec2)
1014 : {
1015 2 : parallel_object_only();
1016 :
1017 70 : this->_restore_array();
1018 :
1019 : // Convert arguments to PetscVector*.
1020 2 : const PetscVector<T> * const vec1_petsc = cast_ptr<const PetscVector<T> *>(&vec1);
1021 2 : const PetscVector<T> * const vec2_petsc = cast_ptr<const PetscVector<T> *>(&vec2);
1022 :
1023 : // Call PETSc function.
1024 70 : LibmeshPetscCall(VecPointwiseDivide(_vec, vec1_petsc->vec(), vec2_petsc->vec()));
1025 :
1026 2 : libmesh_assert(this->comm().verify(int(this->type())));
1027 :
1028 70 : if (this->type() == GHOSTED)
1029 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
1030 :
1031 70 : this->_is_closed = true;
1032 70 : }
1033 :
1034 : template <typename T>
1035 140 : void PetscVector<T>::print_matlab (const std::string & name) const
1036 : {
1037 4 : parallel_object_only();
1038 :
1039 140 : this->_restore_array();
1040 4 : libmesh_assert (this->closed());
1041 :
1042 : // Create an ASCII file containing the matrix
1043 : // if a filename was provided.
1044 140 : if (name != "")
1045 : {
1046 8 : WrappedPetsc<PetscViewer> petsc_viewer;
1047 :
1048 140 : LibmeshPetscCall(PetscViewerASCIIOpen(this->comm().get(),
1049 : name.c_str(),
1050 : petsc_viewer.get()));
1051 :
1052 : #if PETSC_VERSION_LESS_THAN(3,7,0)
1053 : LibmeshPetscCall(PetscViewerSetFormat (petsc_viewer,
1054 : PETSC_VIEWER_ASCII_MATLAB));
1055 : #else
1056 140 : LibmeshPetscCall(PetscViewerPushFormat (petsc_viewer,
1057 : PETSC_VIEWER_ASCII_MATLAB));
1058 : #endif
1059 :
1060 140 : LibmeshPetscCall(VecView (_vec, petsc_viewer));
1061 : }
1062 :
1063 : // Otherwise the matrix will be dumped to the screen.
1064 : else
1065 : {
1066 :
1067 : #if PETSC_VERSION_LESS_THAN(3,7,0)
1068 : LibmeshPetscCall(PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
1069 : PETSC_VIEWER_ASCII_MATLAB));
1070 : #else
1071 0 : LibmeshPetscCall(PetscViewerPushFormat (PETSC_VIEWER_STDOUT_WORLD,
1072 : PETSC_VIEWER_ASCII_MATLAB));
1073 : #endif
1074 :
1075 0 : LibmeshPetscCall(VecView (_vec, PETSC_VIEWER_STDOUT_WORLD));
1076 : }
1077 140 : }
1078 :
1079 :
1080 :
1081 :
1082 :
1083 : template <typename T>
1084 1330 : void PetscVector<T>::create_subvector(NumericVector<T> & subvector,
1085 : const std::vector<numeric_index_type> & rows,
1086 : const bool supplying_global_rows) const
1087 : {
1088 38 : parallel_object_only();
1089 :
1090 1330 : libmesh_error_msg_if(
1091 : subvector.type() == GHOSTED,
1092 : "We do not support scattering parallel information to ghosts for subvectors");
1093 :
1094 1330 : this->_restore_array();
1095 :
1096 : // PETSc data structures
1097 :
1098 : // Make sure the passed in subvector is really a PetscVector
1099 38 : PetscVector<T> * petsc_subvector = cast_ptr<PetscVector<T> *>(&subvector);
1100 :
1101 : // If the petsc_subvector is already initialized, we assume that the
1102 : // user has already allocated the *correct* amount of space for it.
1103 : // If not, we use the appropriate PETSc routines to initialize it.
1104 1330 : if (!petsc_subvector->initialized())
1105 : {
1106 0 : libmesh_assert(petsc_subvector->_type == AUTOMATIC || petsc_subvector->_type == PARALLEL);
1107 :
1108 0 : if (supplying_global_rows)
1109 : // Initialize the petsc_subvector to have enough space to hold
1110 : // the entries which will be scattered into it. Note: such an
1111 : // init() function (where we let PETSc decide the number of local
1112 : // entries) is not currently offered by the PetscVector
1113 : // class. Should we differentiate here between sequential and
1114 : // parallel vector creation based on this->n_processors() ?
1115 0 : LibmeshPetscCall(VecCreateMPI(this->comm().get(),
1116 : PETSC_DECIDE, // n_local
1117 : cast_int<PetscInt>(rows.size()), // n_global
1118 : &(petsc_subvector->_vec)));
1119 : else
1120 0 : LibmeshPetscCall(VecCreateMPI(this->comm().get(),
1121 : cast_int<PetscInt>(rows.size()),
1122 : PETSC_DETERMINE,
1123 : &(petsc_subvector->_vec)));
1124 :
1125 0 : LibmeshPetscCall(VecSetFromOptions (petsc_subvector->_vec));
1126 :
1127 : // We created a parallel vector
1128 0 : petsc_subvector->_type = PARALLEL;
1129 :
1130 : // Mark the subvector as initialized
1131 0 : petsc_subvector->_is_initialized = true;
1132 : }
1133 : else
1134 : {
1135 1330 : petsc_subvector->_restore_array();
1136 : }
1137 :
1138 1406 : std::vector<PetscInt> idx(rows.size());
1139 1330 : if (supplying_global_rows)
1140 0 : std::iota (idx.begin(), idx.end(), 0);
1141 : else
1142 : {
1143 : PetscInt start;
1144 1330 : LibmeshPetscCall(VecGetOwnershipRange(petsc_subvector->_vec, &start, nullptr));
1145 1330 : std::iota (idx.begin(), idx.end(), start);
1146 : }
1147 :
1148 : // Construct index sets
1149 76 : WrappedPetsc<IS> parent_is;
1150 1368 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1151 : cast_int<PetscInt>(rows.size()),
1152 : numeric_petsc_cast(rows.data()),
1153 : PETSC_USE_POINTER,
1154 : parent_is.get()));
1155 :
1156 76 : WrappedPetsc<IS> subvector_is;
1157 1368 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1158 : cast_int<PetscInt>(rows.size()),
1159 : idx.data(),
1160 : PETSC_USE_POINTER,
1161 : subvector_is.get()));
1162 :
1163 : // Construct the scatter object
1164 38 : WrappedPetsc<VecScatter> scatter;
1165 1330 : LibmeshPetscCall(VecScatterCreate(this->_vec,
1166 : parent_is,
1167 : petsc_subvector->_vec,
1168 : subvector_is,
1169 : scatter.get()));
1170 :
1171 : // Actually perform the scatter
1172 1330 : VecScatterBeginEnd(this->comm(), scatter, this->_vec, petsc_subvector->_vec, INSERT_VALUES, SCATTER_FORWARD);
1173 :
1174 1330 : petsc_subvector->_is_closed = true;
1175 1330 : }
1176 :
1177 :
1178 :
1179 : template <typename T>
1180 7177316264 : void PetscVector<T>::_get_array(bool read_only) const
1181 : {
1182 570205297 : libmesh_assert (this->initialized());
1183 :
1184 570205297 : bool initially_array_is_present = _array_is_present.load(std::memory_order_acquire);
1185 :
1186 : // If we already have a read/write array - and we're trying
1187 : // to get a read only array - let's just use the read write
1188 7177316264 : if (initially_array_is_present && read_only && !_values_read_only)
1189 24807602 : _read_only_values = _values;
1190 :
1191 : // If the values have already been retrieved and we're currently
1192 : // trying to get a non-read only view (ie read/write) and the
1193 : // values are currently read only... then we need to restore
1194 : // the array first... and then retrieve it again.
1195 7177316264 : if (initially_array_is_present && !read_only && _values_read_only)
1196 : {
1197 210 : _restore_array();
1198 6 : initially_array_is_present = false;
1199 : }
1200 :
1201 7177316060 : if (!initially_array_is_present)
1202 : {
1203 1421676 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1204 1344973 : if (!_array_is_present.load(std::memory_order_relaxed))
1205 : {
1206 1340120 : if (this->type() != GHOSTED)
1207 : {
1208 32761 : if (read_only)
1209 : {
1210 32135 : LibmeshPetscCall(VecGetArrayRead(_vec, &_read_only_values));
1211 32135 : _values_read_only = true;
1212 : }
1213 : else
1214 : {
1215 626 : LibmeshPetscCall(VecGetArray(_vec, &_values));
1216 626 : _values_read_only = false;
1217 : }
1218 32761 : _local_size = this->local_size();
1219 : }
1220 : else
1221 : {
1222 1307359 : LibmeshPetscCall(VecGhostGetLocalForm (_vec,&_local_form));
1223 :
1224 1307359 : if (read_only)
1225 : {
1226 1307359 : LibmeshPetscCall(VecGetArrayRead(_local_form, &_read_only_values));
1227 1307359 : _values_read_only = true;
1228 : }
1229 : else
1230 : {
1231 0 : LibmeshPetscCall(VecGetArray(_local_form, &_values));
1232 0 : _values_read_only = false;
1233 : }
1234 :
1235 1307359 : PetscInt my_local_size = 0;
1236 1307359 : LibmeshPetscCall(VecGetLocalSize(_local_form, &my_local_size));
1237 1307359 : _local_size = static_cast<numeric_index_type>(my_local_size);
1238 : }
1239 :
1240 : { // cache ownership range
1241 1340120 : PetscInt petsc_first=0, petsc_last=0;
1242 1340120 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1243 1340120 : _first = static_cast<numeric_index_type>(petsc_first);
1244 1340120 : _last = static_cast<numeric_index_type>(petsc_last);
1245 : }
1246 74746 : _array_is_present.store(true, std::memory_order_release);
1247 : }
1248 : }
1249 7177316264 : }
1250 :
1251 :
1252 :
1253 : template <typename T>
1254 353813074 : void PetscVector<T>::_restore_array() const
1255 : {
1256 353813074 : libmesh_error_msg_if(_values_manually_retrieved,
1257 : "PetscVector values were manually retrieved but are being automatically restored!");
1258 :
1259 28827970 : libmesh_assert (this->initialized());
1260 353813074 : if (_array_is_present.load(std::memory_order_acquire))
1261 : {
1262 1414866 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1263 1340120 : if (_array_is_present.load(std::memory_order_relaxed))
1264 : {
1265 1340120 : if (this->type() != GHOSTED)
1266 : {
1267 32761 : if (_values_read_only)
1268 32135 : LibmeshPetscCall(VecRestoreArrayRead (_vec, &_read_only_values));
1269 : else
1270 626 : LibmeshPetscCall(VecRestoreArray (_vec, &_values));
1271 :
1272 32761 : _values = nullptr;
1273 : }
1274 : else
1275 : {
1276 1307359 : if (_values_read_only)
1277 1307359 : LibmeshPetscCall(VecRestoreArrayRead (_local_form, &_read_only_values));
1278 : else
1279 0 : LibmeshPetscCall(VecRestoreArray (_local_form, &_values));
1280 :
1281 1307359 : _values = nullptr;
1282 1307359 : LibmeshPetscCall(VecGhostRestoreLocalForm (_vec,&_local_form));
1283 1307359 : _local_form = nullptr;
1284 1307359 : _local_size = 0;
1285 : }
1286 74746 : _array_is_present.store(false, std::memory_order_release);
1287 : }
1288 : }
1289 353813074 : }
1290 :
1291 : template <typename T>
1292 : std::unique_ptr<NumericVector<T>>
1293 8590 : PetscVector<T>::get_subvector(const std::vector<numeric_index_type> & rows)
1294 : {
1295 : // Construct index set
1296 38 : WrappedPetsc<IS> parent_is;
1297 8628 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1298 : cast_int<PetscInt>(rows.size()),
1299 : numeric_petsc_cast(rows.data()),
1300 : PETSC_USE_POINTER,
1301 : parent_is.get()));
1302 :
1303 : Vec subvec;
1304 8590 : LibmeshPetscCall(VecGetSubVector(_vec, parent_is, &subvec));
1305 :
1306 8590 : this->_is_closed = false;
1307 :
1308 8628 : return std::make_unique<PetscVector<T>>(subvec, this->comm());
1309 : }
1310 :
1311 : template <typename T>
1312 : void
1313 8590 : PetscVector<T>::restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
1314 : const std::vector<numeric_index_type> & rows)
1315 : {
1316 38 : auto * const petsc_subvector = cast_ptr<PetscVector<T> *>(subvector.get());
1317 :
1318 : // Construct index set
1319 38 : WrappedPetsc<IS> parent_is;
1320 8628 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1321 : cast_int<PetscInt>(rows.size()),
1322 : numeric_petsc_cast(rows.data()),
1323 : PETSC_USE_POINTER,
1324 : parent_is.get()));
1325 :
1326 8590 : Vec subvec = petsc_subvector->vec();
1327 8590 : LibmeshPetscCall(VecRestoreSubVector(_vec, parent_is, &subvec));
1328 :
1329 8590 : if (this->type() == GHOSTED)
1330 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
1331 :
1332 8590 : this->_is_closed = true;
1333 8590 : }
1334 :
1335 : //------------------------------------------------------------------
1336 : // Explicit instantiations
1337 : template class LIBMESH_EXPORT PetscVector<Number>;
1338 :
1339 : } // namespace libMesh
1340 :
1341 :
1342 :
1343 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
