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 :
19 :
20 :
21 : #ifndef LIBMESH_PETSC_VECTOR_H
22 : #define LIBMESH_PETSC_VECTOR_H
23 :
24 :
25 : #include "libmesh/libmesh_config.h"
26 :
27 : #ifdef LIBMESH_HAVE_PETSC
28 :
29 : // Local includes
30 : #include "libmesh/numeric_vector.h"
31 : #include "libmesh/petsc_macro.h"
32 : #include "libmesh/int_range.h"
33 : #include "libmesh/libmesh_common.h"
34 : #include "libmesh/petsc_solver_exception.h"
35 : #include "libmesh/parallel_only.h"
36 : #include "libmesh/enum_to_string.h"
37 :
38 : // PETSc include files.
39 : #ifdef I
40 : # define LIBMESH_SAW_I
41 : #endif
42 : #include <petscvec.h>
43 : #ifndef LIBMESH_SAW_I
44 : # undef I // Avoid complex.h contamination
45 : #endif
46 :
47 : // C++ includes
48 : #include <cstddef>
49 : #include <cstring>
50 : #include <vector>
51 : #include <unordered_map>
52 : #include <limits>
53 :
54 : #include <atomic>
55 : #include <mutex>
56 : #include <condition_variable>
57 :
58 : namespace libMesh
59 : {
60 :
61 : // forward declarations
62 : template <typename T> class SparseMatrix;
63 :
64 : /**
65 : * This class provides a nice interface to PETSc's Vec object. All
66 : * overridden virtual functions are documented in numeric_vector.h.
67 : *
68 : * \author Benjamin S. Kirk
69 : * \date 2002
70 : * \brief NumericVector interface to PETSc Vec.
71 : */
72 : template <typename T>
73 : class PetscVector final : public NumericVector<T>
74 : {
75 : public:
76 :
77 : /**
78 : * Dummy-Constructor. Dimension=0
79 : */
80 : explicit
81 : PetscVector (const Parallel::Communicator & comm_in,
82 : const ParallelType type = AUTOMATIC);
83 :
84 : /**
85 : * Constructor. Set dimension to \p n and initialize all elements with zero.
86 : */
87 : explicit
88 : PetscVector (const Parallel::Communicator & comm_in,
89 : const numeric_index_type n,
90 : const ParallelType type = AUTOMATIC);
91 :
92 : /**
93 : * Constructor. Set local dimension to \p n_local, the global dimension
94 : * to \p n, and initialize all elements with zero.
95 : */
96 : PetscVector (const Parallel::Communicator & comm_in,
97 : const numeric_index_type n,
98 : const numeric_index_type n_local,
99 : const ParallelType type = AUTOMATIC);
100 :
101 : /**
102 : * Constructor. Set local dimension to \p n_local, the global
103 : * dimension to \p n, but additionally reserve memory for the
104 : * indices specified by the \p ghost argument.
105 : */
106 : PetscVector (const Parallel::Communicator & comm_in,
107 : const numeric_index_type N,
108 : const numeric_index_type n_local,
109 : const std::vector<numeric_index_type> & ghost,
110 : const ParallelType type = AUTOMATIC);
111 :
112 : /**
113 : * Constructor. Creates a PetscVector assuming you already have a
114 : * valid PETSc Vec object. In this case, \p v is NOT destroyed by the
115 : * PetscVector constructor when this object goes out of scope.
116 : * This allows ownership of \p v to remain with the original creator,
117 : * and to simply provide additional functionality with the PetscVector.
118 : */
119 : PetscVector(Vec v,
120 : const Parallel::Communicator & comm_in);
121 :
122 : /**
123 : * Copy assignment operator.
124 : * Calls VecCopy after performing various checks.
125 : * \returns A reference to *this as the derived type.
126 : */
127 : PetscVector<T> & operator= (const PetscVector<T> & v);
128 :
129 : /**
130 : * This class manages a C-style struct (Vec) manually, so we
131 : * don't want to allow any automatic copy/move functions to be
132 : * generated, and we can't default the destructor.
133 : */
134 : PetscVector (PetscVector &&) = delete;
135 : PetscVector (const PetscVector &) = delete;
136 : PetscVector & operator= (PetscVector &&) = delete;
137 : virtual ~PetscVector ();
138 :
139 : virtual void close () override;
140 :
141 : /**
142 : * clear() is called from the destructor, so it should not throw.
143 : */
144 : virtual void clear () noexcept override;
145 :
146 : virtual void zero () override;
147 :
148 : virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
149 :
150 : virtual std::unique_ptr<NumericVector<T>> clone () const override;
151 :
152 : virtual void init (const numeric_index_type N,
153 : const numeric_index_type n_local,
154 : const bool fast=false,
155 : const ParallelType type=AUTOMATIC) override;
156 :
157 : virtual void init (const numeric_index_type N,
158 : const bool fast=false,
159 : const ParallelType type=AUTOMATIC) override;
160 :
161 : virtual void init (const numeric_index_type N,
162 : const numeric_index_type n_local,
163 : const std::vector<numeric_index_type> & ghost,
164 : const bool fast = false,
165 : const ParallelType = AUTOMATIC) override;
166 :
167 : virtual void init (const NumericVector<T> & other,
168 : const bool fast = false) override;
169 :
170 : virtual NumericVector<T> & operator= (const T s) override;
171 :
172 : virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
173 :
174 : virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
175 :
176 : virtual Real min () const override;
177 :
178 : virtual Real max () const override;
179 :
180 : virtual T sum () const override;
181 :
182 : virtual Real l1_norm () const override;
183 :
184 : virtual Real l2_norm () const override;
185 :
186 : virtual Real linfty_norm () const override;
187 :
188 : virtual numeric_index_type size () const override;
189 :
190 : virtual numeric_index_type local_size() const override;
191 :
192 : virtual numeric_index_type first_local_index() const override;
193 :
194 : virtual numeric_index_type last_local_index() const override;
195 :
196 : /**
197 : * \returns The local index corresponding to global index \p i.
198 : *
199 : * If the index is not a ghost cell, this is done by subtraction the
200 : * number of the first local index. If it is a ghost cell, it has
201 : * to be looked up in the map.
202 : */
203 : numeric_index_type map_global_to_local_index(const numeric_index_type i) const;
204 :
205 : virtual T operator() (const numeric_index_type i) const override;
206 :
207 : virtual void get(const std::vector<numeric_index_type> & index,
208 : T * values) const override;
209 :
210 : /**
211 : * Get read/write access to the raw PETSc Vector data array.
212 : *
213 : * \note This is an advanced interface. In general you should avoid
214 : * using it unless you know what you are doing!
215 : */
216 : PetscScalar * get_array();
217 :
218 : /**
219 : * Get read only access to the raw PETSc Vector data array.
220 : *
221 : * \note This is an advanced interface. In general you should avoid
222 : * using it unless you know what you are doing!
223 : */
224 : const PetscScalar * get_array_read() const;
225 :
226 : /**
227 : * Restore the data array.
228 : *
229 : * \note This MUST be called after get_array() or get_array_read()
230 : * and before using any other interface functions on PetscVector.
231 : */
232 : void restore_array();
233 :
234 : virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
235 :
236 : virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
237 :
238 : virtual void reciprocal() override;
239 :
240 : virtual void conjugate() override;
241 :
242 : virtual void set (const numeric_index_type i,
243 : const T value) override;
244 :
245 : virtual void add (const numeric_index_type i,
246 : const T value) override;
247 :
248 : virtual void add (const T s) override;
249 :
250 : virtual void add (const NumericVector<T> & v) override;
251 :
252 : virtual void add (const T a, const NumericVector<T> & v) override;
253 :
254 : /**
255 : * We override two NumericVector<T>::add_vector() methods but don't
256 : * want to hide the other defaults.
257 : */
258 : using NumericVector<T>::add_vector;
259 :
260 : virtual void add_vector (const T * v,
261 : const std::vector<numeric_index_type> & dof_indices) override;
262 :
263 : virtual void add_vector (const NumericVector<T> & v,
264 : const SparseMatrix<T> & A) override;
265 :
266 : virtual void add_vector_transpose (const NumericVector<T> & v,
267 : const SparseMatrix<T> & A) override;
268 :
269 : /**
270 : * \f$ U \leftarrow U + A^H v \f$.
271 : *
272 : * Adds the product of the conjugate-transpose of \p SparseMatrix \p
273 : * A and a \p NumericVector \p v to \p this.
274 : */
275 : void add_vector_conjugate_transpose (const NumericVector<T> & v,
276 : const SparseMatrix<T> & A);
277 :
278 : /**
279 : * We override one NumericVector<T>::insert() method but don't want
280 : * to hide the other defaults
281 : */
282 : using NumericVector<T>::insert;
283 :
284 : virtual void insert (const T * v,
285 : const std::vector<numeric_index_type> & dof_indices) override;
286 :
287 : virtual void scale (const T factor) override;
288 :
289 : virtual NumericVector<T> & operator *= (const NumericVector<T> & v) override;
290 :
291 : virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
292 :
293 : virtual void abs() override;
294 :
295 : virtual T dot(const NumericVector<T> & v) const override;
296 :
297 : /**
298 : * \returns The dot product of (*this) with the vector \p v.
299 : *
300 : * \note Does *not* use the complex-conjugate of v in the complex-valued case.
301 : */
302 : T indefinite_dot(const NumericVector<T> & v) const;
303 :
304 : virtual void localize (std::vector<T> & v_local) const override;
305 :
306 : virtual void localize (NumericVector<T> & v_local) const override;
307 :
308 : virtual void localize (NumericVector<T> & v_local,
309 : const std::vector<numeric_index_type> & send_list) const override;
310 :
311 : virtual void localize (std::vector<T> & v_local,
312 : const std::vector<numeric_index_type> & indices) const override;
313 :
314 : virtual void localize (const numeric_index_type first_local_idx,
315 : const numeric_index_type last_local_idx,
316 : const std::vector<numeric_index_type> & send_list) override;
317 :
318 : virtual void localize_to_one (std::vector<T> & v_local,
319 : const processor_id_type proc_id=0) const override;
320 :
321 : virtual void pointwise_mult (const NumericVector<T> & vec1,
322 : const NumericVector<T> & vec2) override;
323 :
324 : virtual void pointwise_divide (const NumericVector<T> & vec1,
325 : const NumericVector<T> & vec2) override;
326 :
327 : virtual void print_matlab(const std::string & name = "") const override;
328 :
329 : virtual void create_subvector(NumericVector<T> & subvector,
330 : const std::vector<numeric_index_type> & rows,
331 : bool supplying_global_rows = true) const override;
332 :
333 : virtual void swap (NumericVector<T> & v) override;
334 :
335 : virtual std::size_t max_allowed_id() const override;
336 :
337 : /**
338 : * \returns The raw PETSc Vec pointer.
339 : *
340 : * \note This is generally not required in user-level code.
341 : *
342 : * \note Don't do anything crazy like calling VecDestroy() on
343 : * it, or very bad things will likely happen!
344 : */
345 544429 : Vec vec () { libmesh_assert (_vec); return _vec; }
346 :
347 4361157 : Vec vec () const { libmesh_assert (_vec); return _vec; }
348 :
349 : virtual std::unique_ptr<NumericVector<T>>
350 : get_subvector(const std::vector<numeric_index_type> & rows) override;
351 :
352 : virtual void restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
353 : const std::vector<numeric_index_type> & rows) override;
354 :
355 : private:
356 :
357 : /**
358 : * Actual PETSc vector datatype to hold vector entries.
359 : */
360 : Vec _vec;
361 :
362 : /**
363 : * If \p true, the actual PETSc array of the values of the vector is
364 : * currently accessible. That means that the members \p _local_form
365 : * and \p _values are valid.
366 : */
367 : #ifdef LIBMESH_HAVE_CXX11_THREAD
368 : // We can't use std::atomic_flag here because we need load and store operations.
369 : mutable std::atomic<bool> _array_is_present;
370 : #else
371 : mutable bool _array_is_present;
372 : #endif
373 :
374 : /**
375 : * First local index.
376 : *
377 : * Only valid when _array_is_present
378 : */
379 : mutable numeric_index_type _first;
380 :
381 : /**
382 : * Last local index.
383 : *
384 : * Only valid when _array_is_present
385 : */
386 : mutable numeric_index_type _last;
387 :
388 : /**
389 : * Size of the local values from _get_array()
390 : */
391 : mutable numeric_index_type _local_size;
392 :
393 : /**
394 : * PETSc vector datatype to hold the local form of a ghosted vector.
395 : * The contents of this field are only valid if the vector is
396 : * ghosted and \p _array_is_present is \p true.
397 : */
398 : mutable Vec _local_form;
399 :
400 : /**
401 : * Pointer to the actual PETSc array of the values of the vector.
402 : * This pointer is only valid if \p _array_is_present is \p true.
403 : * We're using PETSc's VecGetArrayRead() function, which requires a
404 : * constant PetscScalar *, but _get_array and _restore_array are
405 : * const member functions, so _values also needs to be mutable
406 : * (otherwise it is a "const PetscScalar * const" in that context).
407 : */
408 : mutable const PetscScalar * _read_only_values;
409 :
410 : /**
411 : * Pointer to the actual PETSc array of the values of the vector.
412 : * This pointer is only valid if \p _array_is_present is \p true.
413 : * We're using PETSc's VecGetArrayRead() function, which requires a
414 : * constant PetscScalar *, but _get_array and _restore_array are
415 : * const member functions, so _values also needs to be mutable
416 : * (otherwise it is a "const PetscScalar * const" in that context).
417 : */
418 : mutable PetscScalar * _values;
419 :
420 : /**
421 : * Mutex for _get_array and _restore_array. This is part of the
422 : * object to keep down thread contention when reading from multiple
423 : * PetscVectors simultaneously
424 : */
425 : mutable std::mutex _petsc_get_restore_array_mutex;
426 :
427 : /**
428 : * Queries the array (and the local form if the vector is ghosted)
429 : * from PETSc.
430 : *
431 : * \param read_only Whether or not a read only copy of the raw data
432 : */
433 : void _get_array(bool read_only) const;
434 :
435 : /**
436 : * Restores the array (and the local form if the vector is ghosted)
437 : * to PETSc.
438 : */
439 : void _restore_array() const;
440 :
441 : /**
442 : * Type for map that maps global to local ghost cells.
443 : */
444 : typedef std::unordered_map<numeric_index_type,numeric_index_type> GlobalToLocalMap;
445 :
446 : /**
447 : * Map that maps global to local ghost cells (will be empty if not
448 : * in ghost cell mode).
449 : */
450 : GlobalToLocalMap _global_to_local_map;
451 :
452 : /**
453 : * This boolean value should only be set to false
454 : * for the constructor which takes a PETSc Vec object.
455 : */
456 : bool _destroy_vec_on_exit;
457 :
458 : /**
459 : * Whether or not the data array has been manually retrieved using
460 : * get_array() or get_array_read()
461 : */
462 : mutable bool _values_manually_retrieved;
463 :
464 : /**
465 : * Whether or not the data array is for read only access
466 : */
467 : mutable bool _values_read_only;
468 :
469 : /**
470 : * \returns A norm of the vector, where the type of norm to compute is
471 : * determined by the template parameter N of the PETSc-defined type NormType.
472 : * The valid template arguments are NORM_1, NORM_2 and NORM_INFINITY, as used
473 : * to define l1_norm(), l2_norm() and linfty_norm().
474 : */
475 : template <NormType N> Real norm () const;
476 : };
477 :
478 :
479 : /*----------------------- Inline functions ----------------------------------*/
480 :
481 :
482 :
483 : template <typename T>
484 : inline
485 2213769 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in, const ParallelType ptype) :
486 : NumericVector<T>(comm_in, ptype),
487 2094879 : _vec(nullptr),
488 : _array_is_present(false),
489 2094879 : _first(0),
490 2094879 : _last(0),
491 2094879 : _local_size(0),
492 2094879 : _local_form(nullptr),
493 2094879 : _read_only_values(nullptr),
494 2094879 : _values(nullptr),
495 2094879 : _global_to_local_map(),
496 2094879 : _destroy_vec_on_exit(true),
497 2094879 : _values_manually_retrieved(false),
498 2213769 : _values_read_only(false)
499 : {
500 2147831 : this->_type = ptype;
501 2213769 : }
502 :
503 :
504 :
505 : template <typename T>
506 : inline
507 0 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
508 : const numeric_index_type n,
509 : const ParallelType ptype) :
510 : NumericVector<T>(comm_in, ptype),
511 0 : _vec(nullptr),
512 : _array_is_present(false),
513 0 : _first(0),
514 0 : _last(0),
515 0 : _local_size(0),
516 0 : _local_form(nullptr),
517 0 : _read_only_values(nullptr),
518 0 : _values(nullptr),
519 0 : _global_to_local_map(),
520 0 : _destroy_vec_on_exit(true),
521 0 : _values_manually_retrieved(false),
522 0 : _values_read_only(false)
523 : {
524 0 : this->init(n, n, false, ptype);
525 0 : }
526 :
527 :
528 :
529 : template <typename T>
530 : inline
531 0 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
532 : const numeric_index_type n,
533 : const numeric_index_type n_local,
534 : const ParallelType ptype) :
535 : NumericVector<T>(comm_in, ptype),
536 0 : _vec(nullptr),
537 : _array_is_present(false),
538 0 : _first(0),
539 0 : _last(0),
540 0 : _local_size(0),
541 0 : _local_form(nullptr),
542 0 : _read_only_values(nullptr),
543 0 : _values(nullptr),
544 0 : _global_to_local_map(),
545 0 : _destroy_vec_on_exit(true),
546 0 : _values_manually_retrieved(false),
547 0 : _values_read_only(false)
548 : {
549 0 : this->init(n, n_local, false, ptype);
550 0 : }
551 :
552 :
553 :
554 : template <typename T>
555 : inline
556 0 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
557 : const numeric_index_type n,
558 : const numeric_index_type n_local,
559 : const std::vector<numeric_index_type> & ghost,
560 : const ParallelType ptype) :
561 : NumericVector<T>(comm_in, ptype),
562 0 : _vec(nullptr),
563 : _array_is_present(false),
564 0 : _first(0),
565 0 : _last(0),
566 0 : _local_size(0),
567 0 : _local_form(nullptr),
568 0 : _read_only_values(nullptr),
569 0 : _values(nullptr),
570 0 : _global_to_local_map(),
571 0 : _destroy_vec_on_exit(true),
572 0 : _values_manually_retrieved(false),
573 0 : _values_read_only(false)
574 : {
575 0 : this->init(n, n_local, ghost, false, ptype);
576 0 : }
577 :
578 :
579 :
580 :
581 :
582 : template <typename T>
583 : inline
584 989439 : PetscVector<T>::PetscVector (Vec v,
585 : const Parallel::Communicator & comm_in) :
586 : NumericVector<T>(comm_in, AUTOMATIC),
587 948479 : _vec(v),
588 : _array_is_present(false),
589 948479 : _first(0),
590 948479 : _last(0),
591 948479 : _local_size(0),
592 948479 : _local_form(nullptr),
593 948479 : _read_only_values(nullptr),
594 948479 : _values(nullptr),
595 948479 : _global_to_local_map(),
596 948479 : _destroy_vec_on_exit(false),
597 948479 : _values_manually_retrieved(false),
598 989439 : _values_read_only(false)
599 : {
600 989439 : this->_is_closed = true;
601 989439 : this->_is_initialized = true;
602 :
603 : /* We need to ask PETSc about the (local to global) ghost value
604 : mapping and create the inverse mapping out of it. */
605 989439 : PetscInt petsc_local_size=0;
606 989439 : LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_local_size));
607 :
608 : // Get the vector type from PETSc.
609 : VecType ptype;
610 989439 : LibmeshPetscCall(VecGetType(_vec, &ptype));
611 :
612 : #if PETSC_RELEASE_GREATER_EQUALS(3, 21, 0)
613 : // Fande only implemented VecGhostGetGhostIS for VECMPI
614 : if (std::strcmp(ptype, VECMPI) == 0)
615 : {
616 : IS ghostis;
617 : LibmeshPetscCall(VecGhostGetGhostIS(_vec, &ghostis));
618 :
619 : Vec localrep;
620 : LibmeshPetscCall(VecGhostGetLocalForm(_vec, &localrep));
621 :
622 : // If is a sparsely stored vector, set up our new mapping
623 : // Only checking mapping is not enough to determine if a vec is ghosted
624 : // We need to check if vec has a local representation
625 : if (ghostis && localrep)
626 : {
627 : PetscInt ghost_size;
628 : LibmeshPetscCall(ISGetSize(ghostis, &ghost_size));
629 :
630 : const PetscInt * indices;
631 : LibmeshPetscCall(ISGetIndices(ghostis, &indices));
632 :
633 : for (const auto i : make_range(ghost_size))
634 : _global_to_local_map[indices[i]] = i;
635 : this->_type = GHOSTED;
636 : LibmeshPetscCall(ISRestoreIndices(ghostis, &indices));
637 : }
638 : }
639 : else if (std::strcmp(ptype,VECSHARED) == 0)
640 : #else
641 989439 : if ((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
642 : #endif
643 : {
644 : ISLocalToGlobalMapping mapping;
645 974895 : LibmeshPetscCall(VecGetLocalToGlobalMapping(_vec, &mapping));
646 :
647 : Vec localrep;
648 974895 : LibmeshPetscCall(VecGhostGetLocalForm(_vec,&localrep));
649 : // If is a sparsely stored vector, set up our new mapping
650 : // Only checking mapping is not enough to determine if a vec is ghosted
651 : // We need to check if vec has a local representation
652 974895 : if (mapping && localrep)
653 : {
654 0 : const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
655 0 : const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
656 : PetscInt n;
657 0 : LibmeshPetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
658 :
659 0 : const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
660 : const PetscInt * indices;
661 0 : LibmeshPetscCall(ISLocalToGlobalMappingGetIndices(mapping,&indices));
662 :
663 0 : for (numeric_index_type i=ghost_begin; i<ghost_end; i++)
664 0 : _global_to_local_map[indices[i]] = i-my_local_size;
665 0 : this->_type = GHOSTED;
666 0 : LibmeshPetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
667 0 : }
668 : else
669 974895 : this->_type = PARALLEL;
670 :
671 974895 : LibmeshPetscCall(VecGhostRestoreLocalForm(_vec,&localrep));
672 933935 : }
673 : else
674 14544 : this->_type = SERIAL;
675 989439 : }
676 :
677 :
678 :
679 :
680 : template <typename T>
681 : inline
682 4506103 : PetscVector<T>::~PetscVector ()
683 : {
684 2751117 : this->clear ();
685 4578803 : }
686 :
687 :
688 :
689 : template <typename T>
690 : inline
691 914472 : void PetscVector<T>::init (const numeric_index_type n,
692 : const numeric_index_type n_local,
693 : const bool fast,
694 : const ParallelType ptype)
695 : {
696 27300 : parallel_object_only();
697 :
698 914472 : PetscInt petsc_n=static_cast<PetscInt>(n);
699 :
700 : // Clear initialized vectors
701 914472 : if (this->initialized())
702 26764 : this->clear();
703 :
704 914472 : if (ptype == AUTOMATIC)
705 : {
706 1190 : if (n == n_local)
707 569 : this->_type = SERIAL;
708 : else
709 621 : this->_type = PARALLEL;
710 : }
711 : else
712 913282 : this->_type = ptype;
713 :
714 27300 : libmesh_assert ((this->_type==SERIAL && n==n_local) ||
715 : this->_type==PARALLEL);
716 :
717 : // create a sequential vector if on only 1 processor
718 914472 : if (this->_type == SERIAL)
719 : {
720 2146 : LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
721 2146 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
722 2146 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetFromOptions (_vec));
723 : }
724 : // otherwise create an MPI-enabled vector
725 912326 : else if (this->_type == PARALLEL)
726 : {
727 : #ifdef LIBMESH_HAVE_MPI
728 27250 : PetscInt petsc_n_local=cast_int<PetscInt>(n_local);
729 27250 : libmesh_assert_less_equal (n_local, n);
730 : // Use more generic function instead of VecCreateSeq/MPI
731 912326 : LibmeshPetscCall(VecCreate(this->comm().get(), &_vec));
732 912326 : LibmeshPetscCall(VecSetSizes(_vec, petsc_n_local, petsc_n));
733 : #else
734 : libmesh_assert_equal_to (n_local, n);
735 : LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
736 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
737 : #endif
738 912326 : LibmeshPetscCall(VecSetFromOptions(_vec));
739 : }
740 : else
741 0 : libmesh_error_msg("Unsupported type " << Utility::enum_to_string(this->_type));
742 :
743 914472 : this->_is_initialized = true;
744 914472 : this->_is_closed = true;
745 :
746 :
747 914472 : if (fast == false)
748 905624 : this->zero ();
749 914472 : }
750 :
751 :
752 :
753 : template <typename T>
754 : inline
755 1593 : void PetscVector<T>::init (const numeric_index_type n,
756 : const bool fast,
757 : const ParallelType ptype)
758 : {
759 2137 : this->init(n,n,fast,ptype);
760 2121 : }
761 :
762 :
763 :
764 : template <typename T>
765 : inline
766 508076 : void PetscVector<T>::init (const numeric_index_type n,
767 : const numeric_index_type n_local,
768 : const std::vector<numeric_index_type> & ghost,
769 : const bool fast,
770 : const ParallelType libmesh_dbg_var(ptype))
771 : {
772 16244 : parallel_object_only();
773 :
774 508076 : PetscInt petsc_n=static_cast<PetscInt>(n);
775 508076 : PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
776 32300 : PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
777 :
778 : // If the mesh is not disjoint, every processor will either have
779 : // all the dofs, none of the dofs, or some non-zero dofs at the
780 : // boundary between processors.
781 : //
782 : // However we can't assert this, because someone might want to
783 : // construct a GHOSTED vector which doesn't include neighbor element
784 : // dofs. Boyce tried to do so in user code, and we're going to want
785 : // to do so in System::project_vector().
786 : //
787 : // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
788 :
789 : libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
790 :
791 508076 : PetscInt * petsc_ghost = ghost.empty() ? LIBMESH_PETSC_NULLPTR :
792 15855 : const_cast<PetscInt *>(reinterpret_cast<const PetscInt *>(ghost.data()));
793 :
794 : // Clear initialized vectors
795 508076 : if (this->initialized())
796 20692 : this->clear();
797 :
798 16244 : libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
799 508076 : this->_type = GHOSTED;
800 :
801 : /* Make the global-to-local ghost cell map. */
802 35174025 : for (auto i : index_range(ghost))
803 : {
804 34665949 : _global_to_local_map[ghost[i]] = i;
805 : }
806 :
807 : /* Create vector. */
808 508076 : LibmeshPetscCall(VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
809 : petsc_n_ghost, petsc_ghost,
810 : &_vec));
811 :
812 : // Add a prefix so that we can at least distinguish a ghosted vector from a
813 : // nonghosted vector when using a petsc option.
814 : // PETSc does not fully support VecGhost on GPU yet. This change allows us to
815 : // trigger a nonghosted vector to use GPU without bothering the ghosted vectors.
816 508076 : LibmeshPetscCall(PetscObjectAppendOptionsPrefix((PetscObject)_vec,"ghost_"));
817 :
818 508076 : LibmeshPetscCall(VecSetFromOptions (_vec));
819 :
820 508076 : this->_is_initialized = true;
821 508076 : this->_is_closed = true;
822 :
823 508076 : if (fast == false)
824 157395 : this->zero ();
825 508076 : }
826 :
827 :
828 :
829 : template <typename T>
830 : inline
831 424818 : void PetscVector<T>::init (const NumericVector<T> & other,
832 : const bool fast)
833 : {
834 12538 : parallel_object_only();
835 :
836 : // Clear initialized vectors
837 424818 : if (this->initialized())
838 32690 : this->clear();
839 :
840 12538 : const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
841 :
842 : // Other vector should restore array.
843 424818 : if (v.initialized())
844 : {
845 424818 : v._restore_array();
846 : }
847 :
848 12538 : this->_global_to_local_map = v._global_to_local_map;
849 :
850 : // Even if we're initializing sizes based on an uninitialized or
851 : // unclosed vector, *this* vector is being initialized now and is
852 : // initially closed.
853 424818 : this->_is_closed = true; // v._is_closed;
854 424818 : this->_is_initialized = true; // v._is_initialized;
855 :
856 424818 : this->_type = v._type;
857 :
858 : // We want to have a valid Vec, even if it's initially of size zero
859 424818 : LibmeshPetscCall(VecDuplicate (v._vec, &this->_vec));
860 :
861 424818 : if (fast == false)
862 53720 : this->zero ();
863 424818 : }
864 :
865 :
866 :
867 : template <typename T>
868 : inline
869 3557038 : void PetscVector<T>::close ()
870 : {
871 97244 : parallel_object_only();
872 :
873 3557038 : this->_restore_array();
874 :
875 3557038 : VecAssemblyBeginEnd(this->comm(), _vec);
876 :
877 3557038 : if (this->type() == GHOSTED)
878 898608 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
879 :
880 3557038 : this->_is_closed = true;
881 3557038 : }
882 :
883 :
884 :
885 : template <typename T>
886 : inline
887 2997195 : void PetscVector<T>::clear () noexcept
888 : {
889 81766 : exceptionless_parallel_object_only();
890 :
891 2997195 : if (this->initialized())
892 2836805 : this->_restore_array();
893 :
894 2997195 : if ((this->initialized()) && (this->_destroy_vec_on_exit))
895 : {
896 : // If we encounter an error here, print a warning but otherwise
897 : // keep going since we may be recovering from an exception.
898 1847366 : PetscErrorCode ierr = VecDestroy(&_vec);
899 : if (ierr)
900 : libmesh_warning("Warning: VecDestroy returned a non-zero error code which we ignored.");
901 : }
902 :
903 2997195 : this->_is_closed = this->_is_initialized = false;
904 :
905 81766 : _global_to_local_map.clear();
906 2997195 : }
907 :
908 :
909 :
910 : template <typename T>
911 : inline
912 5986415 : void PetscVector<T>::zero ()
913 : {
914 171112 : parallel_object_only();
915 :
916 171112 : libmesh_assert(this->closed());
917 :
918 5986415 : this->_restore_array();
919 :
920 171112 : PetscScalar z=0.;
921 :
922 5986415 : if (this->type() != GHOSTED)
923 5774400 : LibmeshPetscCall(VecSet (_vec, z));
924 : else
925 : {
926 : /* Vectors that include ghost values require a special
927 : handling. */
928 : Vec loc_vec;
929 212015 : LibmeshPetscCall(VecGhostGetLocalForm (_vec,&loc_vec));
930 212015 : LibmeshPetscCall(VecSet (loc_vec, z));
931 212015 : LibmeshPetscCall(VecGhostRestoreLocalForm (_vec, &loc_vec));
932 : }
933 5986415 : }
934 :
935 :
936 :
937 : template <typename T>
938 : inline
939 21030 : std::unique_ptr<NumericVector<T>> PetscVector<T>::zero_clone () const
940 : {
941 21030 : std::unique_ptr<NumericVector<T>> cloned_vector =
942 20398 : std::make_unique<PetscVector<T>>(this->comm(), this->type());
943 21030 : cloned_vector->init(*this);
944 21030 : return cloned_vector;
945 0 : }
946 :
947 :
948 :
949 : template <typename T>
950 : inline
951 337988 : std::unique_ptr<NumericVector<T>> PetscVector<T>::clone () const
952 : {
953 337988 : std::unique_ptr<NumericVector<T>> cloned_vector =
954 327962 : std::make_unique<PetscVector<T>>(this->comm(), this->type());
955 337988 : cloned_vector->init(*this, true);
956 337988 : *cloned_vector = *this;
957 337988 : return cloned_vector;
958 0 : }
959 :
960 :
961 :
962 : template <typename T>
963 : inline
964 183541778 : numeric_index_type PetscVector<T>::size () const
965 : {
966 182962819 : libmesh_assert (this->initialized());
967 :
968 183541778 : PetscInt petsc_size=0;
969 :
970 183541778 : if (!this->initialized())
971 0 : return 0;
972 :
973 183541778 : LibmeshPetscCall(VecGetSize(_vec, &petsc_size));
974 :
975 183541778 : return static_cast<numeric_index_type>(petsc_size);
976 : }
977 :
978 :
979 :
980 : template <typename T>
981 : inline
982 1021164 : numeric_index_type PetscVector<T>::local_size () const
983 : {
984 212356 : libmesh_assert (this->initialized());
985 :
986 1021164 : PetscInt petsc_size=0;
987 :
988 1021164 : LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_size));
989 :
990 1021164 : return static_cast<numeric_index_type>(petsc_size);
991 : }
992 :
993 :
994 :
995 : template <typename T>
996 : inline
997 32403114 : numeric_index_type PetscVector<T>::first_local_index () const
998 : {
999 4316297 : libmesh_assert (this->initialized());
1000 :
1001 4316297 : numeric_index_type first = 0;
1002 :
1003 32403114 : if (_array_is_present) // Can we use cached values?
1004 1568641 : first = _first;
1005 : else
1006 : {
1007 30834473 : PetscInt petsc_first=0, petsc_last=0;
1008 30834473 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1009 30834473 : first = static_cast<numeric_index_type>(petsc_first);
1010 : }
1011 :
1012 32403114 : return first;
1013 : }
1014 :
1015 :
1016 :
1017 : template <typename T>
1018 : inline
1019 6315632 : numeric_index_type PetscVector<T>::last_local_index () const
1020 : {
1021 3247717 : libmesh_assert (this->initialized());
1022 :
1023 3247717 : numeric_index_type last = 0;
1024 :
1025 6315632 : if (_array_is_present) // Can we use cached values?
1026 1450354 : last = _last;
1027 : else
1028 : {
1029 4865278 : PetscInt petsc_first=0, petsc_last=0;
1030 4865278 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1031 4865278 : last = static_cast<numeric_index_type>(petsc_last);
1032 : }
1033 :
1034 6315632 : return last;
1035 : }
1036 :
1037 :
1038 :
1039 : template <typename T>
1040 : inline
1041 6963289975 : numeric_index_type PetscVector<T>::map_global_to_local_index (const numeric_index_type i) const
1042 : {
1043 551817845 : libmesh_assert (this->initialized());
1044 :
1045 551817845 : numeric_index_type first=0;
1046 551817845 : numeric_index_type last=0;
1047 :
1048 6963289975 : if (_array_is_present) // Can we use cached values?
1049 : {
1050 6963289975 : first = _first;
1051 6963289975 : last = _last;
1052 : }
1053 : else
1054 : {
1055 0 : PetscInt petsc_first=0, petsc_last=0;
1056 0 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1057 0 : first = static_cast<numeric_index_type>(petsc_first);
1058 0 : last = static_cast<numeric_index_type>(petsc_last);
1059 : }
1060 :
1061 :
1062 6963289975 : if ((i>=first) && (i<last))
1063 : {
1064 6795540296 : return i-first;
1065 : }
1066 :
1067 5076926 : GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1068 : #ifndef NDEBUG
1069 5076926 : const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1070 5076926 : if (it == end)
1071 : {
1072 0 : std::ostringstream error_message;
1073 0 : error_message << "No index " << i << " in ghosted vector.\n"
1074 0 : << "Vector contains [" << first << ',' << last << ")\n";
1075 0 : GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1076 0 : if (b == end)
1077 : {
1078 0 : error_message << "And empty ghost array.\n";
1079 : }
1080 : else
1081 : {
1082 0 : error_message << "And ghost array {" << b->first;
1083 0 : for (++b; b != end; ++b)
1084 0 : error_message << ',' << b->first;
1085 0 : error_message << "}\n";
1086 : }
1087 :
1088 0 : libmesh_error_msg(error_message.str());
1089 : }
1090 5076926 : libmesh_assert (it != _global_to_local_map.end());
1091 : #endif
1092 167749679 : return it->second+last-first;
1093 : }
1094 :
1095 :
1096 :
1097 : template <typename T>
1098 : inline
1099 6506612800 : T PetscVector<T>::operator() (const numeric_index_type i) const
1100 : {
1101 6506612800 : this->_get_array(true);
1102 :
1103 6506612800 : const numeric_index_type local_index = this->map_global_to_local_index(i);
1104 :
1105 : #ifndef NDEBUG
1106 509667656 : if (this->type() == GHOSTED)
1107 : {
1108 505894854 : libmesh_assert_less (local_index, _local_size);
1109 : }
1110 : #endif
1111 :
1112 6506612800 : return static_cast<T>(_read_only_values[local_index]);
1113 : }
1114 :
1115 :
1116 :
1117 : template <typename T>
1118 : inline
1119 88409209 : void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
1120 : T * values) const
1121 : {
1122 88409209 : this->_get_array(true);
1123 :
1124 16571721 : const std::size_t num = index.size();
1125 :
1126 545086384 : for (std::size_t i=0; i<num; i++)
1127 : {
1128 498826900 : const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1129 : #ifndef NDEBUG
1130 42150189 : if (this->type() == GHOSTED)
1131 : {
1132 42140601 : libmesh_assert_less (local_index, _local_size);
1133 : }
1134 : #endif
1135 456677175 : values[i] = static_cast<T>(_read_only_values[local_index]);
1136 : }
1137 88409209 : }
1138 :
1139 :
1140 : template <typename T>
1141 : inline
1142 0 : PetscScalar * PetscVector<T>::get_array()
1143 : {
1144 0 : _get_array(false);
1145 0 : _values_manually_retrieved = true;
1146 :
1147 0 : return _values;
1148 : }
1149 :
1150 :
1151 : template <typename T>
1152 : inline
1153 0 : const PetscScalar * PetscVector<T>::get_array_read() const
1154 : {
1155 0 : _get_array(true);
1156 0 : _values_manually_retrieved = true;
1157 :
1158 0 : return _read_only_values;
1159 : }
1160 :
1161 : template <typename T>
1162 : inline
1163 0 : void PetscVector<T>::restore_array()
1164 : {
1165 : // \note \p _values_manually_retrieved needs to be set to \p false
1166 : // \e before calling \p _restore_array()!
1167 0 : _values_manually_retrieved = false;
1168 0 : _restore_array();
1169 0 : }
1170 :
1171 : template <typename T>
1172 : inline
1173 0 : Real PetscVector<T>::min () const
1174 : {
1175 0 : parallel_object_only();
1176 :
1177 0 : this->_restore_array();
1178 :
1179 0 : PetscInt index=0;
1180 0 : PetscReal returnval=0.;
1181 :
1182 0 : LibmeshPetscCall(VecMin (_vec, &index, &returnval));
1183 :
1184 : // this return value is correct: VecMin returns a PetscReal
1185 0 : return static_cast<Real>(returnval);
1186 : }
1187 :
1188 :
1189 :
1190 : template <typename T>
1191 : inline
1192 0 : Real PetscVector<T>::max() const
1193 : {
1194 0 : parallel_object_only();
1195 :
1196 0 : this->_restore_array();
1197 :
1198 0 : PetscInt index=0;
1199 0 : PetscReal returnval=0.;
1200 :
1201 0 : LibmeshPetscCall(VecMax (_vec, &index, &returnval));
1202 :
1203 : // this return value is correct: VecMax returns a PetscReal
1204 0 : return static_cast<Real>(returnval);
1205 : }
1206 :
1207 :
1208 :
1209 : template <typename T>
1210 : inline
1211 1328136 : void PetscVector<T>::swap (NumericVector<T> & other)
1212 : {
1213 28720 : parallel_object_only();
1214 :
1215 28720 : NumericVector<T>::swap(other);
1216 :
1217 28720 : PetscVector<T> & v = cast_ref<PetscVector<T> &>(other);
1218 :
1219 28720 : std::swap(_vec, v._vec);
1220 28720 : std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1221 28720 : std::swap(_global_to_local_map, v._global_to_local_map);
1222 :
1223 : #ifdef LIBMESH_HAVE_CXX11_THREAD
1224 : // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1225 1355224 : _array_is_present = v._array_is_present.exchange(_array_is_present);
1226 : #else
1227 : std::swap(_array_is_present, v._array_is_present);
1228 : #endif
1229 :
1230 28720 : std::swap(_local_form, v._local_form);
1231 28720 : std::swap(_values, v._values);
1232 1328136 : }
1233 :
1234 :
1235 :
1236 : template <typename T>
1237 : inline
1238 27455 : std::size_t PetscVector<T>::max_allowed_id () const
1239 : {
1240 : // The PetscInt type is used for indexing, it may be either a signed
1241 : // 4-byte or 8-byte integer depending on how PETSc is configured.
1242 27455 : return std::numeric_limits<PetscInt>::max();
1243 : }
1244 :
1245 :
1246 :
1247 : #ifdef LIBMESH_HAVE_CXX11
1248 : static_assert(sizeof(PetscInt) == sizeof(numeric_index_type),
1249 : "PETSc and libMesh integer sizes must match!");
1250 : #endif
1251 :
1252 :
1253 : inline
1254 6986722 : PetscInt * numeric_petsc_cast(const numeric_index_type * p)
1255 : {
1256 6986722 : return reinterpret_cast<PetscInt *>(const_cast<numeric_index_type *>(p));
1257 : }
1258 :
1259 : } // namespace libMesh
1260 :
1261 :
1262 : #endif // #ifdef LIBMESH_HAVE_PETSC
1263 : #endif // LIBMESH_PETSC_VECTOR_H
|
