Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
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 : * Supported assignments based on ParallelType combination (note that we lump ghosted into
125 : * parallel for this method documentation):
126 : * - Assign from parallel to parallel
127 : * - Assign from serial to serial
128 : * - Assign from parallel to serial
129 : * \returns A reference to *this as the derived type.
130 : */
131 : PetscVector<T> & operator= (const PetscVector<T> & v);
132 :
133 : /**
134 : * This class manages a C-style struct (Vec) manually, so we
135 : * don't want to allow any automatic copy/move functions to be
136 : * generated, and we can't default the destructor.
137 : */
138 : PetscVector (PetscVector &&) = delete;
139 : PetscVector (const PetscVector &) = delete;
140 : PetscVector & operator= (PetscVector &&) = delete;
141 : virtual ~PetscVector ();
142 :
143 : virtual void close () override;
144 :
145 : /**
146 : * clear() is called from the destructor, so it should not throw.
147 : */
148 : virtual void clear () noexcept override;
149 :
150 : virtual void zero () override;
151 :
152 : virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
153 :
154 : virtual std::unique_ptr<NumericVector<T>> clone () const override;
155 :
156 : virtual void init (const numeric_index_type N,
157 : const numeric_index_type n_local,
158 : const bool fast=false,
159 : const ParallelType type=AUTOMATIC) override;
160 :
161 : virtual void init (const numeric_index_type N,
162 : const bool fast=false,
163 : const ParallelType type=AUTOMATIC) override;
164 :
165 : virtual void init (const numeric_index_type N,
166 : const numeric_index_type n_local,
167 : const std::vector<numeric_index_type> & ghost,
168 : const bool fast = false,
169 : const ParallelType = AUTOMATIC) override;
170 :
171 : virtual void init (const NumericVector<T> & other,
172 : const bool fast = false) override;
173 :
174 : virtual NumericVector<T> & operator= (const T s) override;
175 :
176 : virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
177 :
178 : virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
179 :
180 : virtual Real min () const override;
181 :
182 : virtual Real max () const override;
183 :
184 : virtual T sum () const override;
185 :
186 : virtual Real l1_norm () const override;
187 :
188 : virtual Real l2_norm () const override;
189 :
190 : virtual Real linfty_norm () const override;
191 :
192 : virtual numeric_index_type size () const override;
193 :
194 : virtual numeric_index_type local_size() const override;
195 :
196 : virtual numeric_index_type first_local_index() const override;
197 :
198 : virtual numeric_index_type last_local_index() const override;
199 :
200 : /**
201 : * \returns The local index corresponding to global index \p i.
202 : *
203 : * If the index is not a ghost cell, this is done by subtraction the
204 : * number of the first local index. If it is a ghost cell, it has
205 : * to be looked up in the map.
206 : */
207 : numeric_index_type map_global_to_local_index(const numeric_index_type i) const;
208 :
209 : virtual T operator() (const numeric_index_type i) const override;
210 :
211 : virtual void get(const std::vector<numeric_index_type> & index,
212 : T * values) const override;
213 :
214 : /**
215 : * Get read/write access to the raw PETSc Vector data array.
216 : *
217 : * \note This is an advanced interface. In general you should avoid
218 : * using it unless you know what you are doing!
219 : */
220 : PetscScalar * get_array();
221 :
222 : /**
223 : * Get read only access to the raw PETSc Vector data array.
224 : *
225 : * \note This is an advanced interface. In general you should avoid
226 : * using it unless you know what you are doing!
227 : */
228 : const PetscScalar * get_array_read() const;
229 :
230 : /**
231 : * Restore the data array.
232 : *
233 : * \note This MUST be called after get_array() or get_array_read()
234 : * and before using any other interface functions on PetscVector.
235 : */
236 : void restore_array();
237 :
238 : virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
239 :
240 : virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
241 :
242 : virtual void reciprocal() override;
243 :
244 : virtual void conjugate() override;
245 :
246 : virtual void set (const numeric_index_type i,
247 : const T value) override;
248 :
249 : virtual void add (const numeric_index_type i,
250 : const T value) override;
251 :
252 : virtual void add (const T s) override;
253 :
254 : virtual void add (const NumericVector<T> & v) override;
255 :
256 : virtual void add (const T a, const NumericVector<T> & v) override;
257 :
258 : /**
259 : * We override two NumericVector<T>::add_vector() methods but don't
260 : * want to hide the other defaults.
261 : */
262 : using NumericVector<T>::add_vector;
263 :
264 : virtual void add_vector (const T * v,
265 : const std::vector<numeric_index_type> & dof_indices) override;
266 :
267 : virtual void add_vector (const NumericVector<T> & v,
268 : const SparseMatrix<T> & A) override;
269 :
270 : virtual void add_vector_transpose (const NumericVector<T> & v,
271 : const SparseMatrix<T> & A) override;
272 :
273 : /**
274 : * \f$ U \leftarrow U + A^H v \f$.
275 : *
276 : * Adds the product of the conjugate-transpose of \p SparseMatrix \p
277 : * A and a \p NumericVector \p v to \p this.
278 : */
279 : void add_vector_conjugate_transpose (const NumericVector<T> & v,
280 : const SparseMatrix<T> & A);
281 :
282 : /**
283 : * We override one NumericVector<T>::insert() method but don't want
284 : * to hide the other defaults
285 : */
286 : using NumericVector<T>::insert;
287 :
288 : virtual void insert (const T * v,
289 : const std::vector<numeric_index_type> & dof_indices) override;
290 :
291 : virtual void scale (const T factor) override;
292 :
293 : virtual NumericVector<T> & operator *= (const NumericVector<T> & v) override;
294 :
295 : virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
296 :
297 : virtual void abs() override;
298 :
299 : virtual T dot(const NumericVector<T> & v) const override;
300 :
301 : /**
302 : * \returns The dot product of (*this) with the vector \p v.
303 : *
304 : * \note Does *not* use the complex-conjugate of v in the complex-valued case.
305 : */
306 : T indefinite_dot(const NumericVector<T> & v) const;
307 :
308 : virtual void localize (std::vector<T> & v_local) const override;
309 :
310 : virtual void localize (NumericVector<T> & v_local) const override;
311 :
312 : virtual void localize (NumericVector<T> & v_local,
313 : const std::vector<numeric_index_type> & send_list) const override;
314 :
315 : virtual void localize (std::vector<T> & v_local,
316 : const std::vector<numeric_index_type> & indices) const override;
317 :
318 : virtual void localize (const numeric_index_type first_local_idx,
319 : const numeric_index_type last_local_idx,
320 : const std::vector<numeric_index_type> & send_list) override;
321 :
322 : virtual void localize_to_one (std::vector<T> & v_local,
323 : const processor_id_type proc_id=0) const override;
324 :
325 : virtual void pointwise_mult (const NumericVector<T> & vec1,
326 : const NumericVector<T> & vec2) override;
327 :
328 : virtual void pointwise_divide (const NumericVector<T> & vec1,
329 : const NumericVector<T> & vec2) override;
330 :
331 : virtual void print_matlab(const std::string & name = "") const override;
332 :
333 : virtual void create_subvector(NumericVector<T> & subvector,
334 : const std::vector<numeric_index_type> & rows,
335 : bool supplying_global_rows = true) const override;
336 :
337 : virtual void swap (NumericVector<T> & v) override;
338 :
339 : virtual std::size_t max_allowed_id() const override;
340 :
341 : /**
342 : * \returns The raw PETSc Vec pointer.
343 : *
344 : * \note This is generally not required in user-level code.
345 : *
346 : * \note Don't do anything crazy like calling VecDestroy() on
347 : * it, or very bad things will likely happen!
348 : */
349 571551 : Vec vec () { libmesh_assert (_vec); return _vec; }
350 :
351 4369175 : Vec vec () const { libmesh_assert (_vec); return _vec; }
352 :
353 : virtual std::unique_ptr<NumericVector<T>>
354 : get_subvector(const std::vector<numeric_index_type> & rows) override;
355 :
356 : virtual void restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
357 : const std::vector<numeric_index_type> & rows) override;
358 :
359 : private:
360 :
361 : /**
362 : * Actual PETSc vector datatype to hold vector entries.
363 : */
364 : Vec _vec;
365 :
366 : /**
367 : * If \p true, the actual PETSc array of the values of the vector is
368 : * currently accessible. That means that the members \p _local_form
369 : * and \p _values are valid.
370 : */
371 : #ifdef LIBMESH_HAVE_CXX11_THREAD
372 : // We can't use std::atomic_flag here because we need load and store operations.
373 : mutable std::atomic<bool> _array_is_present;
374 : #else
375 : mutable bool _array_is_present;
376 : #endif
377 :
378 : /**
379 : * First local index.
380 : *
381 : * Only valid when _array_is_present
382 : */
383 : mutable numeric_index_type _first;
384 :
385 : /**
386 : * Last local index.
387 : *
388 : * Only valid when _array_is_present
389 : */
390 : mutable numeric_index_type _last;
391 :
392 : /**
393 : * Size of the local values from _get_array()
394 : */
395 : mutable numeric_index_type _local_size;
396 :
397 : /**
398 : * PETSc vector datatype to hold the local form of a ghosted vector.
399 : * The contents of this field are only valid if the vector is
400 : * ghosted and \p _array_is_present is \p true.
401 : */
402 : mutable Vec _local_form;
403 :
404 : /**
405 : * Pointer to the actual PETSc array of the values of the vector.
406 : * This pointer is only valid if \p _array_is_present is \p true.
407 : * We're using PETSc's VecGetArrayRead() function, which requires a
408 : * constant PetscScalar *, but _get_array and _restore_array are
409 : * const member functions, so _values also needs to be mutable
410 : * (otherwise it is a "const PetscScalar * const" in that context).
411 : */
412 : mutable const PetscScalar * _read_only_values;
413 :
414 : /**
415 : * Pointer to the actual PETSc array of the values of the vector.
416 : * This pointer is only valid if \p _array_is_present is \p true.
417 : * We're using PETSc's VecGetArrayRead() function, which requires a
418 : * constant PetscScalar *, but _get_array and _restore_array are
419 : * const member functions, so _values also needs to be mutable
420 : * (otherwise it is a "const PetscScalar * const" in that context).
421 : */
422 : mutable PetscScalar * _values;
423 :
424 : /**
425 : * Mutex for _get_array and _restore_array. This is part of the
426 : * object to keep down thread contention when reading from multiple
427 : * PetscVectors simultaneously
428 : */
429 : mutable std::mutex _petsc_get_restore_array_mutex;
430 :
431 : /**
432 : * Queries the array (and the local form if the vector is ghosted)
433 : * from PETSc.
434 : *
435 : * \param read_only Whether or not a read only copy of the raw data
436 : */
437 : void _get_array(bool read_only) const;
438 :
439 : /**
440 : * Restores the array (and the local form if the vector is ghosted)
441 : * to PETSc.
442 : */
443 : void _restore_array() const;
444 :
445 : /**
446 : * Type for map that maps global to local ghost cells.
447 : */
448 : typedef std::unordered_map<numeric_index_type,numeric_index_type> GlobalToLocalMap;
449 :
450 : /**
451 : * Map that maps global to local ghost cells (will be empty if not
452 : * in ghost cell mode).
453 : */
454 : GlobalToLocalMap _global_to_local_map;
455 :
456 : /**
457 : * This boolean value should only be set to false
458 : * for the constructor which takes a PETSc Vec object.
459 : */
460 : bool _destroy_vec_on_exit;
461 :
462 : /**
463 : * Whether or not the data array has been manually retrieved using
464 : * get_array() or get_array_read()
465 : */
466 : mutable bool _values_manually_retrieved;
467 :
468 : /**
469 : * Whether or not the data array is for read only access
470 : */
471 : mutable bool _values_read_only;
472 :
473 : /**
474 : * \returns A norm of the vector, where the type of norm to compute is
475 : * determined by the template parameter N of the PETSc-defined type NormType.
476 : * The valid template arguments are NORM_1, NORM_2 and NORM_INFINITY, as used
477 : * to define l1_norm(), l2_norm() and linfty_norm().
478 : */
479 : template <NormType N> Real norm () const;
480 : };
481 :
482 :
483 : /*----------------------- Inline functions ----------------------------------*/
484 :
485 :
486 :
487 : template <typename T>
488 : inline
489 2296304 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in, const ParallelType ptype) :
490 : NumericVector<T>(comm_in, ptype),
491 2173572 : _vec(nullptr),
492 : _array_is_present(false),
493 2173572 : _first(0),
494 2173572 : _last(0),
495 2173572 : _local_size(0),
496 2173572 : _local_form(nullptr),
497 2173572 : _read_only_values(nullptr),
498 2173572 : _values(nullptr),
499 2173572 : _global_to_local_map(),
500 2173572 : _destroy_vec_on_exit(true),
501 2173572 : _values_manually_retrieved(false),
502 2296304 : _values_read_only(false)
503 : {
504 2228160 : this->_type = ptype;
505 2296304 : }
506 :
507 :
508 :
509 : template <typename T>
510 : inline
511 0 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
512 : const numeric_index_type n,
513 : const ParallelType ptype) :
514 : NumericVector<T>(comm_in, ptype),
515 0 : _vec(nullptr),
516 : _array_is_present(false),
517 0 : _first(0),
518 0 : _last(0),
519 0 : _local_size(0),
520 0 : _local_form(nullptr),
521 0 : _read_only_values(nullptr),
522 0 : _values(nullptr),
523 0 : _global_to_local_map(),
524 0 : _destroy_vec_on_exit(true),
525 0 : _values_manually_retrieved(false),
526 0 : _values_read_only(false)
527 : {
528 0 : this->init(n, n, false, ptype);
529 0 : }
530 :
531 :
532 :
533 : template <typename T>
534 : inline
535 0 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
536 : const numeric_index_type n,
537 : const numeric_index_type n_local,
538 : const ParallelType ptype) :
539 : NumericVector<T>(comm_in, ptype),
540 0 : _vec(nullptr),
541 : _array_is_present(false),
542 0 : _first(0),
543 0 : _last(0),
544 0 : _local_size(0),
545 0 : _local_form(nullptr),
546 0 : _read_only_values(nullptr),
547 0 : _values(nullptr),
548 0 : _global_to_local_map(),
549 0 : _destroy_vec_on_exit(true),
550 0 : _values_manually_retrieved(false),
551 0 : _values_read_only(false)
552 : {
553 0 : this->init(n, n_local, false, ptype);
554 0 : }
555 :
556 :
557 :
558 : template <typename T>
559 : inline
560 0 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
561 : const numeric_index_type n,
562 : const numeric_index_type n_local,
563 : const std::vector<numeric_index_type> & ghost,
564 : const ParallelType ptype) :
565 : NumericVector<T>(comm_in, ptype),
566 0 : _vec(nullptr),
567 : _array_is_present(false),
568 0 : _first(0),
569 0 : _last(0),
570 0 : _local_size(0),
571 0 : _local_form(nullptr),
572 0 : _read_only_values(nullptr),
573 0 : _values(nullptr),
574 0 : _global_to_local_map(),
575 0 : _destroy_vec_on_exit(true),
576 0 : _values_manually_retrieved(false),
577 0 : _values_read_only(false)
578 : {
579 0 : this->init(n, n_local, ghost, false, ptype);
580 0 : }
581 :
582 :
583 :
584 :
585 :
586 : template <typename T>
587 : inline
588 1010618 : PetscVector<T>::PetscVector (Vec v,
589 : const Parallel::Communicator & comm_in) :
590 : NumericVector<T>(comm_in, AUTOMATIC),
591 968450 : _vec(v),
592 : _array_is_present(false),
593 968450 : _first(0),
594 968450 : _last(0),
595 968450 : _local_size(0),
596 968450 : _local_form(nullptr),
597 968450 : _read_only_values(nullptr),
598 968450 : _values(nullptr),
599 968450 : _global_to_local_map(),
600 968450 : _destroy_vec_on_exit(false),
601 968450 : _values_manually_retrieved(false),
602 1010618 : _values_read_only(false)
603 : {
604 1010618 : this->_is_closed = true;
605 1010618 : this->_is_initialized = true;
606 :
607 : /* We need to ask PETSc about the (local to global) ghost value
608 : mapping and create the inverse mapping out of it. */
609 1010618 : PetscInt petsc_local_size=0;
610 1010618 : LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_local_size));
611 :
612 : // Get the vector type from PETSc.
613 : VecType ptype;
614 1010618 : LibmeshPetscCall(VecGetType(_vec, &ptype));
615 :
616 : // PETSc supports some exotic types nowadays. Just use Vec sizes to
617 : // determine whether we're SERIAL vs PARALLEL-or-GHOSTED.
618 1010618 : PetscInt petsc_global_size = 0;
619 1010618 : LibmeshPetscCall(VecGetSize(_vec, &petsc_global_size));
620 :
621 : // If we have a parallel Vec with all its DoFs on one processor, we
622 : // might be unable to tell on that processor that it's not a serial
623 : // vector unless we communicate.
624 1010618 : bool is_serial = (petsc_local_size == petsc_global_size);
625 1010618 : comm_in.min(is_serial);
626 :
627 : #ifdef DEBUG
628 21696 : dof_id_type sum_of_local_sizes = petsc_local_size;
629 21696 : comm_in.sum(sum_of_local_sizes);
630 21696 : libmesh_assert_equal_to(sum_of_local_sizes,
631 : cast_int<dof_id_type>(petsc_global_size));
632 : #endif
633 :
634 : #if PETSC_RELEASE_GREATER_EQUALS(3, 21, 0)
635 : // Fande only implemented VecGhostGetGhostIS for VECMPI
636 : if (std::strcmp(ptype, VECMPI) == 0)
637 : {
638 : IS ghostis;
639 : LibmeshPetscCall(VecGhostGetGhostIS(_vec, &ghostis));
640 :
641 : Vec localrep;
642 : LibmeshPetscCall(VecGhostGetLocalForm(_vec, &localrep));
643 :
644 : // If is a sparsely stored vector, set up our new mapping
645 : // Only checking mapping is not enough to determine if a vec is ghosted
646 : // We need to check if vec has a local representation
647 : if (ghostis && localrep)
648 : {
649 : PetscInt ghost_size;
650 : LibmeshPetscCall(ISGetSize(ghostis, &ghost_size));
651 :
652 : const PetscInt * indices;
653 : LibmeshPetscCall(ISGetIndices(ghostis, &indices));
654 :
655 : for (const auto i : make_range(ghost_size))
656 : _global_to_local_map[indices[i]] = i;
657 : this->_type = GHOSTED;
658 : LibmeshPetscCall(ISRestoreIndices(ghostis, &indices));
659 : }
660 : else
661 : this->_type = PARALLEL;
662 : }
663 : else if (!is_serial)
664 : #else // PETSc < 3.21.0
665 1010618 : if (!is_serial)
666 : #endif
667 : {
668 995772 : ISLocalToGlobalMapping mapping = nullptr;
669 995772 : Vec localrep = nullptr;
670 :
671 : // Non-VECMPI types won't even let us call VecGetLocalToGlobalMapping;
672 : // we'll just assume for now that no such type is ghosted
673 995772 : if (std::strcmp(ptype, VECMPI) == 0)
674 : {
675 995772 : LibmeshPetscCall(VecGetLocalToGlobalMapping(_vec, &mapping));
676 995772 : LibmeshPetscCall(VecGhostGetLocalForm(_vec,&localrep));
677 : }
678 :
679 : // If is a sparsely stored vector, set up our new mapping
680 : // Only checking mapping is not enough to determine if a vec is ghosted
681 : // We need to check if vec has a local representation
682 995772 : if (mapping && localrep)
683 : {
684 0 : const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
685 0 : const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
686 : PetscInt n;
687 0 : LibmeshPetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
688 :
689 0 : const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
690 : const PetscInt * indices;
691 0 : LibmeshPetscCall(ISLocalToGlobalMappingGetIndices(mapping,&indices));
692 :
693 0 : for (numeric_index_type i=ghost_begin; i<ghost_end; i++)
694 0 : _global_to_local_map[indices[i]] = i-my_local_size;
695 0 : this->_type = GHOSTED;
696 0 : LibmeshPetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
697 0 : }
698 : else
699 995772 : this->_type = PARALLEL;
700 :
701 995772 : LibmeshPetscCall(VecGhostRestoreLocalForm(_vec,&localrep));
702 : }
703 : else
704 14846 : this->_type = SERIAL;
705 1010618 : }
706 :
707 :
708 :
709 :
710 : template <typename T>
711 : inline
712 4653532 : PetscVector<T>::~PetscVector ()
713 : {
714 2835301 : this->clear ();
715 4728468 : }
716 :
717 :
718 :
719 : template <typename T>
720 : inline
721 945867 : void PetscVector<T>::init (const numeric_index_type n,
722 : const numeric_index_type n_local,
723 : const bool fast,
724 : const ParallelType ptype)
725 : {
726 27860 : parallel_object_only();
727 :
728 945867 : PetscInt petsc_n=static_cast<PetscInt>(n);
729 :
730 : // Clear initialized vectors
731 945867 : if (this->initialized())
732 31058 : this->clear();
733 :
734 945867 : if (ptype == AUTOMATIC)
735 : {
736 4830 : if (n == n_local)
737 621 : this->_type = SERIAL;
738 : else
739 4209 : this->_type = PARALLEL;
740 : }
741 : else
742 941037 : this->_type = ptype;
743 :
744 : // We should have been given consistent settings
745 27860 : libmesh_assert ((this->_type==SERIAL && n==n_local) ||
746 : (this->n_processors()==1 && n==n_local) ||
747 : this->_type==PARALLEL);
748 :
749 : // create a sequential vector if on only 1 processor, or if asked
750 945867 : if (this->_type == SERIAL || (this->n_processors() == 1))
751 : {
752 23258 : LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
753 23258 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
754 23258 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetFromOptions (_vec));
755 : }
756 : // or create an MPI-enabled PARALLEL vector w/o ghosting if asked
757 922609 : else if (this->_type == PARALLEL)
758 : {
759 : #ifdef LIBMESH_HAVE_MPI
760 27810 : PetscInt petsc_n_local=cast_int<PetscInt>(n_local);
761 27810 : libmesh_assert_less_equal (n_local, n);
762 : // Use more generic function instead of VecCreateSeq/MPI
763 922609 : LibmeshPetscCall(VecCreate(this->comm().get(), &_vec));
764 922609 : LibmeshPetscCall(VecSetSizes(_vec, petsc_n_local, petsc_n));
765 : #else
766 : libmesh_assert_equal_to (n_local, n);
767 : LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
768 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
769 : #endif
770 922609 : LibmeshPetscCall(VecSetFromOptions(_vec));
771 : }
772 : // or yell because we don't know what to do
773 : else
774 0 : libmesh_error_msg("Unsupported type " <<
775 : Utility::enum_to_string(this->_type) <<
776 : " for parallel init with no ghost indices supplied");
777 :
778 945867 : this->_is_initialized = true;
779 945867 : this->_is_closed = true;
780 :
781 :
782 945867 : if (fast == false)
783 924958 : this->zero ();
784 945867 : }
785 :
786 :
787 :
788 : template <typename T>
789 : inline
790 1606 : void PetscVector<T>::init (const numeric_index_type n,
791 : const bool fast,
792 : const ParallelType ptype)
793 : {
794 2150 : this->init(n,n,fast,ptype);
795 2134 : }
796 :
797 :
798 :
799 : template <typename T>
800 : inline
801 555106 : void PetscVector<T>::init (const numeric_index_type n,
802 : const numeric_index_type n_local,
803 : const std::vector<numeric_index_type> & ghost,
804 : const bool fast,
805 : const ParallelType ptype)
806 : {
807 17542 : parallel_object_only();
808 :
809 555106 : if (this->comm().size() == 1)
810 : {
811 0 : libmesh_assert(ghost.empty());
812 7906 : this->init(n, n_local, fast, ptype);
813 7906 : return;
814 : }
815 :
816 547200 : PetscInt petsc_n=static_cast<PetscInt>(n);
817 547200 : PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
818 34896 : PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
819 :
820 : // If the mesh is not disjoint, every processor will either have
821 : // all the dofs, none of the dofs, or some non-zero dofs at the
822 : // boundary between processors.
823 : //
824 : // However we can't assert this, because someone might want to
825 : // construct a GHOSTED vector which doesn't include neighbor element
826 : // dofs. Boyce tried to do so in user code, and we're going to want
827 : // to do so in System::project_vector().
828 : //
829 : // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
830 :
831 : libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
832 :
833 547200 : PetscInt * petsc_ghost = ghost.empty() ? LIBMESH_PETSC_NULLPTR :
834 16830 : const_cast<PetscInt *>(reinterpret_cast<const PetscInt *>(ghost.data()));
835 :
836 : // Clear initialized vectors
837 547200 : if (this->initialized())
838 27690 : this->clear();
839 :
840 17542 : libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
841 547200 : this->_type = GHOSTED;
842 :
843 : /* Make the global-to-local ghost cell map. */
844 42549329 : for (auto i : index_range(ghost))
845 : {
846 42002129 : _global_to_local_map[ghost[i]] = i;
847 : }
848 :
849 : /* Create vector. */
850 547200 : LibmeshPetscCall(VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
851 : petsc_n_ghost, petsc_ghost,
852 : &_vec));
853 :
854 : // Add a prefix so that we can at least distinguish a ghosted vector from a
855 : // nonghosted vector when using a petsc option.
856 : // PETSc does not fully support VecGhost on GPU yet. This change allows us to
857 : // trigger a nonghosted vector to use GPU without bothering the ghosted vectors.
858 547200 : LibmeshPetscCall(PetscObjectAppendOptionsPrefix((PetscObject)_vec,"ghost_"));
859 :
860 547200 : LibmeshPetscCall(VecSetFromOptions (_vec));
861 :
862 547200 : this->_is_initialized = true;
863 547200 : this->_is_closed = true;
864 :
865 547200 : if (fast == false)
866 191040 : this->zero ();
867 : }
868 :
869 :
870 :
871 : template <typename T>
872 : inline
873 442504 : void PetscVector<T>::init (const NumericVector<T> & other,
874 : const bool fast)
875 : {
876 13032 : parallel_object_only();
877 :
878 : // Clear initialized vectors
879 442504 : if (this->initialized())
880 32690 : this->clear();
881 :
882 13032 : const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
883 :
884 : // Other vector should restore array.
885 442504 : if (v.initialized())
886 : {
887 442504 : v._restore_array();
888 : }
889 :
890 13032 : this->_global_to_local_map = v._global_to_local_map;
891 :
892 : // Even if we're initializing sizes based on an uninitialized or
893 : // unclosed vector, *this* vector is being initialized now and is
894 : // initially closed.
895 442504 : this->_is_closed = true; // v._is_closed;
896 442504 : this->_is_initialized = true; // v._is_initialized;
897 :
898 442504 : this->_type = v._type;
899 :
900 : // We want to have a valid Vec, even if it's initially of size zero
901 442504 : LibmeshPetscCall(VecDuplicate (v._vec, &this->_vec));
902 :
903 442504 : if (fast == false)
904 53766 : this->zero ();
905 442504 : }
906 :
907 :
908 :
909 : template <typename T>
910 : inline
911 3544041 : void PetscVector<T>::close ()
912 : {
913 96392 : parallel_object_only();
914 :
915 3544041 : this->_restore_array();
916 :
917 3544041 : VecAssemblyBeginEnd(this->comm(), _vec);
918 :
919 96392 : if (this->is_effectively_ghosted())
920 801770 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
921 :
922 3544041 : this->_is_closed = true;
923 3544041 : }
924 :
925 :
926 :
927 : template <typename T>
928 : inline
929 3106811 : void PetscVector<T>::clear () noexcept
930 : {
931 84706 : exceptionless_parallel_object_only();
932 :
933 3106811 : if (this->initialized())
934 2946049 : this->_restore_array();
935 :
936 3106811 : if ((this->initialized()) && (this->_destroy_vec_on_exit))
937 : {
938 : // If we encounter an error here, print a warning but otherwise
939 : // keep going since we may be recovering from an exception.
940 1935571 : PetscErrorCode ierr = VecDestroy(&_vec);
941 : if (ierr)
942 : libmesh_warning("Warning: VecDestroy returned a non-zero error code which we ignored.");
943 : }
944 :
945 3106811 : this->_is_closed = this->_is_initialized = false;
946 :
947 84706 : _global_to_local_map.clear();
948 3106811 : }
949 :
950 :
951 :
952 : template <typename T>
953 : inline
954 6053370 : void PetscVector<T>::zero ()
955 : {
956 172808 : parallel_object_only();
957 :
958 172808 : libmesh_assert(this->closed());
959 :
960 6053370 : this->_restore_array();
961 :
962 172808 : PetscScalar z=0.;
963 :
964 172808 : if (!this->is_effectively_ghosted())
965 5808924 : LibmeshPetscCall(VecSet (_vec, z));
966 : else
967 : {
968 : /* Vectors that include ghost values require a special
969 : handling. */
970 : Vec loc_vec;
971 244446 : LibmeshPetscCall(VecGhostGetLocalForm (_vec,&loc_vec));
972 244446 : LibmeshPetscCall(VecSet (loc_vec, z));
973 244446 : LibmeshPetscCall(VecGhostRestoreLocalForm (_vec, &loc_vec));
974 : }
975 6053370 : }
976 :
977 :
978 :
979 : template <typename T>
980 : inline
981 21076 : std::unique_ptr<NumericVector<T>> PetscVector<T>::zero_clone () const
982 : {
983 21076 : std::unique_ptr<NumericVector<T>> cloned_vector =
984 20454 : std::make_unique<PetscVector<T>>(this->comm(), this->type());
985 21076 : cloned_vector->init(*this);
986 21076 : return cloned_vector;
987 0 : }
988 :
989 :
990 :
991 : template <typename T>
992 : inline
993 355628 : std::unique_ptr<NumericVector<T>> PetscVector<T>::clone () const
994 : {
995 355628 : std::unique_ptr<NumericVector<T>> cloned_vector =
996 345098 : std::make_unique<PetscVector<T>>(this->comm(), this->type());
997 355628 : cloned_vector->init(*this, true);
998 355628 : *cloned_vector = *this;
999 355628 : return cloned_vector;
1000 0 : }
1001 :
1002 :
1003 :
1004 : template <typename T>
1005 : inline
1006 190760495 : numeric_index_type PetscVector<T>::size () const
1007 : {
1008 190152142 : libmesh_assert (this->initialized());
1009 :
1010 190760495 : PetscInt petsc_size=0;
1011 :
1012 190760495 : if (!this->initialized())
1013 0 : return 0;
1014 :
1015 190760495 : LibmeshPetscCall(VecGetSize(_vec, &petsc_size));
1016 :
1017 190760495 : return static_cast<numeric_index_type>(petsc_size);
1018 : }
1019 :
1020 :
1021 :
1022 : template <typename T>
1023 : inline
1024 8356631 : numeric_index_type PetscVector<T>::local_size () const
1025 : {
1026 213190 : libmesh_assert (this->initialized());
1027 :
1028 8356631 : PetscInt petsc_size=0;
1029 :
1030 8356631 : LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_size));
1031 :
1032 8356631 : return static_cast<numeric_index_type>(petsc_size);
1033 : }
1034 :
1035 :
1036 :
1037 : template <typename T>
1038 : inline
1039 7119691 : numeric_index_type PetscVector<T>::first_local_index () const
1040 : {
1041 3321756 : libmesh_assert (this->initialized());
1042 :
1043 3321756 : numeric_index_type first = 0;
1044 :
1045 7119691 : if (_array_is_present) // Can we use cached values?
1046 1572489 : first = _first;
1047 : else
1048 : {
1049 5547202 : PetscInt petsc_first=0, petsc_last=0;
1050 5547202 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1051 5547202 : first = static_cast<numeric_index_type>(petsc_first);
1052 : }
1053 :
1054 7119691 : return first;
1055 : }
1056 :
1057 :
1058 :
1059 : template <typename T>
1060 : inline
1061 6979349 : numeric_index_type PetscVector<T>::last_local_index () const
1062 : {
1063 3311104 : libmesh_assert (this->initialized());
1064 :
1065 3311104 : numeric_index_type last = 0;
1066 :
1067 6979349 : if (_array_is_present) // Can we use cached values?
1068 1454202 : last = _last;
1069 : else
1070 : {
1071 5525147 : PetscInt petsc_first=0, petsc_last=0;
1072 5525147 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1073 5525147 : last = static_cast<numeric_index_type>(petsc_last);
1074 : }
1075 :
1076 6979349 : return last;
1077 : }
1078 :
1079 :
1080 :
1081 : template <typename T>
1082 : inline
1083 7264366228 : numeric_index_type PetscVector<T>::map_global_to_local_index (const numeric_index_type i) const
1084 : {
1085 573656336 : libmesh_assert (this->initialized());
1086 :
1087 573656336 : numeric_index_type first=0;
1088 573656336 : numeric_index_type last=0;
1089 :
1090 7264366228 : if (_array_is_present) // Can we use cached values?
1091 : {
1092 7264366228 : first = _first;
1093 7264366228 : last = _last;
1094 : }
1095 : else
1096 : {
1097 0 : PetscInt petsc_first=0, petsc_last=0;
1098 0 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1099 0 : first = static_cast<numeric_index_type>(petsc_first);
1100 0 : last = static_cast<numeric_index_type>(petsc_last);
1101 : }
1102 :
1103 :
1104 7264366228 : if ((i>=first) && (i<last))
1105 : {
1106 7072572232 : return i-first;
1107 : }
1108 :
1109 5915587 : GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1110 : #ifndef NDEBUG
1111 5915587 : const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1112 5915587 : if (it == end)
1113 : {
1114 0 : std::ostringstream error_message;
1115 0 : error_message << "No index " << i << " in ghosted vector.\n"
1116 0 : << "Vector contains [" << first << ',' << last << ")\n";
1117 0 : GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1118 0 : if (b == end)
1119 : {
1120 0 : error_message << "And empty ghost array.\n";
1121 : }
1122 : else
1123 : {
1124 0 : error_message << "And ghost array {" << b->first;
1125 0 : for (++b; b != end; ++b)
1126 0 : error_message << ',' << b->first;
1127 0 : error_message << "}\n";
1128 : }
1129 :
1130 0 : libmesh_error_msg(error_message.str());
1131 : }
1132 5915587 : libmesh_assert (it != _global_to_local_map.end());
1133 : #endif
1134 191793996 : return it->second+last-first;
1135 : }
1136 :
1137 :
1138 :
1139 : template <typename T>
1140 : inline
1141 6783813454 : T PetscVector<T>::operator() (const numeric_index_type i) const
1142 : {
1143 6783813454 : this->_get_array(true);
1144 :
1145 6783813454 : const numeric_index_type local_index = this->map_global_to_local_index(i);
1146 :
1147 : #ifndef NDEBUG
1148 529618523 : if (this->is_effectively_ghosted())
1149 525841875 : libmesh_assert_less (local_index, _local_size);
1150 : #endif
1151 :
1152 6783813454 : return static_cast<T>(_read_only_values[local_index]);
1153 : }
1154 :
1155 :
1156 :
1157 : template <typename T>
1158 : inline
1159 92409442 : void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
1160 : T * values) const
1161 : {
1162 92409442 : this->_get_array(true);
1163 :
1164 17171828 : const std::size_t num = index.size();
1165 :
1166 572962216 : for (std::size_t i=0; i<num; i++)
1167 : {
1168 524587774 : const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1169 : #ifndef NDEBUG
1170 44037813 : if (this->is_effectively_ghosted())
1171 43867417 : libmesh_assert_less (local_index, _local_size);
1172 : #endif
1173 480552774 : values[i] = static_cast<T>(_read_only_values[local_index]);
1174 : }
1175 92409442 : }
1176 :
1177 :
1178 : template <typename T>
1179 : inline
1180 0 : PetscScalar * PetscVector<T>::get_array()
1181 : {
1182 0 : _get_array(false);
1183 0 : _values_manually_retrieved = true;
1184 :
1185 0 : return _values;
1186 : }
1187 :
1188 :
1189 : template <typename T>
1190 : inline
1191 0 : const PetscScalar * PetscVector<T>::get_array_read() const
1192 : {
1193 0 : _get_array(true);
1194 0 : _values_manually_retrieved = true;
1195 :
1196 0 : return _read_only_values;
1197 : }
1198 :
1199 : template <typename T>
1200 : inline
1201 0 : void PetscVector<T>::restore_array()
1202 : {
1203 : // \note \p _values_manually_retrieved needs to be set to \p false
1204 : // \e before calling \p _restore_array()!
1205 0 : _values_manually_retrieved = false;
1206 0 : _restore_array();
1207 0 : }
1208 :
1209 : template <typename T>
1210 : inline
1211 0 : Real PetscVector<T>::min () const
1212 : {
1213 0 : parallel_object_only();
1214 :
1215 0 : this->_restore_array();
1216 :
1217 0 : PetscInt index=0;
1218 0 : PetscReal returnval=0.;
1219 :
1220 0 : LibmeshPetscCall(VecMin (_vec, &index, &returnval));
1221 :
1222 : // this return value is correct: VecMin returns a PetscReal
1223 0 : return static_cast<Real>(returnval);
1224 : }
1225 :
1226 :
1227 :
1228 : template <typename T>
1229 : inline
1230 0 : Real PetscVector<T>::max() const
1231 : {
1232 0 : parallel_object_only();
1233 :
1234 0 : this->_restore_array();
1235 :
1236 0 : PetscInt index=0;
1237 0 : PetscReal returnval=0.;
1238 :
1239 0 : LibmeshPetscCall(VecMax (_vec, &index, &returnval));
1240 :
1241 : // this return value is correct: VecMax returns a PetscReal
1242 0 : return static_cast<Real>(returnval);
1243 : }
1244 :
1245 :
1246 :
1247 : template <typename T>
1248 : inline
1249 1328188 : void PetscVector<T>::swap (NumericVector<T> & other)
1250 : {
1251 28720 : parallel_object_only();
1252 :
1253 28720 : NumericVector<T>::swap(other);
1254 :
1255 28720 : PetscVector<T> & v = cast_ref<PetscVector<T> &>(other);
1256 :
1257 28720 : std::swap(_vec, v._vec);
1258 28720 : std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1259 28720 : std::swap(_global_to_local_map, v._global_to_local_map);
1260 :
1261 : #ifdef LIBMESH_HAVE_CXX11_THREAD
1262 : // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1263 1355276 : _array_is_present = v._array_is_present.exchange(_array_is_present);
1264 : #else
1265 : std::swap(_array_is_present, v._array_is_present);
1266 : #endif
1267 :
1268 28720 : std::swap(_local_form, v._local_form);
1269 28720 : std::swap(_values, v._values);
1270 1328188 : }
1271 :
1272 :
1273 :
1274 : template <typename T>
1275 : inline
1276 28565 : std::size_t PetscVector<T>::max_allowed_id () const
1277 : {
1278 : // The PetscInt type is used for indexing, it may be either a signed
1279 : // 4-byte or 8-byte integer depending on how PETSc is configured.
1280 28565 : return std::numeric_limits<PetscInt>::max();
1281 : }
1282 :
1283 :
1284 :
1285 : #ifdef LIBMESH_HAVE_CXX11
1286 : static_assert(sizeof(PetscInt) == sizeof(numeric_index_type),
1287 : "PETSc and libMesh integer sizes must match!");
1288 : #endif
1289 :
1290 :
1291 : inline
1292 7067359 : PetscInt * numeric_petsc_cast(const numeric_index_type * p)
1293 : {
1294 7067359 : return reinterpret_cast<PetscInt *>(const_cast<numeric_index_type *>(p));
1295 : }
1296 :
1297 : } // namespace libMesh
1298 :
1299 :
1300 : #endif // #ifdef LIBMESH_HAVE_PETSC
1301 : #endif // LIBMESH_PETSC_VECTOR_H
|
