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 233809 : Real PetscVector<T>::norm () const
61 : {
62 6862 : parallel_object_only();
63 :
64 233809 : this->_restore_array();
65 6862 : libmesh_assert(this->closed());
66 :
67 233809 : PetscReal value=0.;
68 :
69 233809 : LibmeshPetscCall(VecNorm (_vec, N, &value));
70 :
71 233809 : 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 232059 : Real PetscVector<T>::l2_norm () const
80 : {
81 232059 : 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 39082656 : void PetscVector<T>::set (const numeric_index_type i, const T value)
125 : {
126 39082656 : this->_restore_array();
127 3671333 : libmesh_assert_less (i, size());
128 :
129 39082656 : PetscInt i_val = static_cast<PetscInt>(i);
130 39082656 : PetscScalar petsc_value = PS(value);
131 :
132 39082656 : std::scoped_lock lock(this->_numeric_vector_mutex);
133 39082656 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES));
134 :
135 39082656 : this->_is_closed = false;
136 39082656 : }
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 190919454 : void PetscVector<T>::add (const numeric_index_type i, const T value)
166 : {
167 190919454 : this->_restore_array();
168 17728297 : libmesh_assert_less (i, size());
169 :
170 190919454 : PetscInt i_val = static_cast<PetscInt>(i);
171 190919454 : PetscScalar petsc_value = PS(value);
172 :
173 190919454 : std::scoped_lock lock(this->_numeric_vector_mutex);
174 190919454 : LibmeshPetscCall(VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES));
175 :
176 190919454 : this->_is_closed = false;
177 190919454 : }
178 :
179 :
180 :
181 : template <typename T>
182 73914189 : 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 73914189 : if (dof_indices.empty())
187 0 : return;
188 :
189 73914189 : this->_restore_array();
190 :
191 11234253 : const PetscInt * i_val = reinterpret_cast<const PetscInt *>(dof_indices.data());
192 5721574 : const PetscScalar * petsc_value = pPS(v);
193 :
194 73914189 : std::scoped_lock lock(this->_numeric_vector_mutex);
195 79426868 : LibmeshPetscCall(VecSetValues (_vec, cast_int<PetscInt>(dof_indices.size()),
196 : i_val, petsc_value, ADD_VALUES));
197 :
198 73914189 : 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 4361153 : void PetscVector<T>::add (const T a_in, const NumericVector<T> & v_in)
316 : {
317 124770 : parallel_object_only();
318 :
319 4361153 : this->_restore_array();
320 :
321 : // VecAXPY doesn't support &x==&y
322 4361153 : if (this == &v_in)
323 : {
324 140 : this->scale(a_in+1);
325 140 : return;
326 : }
327 :
328 124766 : PetscScalar a = PS(a_in);
329 :
330 : // Make sure the NumericVector passed in is really a PetscVector
331 124766 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
332 4361013 : v->_restore_array();
333 :
334 124766 : libmesh_assert_equal_to (this->size(), v->size());
335 :
336 4361013 : LibmeshPetscCall(VecAXPY(_vec, a, v->vec()));
337 :
338 124766 : libmesh_assert(this->comm().verify(int(this->type())));
339 :
340 4361013 : if (this->type() == GHOSTED)
341 97440 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
342 :
343 4361013 : 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 628490 : void PetscVector<T>::scale (const T factor_in)
369 : {
370 17956 : parallel_object_only();
371 :
372 628490 : this->_restore_array();
373 :
374 17956 : PetscScalar factor = PS(factor_in);
375 :
376 628490 : LibmeshPetscCall(VecScale(_vec, factor));
377 :
378 17956 : libmesh_assert(this->comm().verify(int(this->type())));
379 :
380 628490 : if (this->type() == GHOSTED)
381 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
382 628490 : }
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 3689309 : PetscVector<T>::operator = (const NumericVector<T> & v_in)
497 : {
498 103198 : parallel_object_only();
499 :
500 : // Make sure the NumericVector passed in is really a PetscVector
501 103198 : const PetscVector<T> * v = cast_ptr<const PetscVector<T> *>(&v_in);
502 :
503 3689309 : *this = *v;
504 :
505 3689309 : return *this;
506 : }
507 :
508 :
509 :
510 : template <typename T>
511 : PetscVector<T> &
512 3689379 : PetscVector<T>::operator = (const PetscVector<T> & v)
513 : {
514 103200 : parallel_object_only();
515 :
516 3689379 : this->_restore_array();
517 3689379 : v._restore_array();
518 :
519 103200 : libmesh_assert_equal_to (this->size(), v.size());
520 103200 : libmesh_assert_equal_to (this->local_size(), v.local_size());
521 103200 : libmesh_assert (v.closed());
522 :
523 :
524 3689379 : LibmeshPetscCall(VecCopy (v._vec, this->_vec));
525 :
526 103200 : libmesh_assert(this->comm().verify(int(this->type())));
527 :
528 3689379 : if (this->type() == GHOSTED)
529 2684711 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
530 :
531 3689379 : this->_is_closed = true;
532 :
533 3689379 : 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 2137 : void PetscVector<T>::localize (NumericVector<T> & v_local_in) const
579 : {
580 50 : parallel_object_only();
581 :
582 2137 : 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 2137 : if (v_local->type() == SERIAL && this->size() == v_local->local_size())
599 : {
600 48 : WrappedPetsc<VecScatter> scatter;
601 2067 : LibmeshPetscCall(VecScatterCreateToAll(_vec, scatter.get(), nullptr));
602 2067 : 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 2137 : if (v_local->type() == GHOSTED)
614 70 : VecGhostUpdateBeginEnd(this->comm(), v_local->_vec, INSERT_VALUES, SCATTER_FORWARD);
615 2137 : }
616 :
617 :
618 :
619 : template <typename T>
620 2207330 : void PetscVector<T>::localize (NumericVector<T> & v_local_in,
621 : const std::vector<numeric_index_type> & send_list) const
622 : {
623 60542 : parallel_object_only();
624 :
625 60542 : libmesh_assert(this->comm().verify(int(this->type())));
626 60542 : 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 2267872 : if (v_local_in.type() == GHOSTED &&
632 120080 : this->type() == PARALLEL)
633 : {
634 2075994 : v_local_in = *this;
635 2075994 : 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 27338900 : for (PetscInt i=0; i<n; i++)
865 30179110 : 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 :
1043 8 : WrappedPetsc<PetscViewer> petsc_viewer;
1044 140 : LibmeshPetscCall(PetscViewerCreate (this->comm().get(), petsc_viewer.get()));
1045 :
1046 : // Create an ASCII file containing the matrix
1047 : // if a filename was provided.
1048 140 : if (name != "")
1049 : {
1050 140 : LibmeshPetscCall(PetscViewerASCIIOpen(this->comm().get(),
1051 : name.c_str(),
1052 : petsc_viewer.get()));
1053 :
1054 : #if PETSC_VERSION_LESS_THAN(3,7,0)
1055 : LibmeshPetscCall(PetscViewerSetFormat (petsc_viewer,
1056 : PETSC_VIEWER_ASCII_MATLAB));
1057 : #else
1058 140 : LibmeshPetscCall(PetscViewerPushFormat (petsc_viewer,
1059 : PETSC_VIEWER_ASCII_MATLAB));
1060 : #endif
1061 :
1062 140 : LibmeshPetscCall(VecView (_vec, petsc_viewer));
1063 : }
1064 :
1065 : // Otherwise the matrix will be dumped to the screen.
1066 : else
1067 : {
1068 :
1069 : #if PETSC_VERSION_LESS_THAN(3,7,0)
1070 : LibmeshPetscCall(PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
1071 : PETSC_VIEWER_ASCII_MATLAB));
1072 : #else
1073 0 : LibmeshPetscCall(PetscViewerPushFormat (PETSC_VIEWER_STDOUT_WORLD,
1074 : PETSC_VIEWER_ASCII_MATLAB));
1075 : #endif
1076 :
1077 0 : LibmeshPetscCall(VecView (_vec, PETSC_VIEWER_STDOUT_WORLD));
1078 : }
1079 140 : }
1080 :
1081 :
1082 :
1083 :
1084 :
1085 : template <typename T>
1086 1330 : void PetscVector<T>::create_subvector(NumericVector<T> & subvector,
1087 : const std::vector<numeric_index_type> & rows,
1088 : const bool supplying_global_rows) const
1089 : {
1090 38 : parallel_object_only();
1091 :
1092 1330 : libmesh_error_msg_if(
1093 : subvector.type() == GHOSTED,
1094 : "We do not support scattering parallel information to ghosts for subvectors");
1095 :
1096 1330 : this->_restore_array();
1097 :
1098 : // PETSc data structures
1099 :
1100 : // Make sure the passed in subvector is really a PetscVector
1101 38 : PetscVector<T> * petsc_subvector = cast_ptr<PetscVector<T> *>(&subvector);
1102 :
1103 : // If the petsc_subvector is already initialized, we assume that the
1104 : // user has already allocated the *correct* amount of space for it.
1105 : // If not, we use the appropriate PETSc routines to initialize it.
1106 1330 : if (!petsc_subvector->initialized())
1107 : {
1108 0 : libmesh_assert(petsc_subvector->_type == AUTOMATIC || petsc_subvector->_type == PARALLEL);
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 0 : LibmeshPetscCall(VecSetFromOptions (petsc_subvector->_vec));
1128 :
1129 : // We created a parallel vector
1130 0 : petsc_subvector->_type = PARALLEL;
1131 :
1132 : // Mark the subvector as initialized
1133 0 : petsc_subvector->_is_initialized = true;
1134 : }
1135 : else
1136 : {
1137 1330 : petsc_subvector->_restore_array();
1138 : }
1139 :
1140 1406 : std::vector<PetscInt> idx(rows.size());
1141 1330 : if (supplying_global_rows)
1142 0 : std::iota (idx.begin(), idx.end(), 0);
1143 : else
1144 : {
1145 : PetscInt start;
1146 1330 : LibmeshPetscCall(VecGetOwnershipRange(petsc_subvector->_vec, &start, nullptr));
1147 1330 : std::iota (idx.begin(), idx.end(), start);
1148 : }
1149 :
1150 : // Construct index sets
1151 76 : WrappedPetsc<IS> parent_is;
1152 1368 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1153 : cast_int<PetscInt>(rows.size()),
1154 : numeric_petsc_cast(rows.data()),
1155 : PETSC_USE_POINTER,
1156 : parent_is.get()));
1157 :
1158 76 : WrappedPetsc<IS> subvector_is;
1159 1368 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1160 : cast_int<PetscInt>(rows.size()),
1161 : idx.data(),
1162 : PETSC_USE_POINTER,
1163 : subvector_is.get()));
1164 :
1165 : // Construct the scatter object
1166 38 : WrappedPetsc<VecScatter> scatter;
1167 1330 : LibmeshPetscCall(VecScatterCreate(this->_vec,
1168 : parent_is,
1169 : petsc_subvector->_vec,
1170 : subvector_is,
1171 : scatter.get()));
1172 :
1173 : // Actually perform the scatter
1174 1330 : VecScatterBeginEnd(this->comm(), scatter, this->_vec, petsc_subvector->_vec, INSERT_VALUES, SCATTER_FORWARD);
1175 :
1176 1330 : petsc_subvector->_is_closed = true;
1177 1330 : }
1178 :
1179 :
1180 :
1181 : template <typename T>
1182 7175880747 : void PetscVector<T>::_get_array(bool read_only) const
1183 : {
1184 570195256 : libmesh_assert (this->initialized());
1185 :
1186 570195256 : bool initially_array_is_present = _array_is_present.load(std::memory_order_acquire);
1187 :
1188 : // If we already have a read/write array - and we're trying
1189 : // to get a read only array - let's just use the read write
1190 7175880747 : if (initially_array_is_present && read_only && !_values_read_only)
1191 24807602 : _read_only_values = _values;
1192 :
1193 : // If the values have already been retrieved and we're currently
1194 : // trying to get a non-read only view (ie read/write) and the
1195 : // values are currently read only... then we need to restore
1196 : // the array first... and then retrieve it again.
1197 7175880747 : if (initially_array_is_present && !read_only && _values_read_only)
1198 : {
1199 210 : _restore_array();
1200 6 : initially_array_is_present = false;
1201 : }
1202 :
1203 7175880543 : if (!initially_array_is_present)
1204 : {
1205 1399973 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1206 1324405 : if (!_array_is_present.load(std::memory_order_relaxed))
1207 : {
1208 1319352 : if (this->type() != GHOSTED)
1209 : {
1210 32748 : if (read_only)
1211 : {
1212 32122 : LibmeshPetscCall(VecGetArrayRead(_vec, &_read_only_values));
1213 32122 : _values_read_only = true;
1214 : }
1215 : else
1216 : {
1217 626 : LibmeshPetscCall(VecGetArray(_vec, &_values));
1218 626 : _values_read_only = false;
1219 : }
1220 32748 : _local_size = this->local_size();
1221 : }
1222 : else
1223 : {
1224 1286604 : LibmeshPetscCall(VecGhostGetLocalForm (_vec,&_local_form));
1225 :
1226 1286604 : if (read_only)
1227 : {
1228 1286604 : LibmeshPetscCall(VecGetArrayRead(_local_form, &_read_only_values));
1229 1286604 : _values_read_only = true;
1230 : }
1231 : else
1232 : {
1233 0 : LibmeshPetscCall(VecGetArray(_local_form, &_values));
1234 0 : _values_read_only = false;
1235 : }
1236 :
1237 1286604 : PetscInt my_local_size = 0;
1238 1286604 : LibmeshPetscCall(VecGetLocalSize(_local_form, &my_local_size));
1239 1286604 : _local_size = static_cast<numeric_index_type>(my_local_size);
1240 : }
1241 :
1242 : { // cache ownership range
1243 1319352 : PetscInt petsc_first=0, petsc_last=0;
1244 1319352 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1245 1319352 : _first = static_cast<numeric_index_type>(petsc_first);
1246 1319352 : _last = static_cast<numeric_index_type>(petsc_last);
1247 : }
1248 73910 : _array_is_present.store(true, std::memory_order_release);
1249 : }
1250 : }
1251 7175880747 : }
1252 :
1253 :
1254 :
1255 : template <typename T>
1256 353609484 : void PetscVector<T>::_restore_array() const
1257 : {
1258 353609484 : libmesh_error_msg_if(_values_manually_retrieved,
1259 : "PetscVector values were manually retrieved but are being automatically restored!");
1260 :
1261 28826540 : libmesh_assert (this->initialized());
1262 353609484 : if (_array_is_present.load(std::memory_order_acquire))
1263 : {
1264 1393262 : std::scoped_lock lock(_petsc_get_restore_array_mutex);
1265 1319352 : if (_array_is_present.load(std::memory_order_relaxed))
1266 : {
1267 1319352 : if (this->type() != GHOSTED)
1268 : {
1269 32748 : if (_values_read_only)
1270 32122 : LibmeshPetscCall(VecRestoreArrayRead (_vec, &_read_only_values));
1271 : else
1272 626 : LibmeshPetscCall(VecRestoreArray (_vec, &_values));
1273 :
1274 32748 : _values = nullptr;
1275 : }
1276 : else
1277 : {
1278 1286604 : if (_values_read_only)
1279 1286604 : LibmeshPetscCall(VecRestoreArrayRead (_local_form, &_read_only_values));
1280 : else
1281 0 : LibmeshPetscCall(VecRestoreArray (_local_form, &_values));
1282 :
1283 1286604 : _values = nullptr;
1284 1286604 : LibmeshPetscCall(VecGhostRestoreLocalForm (_vec,&_local_form));
1285 1286604 : _local_form = nullptr;
1286 1286604 : _local_size = 0;
1287 : }
1288 73910 : _array_is_present.store(false, std::memory_order_release);
1289 : }
1290 : }
1291 353609484 : }
1292 :
1293 : template <typename T>
1294 : std::unique_ptr<NumericVector<T>>
1295 8590 : PetscVector<T>::get_subvector(const std::vector<numeric_index_type> & rows)
1296 : {
1297 : // Construct index set
1298 38 : WrappedPetsc<IS> parent_is;
1299 8628 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1300 : cast_int<PetscInt>(rows.size()),
1301 : numeric_petsc_cast(rows.data()),
1302 : PETSC_USE_POINTER,
1303 : parent_is.get()));
1304 :
1305 : Vec subvec;
1306 8590 : LibmeshPetscCall(VecGetSubVector(_vec, parent_is, &subvec));
1307 :
1308 8590 : this->_is_closed = false;
1309 :
1310 8628 : return std::make_unique<PetscVector<T>>(subvec, this->comm());
1311 : }
1312 :
1313 : template <typename T>
1314 : void
1315 8590 : PetscVector<T>::restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
1316 : const std::vector<numeric_index_type> & rows)
1317 : {
1318 38 : auto * const petsc_subvector = cast_ptr<PetscVector<T> *>(subvector.get());
1319 :
1320 : // Construct index set
1321 38 : WrappedPetsc<IS> parent_is;
1322 8628 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
1323 : cast_int<PetscInt>(rows.size()),
1324 : numeric_petsc_cast(rows.data()),
1325 : PETSC_USE_POINTER,
1326 : parent_is.get()));
1327 :
1328 8590 : Vec subvec = petsc_subvector->vec();
1329 8590 : LibmeshPetscCall(VecRestoreSubVector(_vec, parent_is, &subvec));
1330 :
1331 8590 : if (this->type() == GHOSTED)
1332 0 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
1333 :
1334 8590 : this->_is_closed = true;
1335 8590 : }
1336 :
1337 : //------------------------------------------------------------------
1338 : // Explicit instantiations
1339 : template class LIBMESH_EXPORT PetscVector<Number>;
1340 :
1341 : } // namespace libMesh
1342 :
1343 :
1344 :
1345 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
