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 552766 : Vec vec () { libmesh_assert (_vec); return _vec; }
346 :
347 4370946 : 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 2224533 : PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in, const ParallelType ptype) :
486 : NumericVector<T>(comm_in, ptype),
487 2105327 : _vec(nullptr),
488 : _array_is_present(false),
489 2105327 : _first(0),
490 2105327 : _last(0),
491 2105327 : _local_size(0),
492 2105327 : _local_form(nullptr),
493 2105327 : _read_only_values(nullptr),
494 2105327 : _values(nullptr),
495 2105327 : _global_to_local_map(),
496 2105327 : _destroy_vec_on_exit(true),
497 2105327 : _values_manually_retrieved(false),
498 2224533 : _values_read_only(false)
499 : {
500 2158279 : this->_type = ptype;
501 2224533 : }
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 989478 : PetscVector<T>::PetscVector (Vec v,
585 : const Parallel::Communicator & comm_in) :
586 : NumericVector<T>(comm_in, AUTOMATIC),
587 948518 : _vec(v),
588 : _array_is_present(false),
589 948518 : _first(0),
590 948518 : _last(0),
591 948518 : _local_size(0),
592 948518 : _local_form(nullptr),
593 948518 : _read_only_values(nullptr),
594 948518 : _values(nullptr),
595 948518 : _global_to_local_map(),
596 948518 : _destroy_vec_on_exit(false),
597 948518 : _values_manually_retrieved(false),
598 989478 : _values_read_only(false)
599 : {
600 989478 : this->_is_closed = true;
601 989478 : 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 989478 : PetscInt petsc_local_size=0;
606 989478 : LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_local_size));
607 :
608 : // Get the vector type from PETSc.
609 : VecType ptype;
610 989478 : 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 : else
639 : this->_type = PARALLEL;
640 : }
641 : else if (std::strcmp(ptype,VECSHARED) == 0)
642 : #else // PETSc < 3.21.0
643 989478 : if ((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
644 : #endif
645 : {
646 : ISLocalToGlobalMapping mapping;
647 974934 : LibmeshPetscCall(VecGetLocalToGlobalMapping(_vec, &mapping));
648 :
649 : Vec localrep;
650 974934 : LibmeshPetscCall(VecGhostGetLocalForm(_vec,&localrep));
651 : // If is a sparsely stored vector, set up our new mapping
652 : // Only checking mapping is not enough to determine if a vec is ghosted
653 : // We need to check if vec has a local representation
654 974934 : if (mapping && localrep)
655 : {
656 0 : const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
657 0 : const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
658 : PetscInt n;
659 0 : LibmeshPetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
660 :
661 0 : const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
662 : const PetscInt * indices;
663 0 : LibmeshPetscCall(ISLocalToGlobalMappingGetIndices(mapping,&indices));
664 :
665 0 : for (numeric_index_type i=ghost_begin; i<ghost_end; i++)
666 0 : _global_to_local_map[indices[i]] = i-my_local_size;
667 0 : this->_type = GHOSTED;
668 0 : LibmeshPetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
669 0 : }
670 : else
671 974934 : this->_type = PARALLEL;
672 :
673 974934 : LibmeshPetscCall(VecGhostRestoreLocalForm(_vec,&localrep));
674 933974 : }
675 : else
676 14544 : this->_type = SERIAL;
677 989478 : }
678 :
679 :
680 :
681 :
682 : template <typename T>
683 : inline
684 4506168 : PetscVector<T>::~PetscVector ()
685 : {
686 2751169 : this->clear ();
687 4578868 : }
688 :
689 :
690 :
691 : template <typename T>
692 : inline
693 914485 : void PetscVector<T>::init (const numeric_index_type n,
694 : const numeric_index_type n_local,
695 : const bool fast,
696 : const ParallelType ptype)
697 : {
698 27300 : parallel_object_only();
699 :
700 914485 : PetscInt petsc_n=static_cast<PetscInt>(n);
701 :
702 : // Clear initialized vectors
703 914485 : if (this->initialized())
704 26764 : this->clear();
705 :
706 914485 : if (ptype == AUTOMATIC)
707 : {
708 1190 : if (n == n_local)
709 569 : this->_type = SERIAL;
710 : else
711 621 : this->_type = PARALLEL;
712 : }
713 : else
714 913295 : this->_type = ptype;
715 :
716 27300 : libmesh_assert ((this->_type==SERIAL && n==n_local) ||
717 : this->_type==PARALLEL);
718 :
719 : // create a sequential vector if on only 1 processor
720 914485 : if (this->_type == SERIAL)
721 : {
722 2159 : LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
723 2159 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
724 2159 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetFromOptions (_vec));
725 : }
726 : // otherwise create an MPI-enabled vector
727 912326 : else if (this->_type == PARALLEL)
728 : {
729 : #ifdef LIBMESH_HAVE_MPI
730 27250 : PetscInt petsc_n_local=cast_int<PetscInt>(n_local);
731 27250 : libmesh_assert_less_equal (n_local, n);
732 : // Use more generic function instead of VecCreateSeq/MPI
733 912326 : LibmeshPetscCall(VecCreate(this->comm().get(), &_vec));
734 912326 : LibmeshPetscCall(VecSetSizes(_vec, petsc_n_local, petsc_n));
735 : #else
736 : libmesh_assert_equal_to (n_local, n);
737 : LibmeshPetscCallA(PETSC_COMM_SELF, VecCreate(PETSC_COMM_SELF, &_vec));
738 : LibmeshPetscCallA(PETSC_COMM_SELF, VecSetSizes(_vec, petsc_n, petsc_n));
739 : #endif
740 912326 : LibmeshPetscCall(VecSetFromOptions(_vec));
741 : }
742 : else
743 0 : libmesh_error_msg("Unsupported type " << Utility::enum_to_string(this->_type));
744 :
745 914485 : this->_is_initialized = true;
746 914485 : this->_is_closed = true;
747 :
748 :
749 914485 : if (fast == false)
750 905637 : this->zero ();
751 914485 : }
752 :
753 :
754 :
755 : template <typename T>
756 : inline
757 1606 : void PetscVector<T>::init (const numeric_index_type n,
758 : const bool fast,
759 : const ParallelType ptype)
760 : {
761 2150 : this->init(n,n,fast,ptype);
762 2134 : }
763 :
764 :
765 :
766 : template <typename T>
767 : inline
768 508076 : void PetscVector<T>::init (const numeric_index_type n,
769 : const numeric_index_type n_local,
770 : const std::vector<numeric_index_type> & ghost,
771 : const bool fast,
772 : const ParallelType libmesh_dbg_var(ptype))
773 : {
774 16244 : parallel_object_only();
775 :
776 508076 : PetscInt petsc_n=static_cast<PetscInt>(n);
777 508076 : PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
778 32300 : PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
779 :
780 : // If the mesh is not disjoint, every processor will either have
781 : // all the dofs, none of the dofs, or some non-zero dofs at the
782 : // boundary between processors.
783 : //
784 : // However we can't assert this, because someone might want to
785 : // construct a GHOSTED vector which doesn't include neighbor element
786 : // dofs. Boyce tried to do so in user code, and we're going to want
787 : // to do so in System::project_vector().
788 : //
789 : // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
790 :
791 : libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
792 :
793 508076 : PetscInt * petsc_ghost = ghost.empty() ? LIBMESH_PETSC_NULLPTR :
794 15855 : const_cast<PetscInt *>(reinterpret_cast<const PetscInt *>(ghost.data()));
795 :
796 : // Clear initialized vectors
797 508076 : if (this->initialized())
798 20692 : this->clear();
799 :
800 16244 : libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
801 508076 : this->_type = GHOSTED;
802 :
803 : /* Make the global-to-local ghost cell map. */
804 35173354 : for (auto i : index_range(ghost))
805 : {
806 34665278 : _global_to_local_map[ghost[i]] = i;
807 : }
808 :
809 : /* Create vector. */
810 508076 : LibmeshPetscCall(VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
811 : petsc_n_ghost, petsc_ghost,
812 : &_vec));
813 :
814 : // Add a prefix so that we can at least distinguish a ghosted vector from a
815 : // nonghosted vector when using a petsc option.
816 : // PETSc does not fully support VecGhost on GPU yet. This change allows us to
817 : // trigger a nonghosted vector to use GPU without bothering the ghosted vectors.
818 508076 : LibmeshPetscCall(PetscObjectAppendOptionsPrefix((PetscObject)_vec,"ghost_"));
819 :
820 508076 : LibmeshPetscCall(VecSetFromOptions (_vec));
821 :
822 508076 : this->_is_initialized = true;
823 508076 : this->_is_closed = true;
824 :
825 508076 : if (fast == false)
826 157395 : this->zero ();
827 508076 : }
828 :
829 :
830 :
831 : template <typename T>
832 : inline
833 424818 : void PetscVector<T>::init (const NumericVector<T> & other,
834 : const bool fast)
835 : {
836 12538 : parallel_object_only();
837 :
838 : // Clear initialized vectors
839 424818 : if (this->initialized())
840 32690 : this->clear();
841 :
842 12538 : const PetscVector<T> & v = cast_ref<const PetscVector<T> &>(other);
843 :
844 : // Other vector should restore array.
845 424818 : if (v.initialized())
846 : {
847 424818 : v._restore_array();
848 : }
849 :
850 12538 : this->_global_to_local_map = v._global_to_local_map;
851 :
852 : // Even if we're initializing sizes based on an uninitialized or
853 : // unclosed vector, *this* vector is being initialized now and is
854 : // initially closed.
855 424818 : this->_is_closed = true; // v._is_closed;
856 424818 : this->_is_initialized = true; // v._is_initialized;
857 :
858 424818 : this->_type = v._type;
859 :
860 : // We want to have a valid Vec, even if it's initially of size zero
861 424818 : LibmeshPetscCall(VecDuplicate (v._vec, &this->_vec));
862 :
863 424818 : if (fast == false)
864 53720 : this->zero ();
865 424818 : }
866 :
867 :
868 :
869 : template <typename T>
870 : inline
871 3572799 : void PetscVector<T>::close ()
872 : {
873 97244 : parallel_object_only();
874 :
875 3572799 : this->_restore_array();
876 :
877 3572799 : VecAssemblyBeginEnd(this->comm(), _vec);
878 :
879 3572799 : if (this->type() == GHOSTED)
880 898621 : VecGhostUpdateBeginEnd(this->comm(), _vec, INSERT_VALUES, SCATTER_FORWARD);
881 :
882 3572799 : this->_is_closed = true;
883 3572799 : }
884 :
885 :
886 :
887 : template <typename T>
888 : inline
889 2997247 : void PetscVector<T>::clear () noexcept
890 : {
891 81766 : exceptionless_parallel_object_only();
892 :
893 2997247 : if (this->initialized())
894 2836857 : this->_restore_array();
895 :
896 2997247 : if ((this->initialized()) && (this->_destroy_vec_on_exit))
897 : {
898 : // If we encounter an error here, print a warning but otherwise
899 : // keep going since we may be recovering from an exception.
900 1847379 : PetscErrorCode ierr = VecDestroy(&_vec);
901 : if (ierr)
902 : libmesh_warning("Warning: VecDestroy returned a non-zero error code which we ignored.");
903 : }
904 :
905 2997247 : this->_is_closed = this->_is_initialized = false;
906 :
907 81766 : _global_to_local_map.clear();
908 2997247 : }
909 :
910 :
911 :
912 : template <typename T>
913 : inline
914 5986454 : void PetscVector<T>::zero ()
915 : {
916 171112 : parallel_object_only();
917 :
918 171112 : libmesh_assert(this->closed());
919 :
920 5986454 : this->_restore_array();
921 :
922 171112 : PetscScalar z=0.;
923 :
924 5986454 : if (this->type() != GHOSTED)
925 5774439 : LibmeshPetscCall(VecSet (_vec, z));
926 : else
927 : {
928 : /* Vectors that include ghost values require a special
929 : handling. */
930 : Vec loc_vec;
931 212015 : LibmeshPetscCall(VecGhostGetLocalForm (_vec,&loc_vec));
932 212015 : LibmeshPetscCall(VecSet (loc_vec, z));
933 212015 : LibmeshPetscCall(VecGhostRestoreLocalForm (_vec, &loc_vec));
934 : }
935 5986454 : }
936 :
937 :
938 :
939 : template <typename T>
940 : inline
941 21030 : std::unique_ptr<NumericVector<T>> PetscVector<T>::zero_clone () const
942 : {
943 21030 : std::unique_ptr<NumericVector<T>> cloned_vector =
944 20398 : std::make_unique<PetscVector<T>>(this->comm(), this->type());
945 21030 : cloned_vector->init(*this);
946 21030 : return cloned_vector;
947 0 : }
948 :
949 :
950 :
951 : template <typename T>
952 : inline
953 337988 : std::unique_ptr<NumericVector<T>> PetscVector<T>::clone () const
954 : {
955 337988 : std::unique_ptr<NumericVector<T>> cloned_vector =
956 327962 : std::make_unique<PetscVector<T>>(this->comm(), this->type());
957 337988 : cloned_vector->init(*this, true);
958 337988 : *cloned_vector = *this;
959 337988 : return cloned_vector;
960 0 : }
961 :
962 :
963 :
964 : template <typename T>
965 : inline
966 183588533 : numeric_index_type PetscVector<T>::size () const
967 : {
968 183009556 : libmesh_assert (this->initialized());
969 :
970 183588533 : PetscInt petsc_size=0;
971 :
972 183588533 : if (!this->initialized())
973 0 : return 0;
974 :
975 183588533 : LibmeshPetscCall(VecGetSize(_vec, &petsc_size));
976 :
977 183588533 : return static_cast<numeric_index_type>(petsc_size);
978 : }
979 :
980 :
981 :
982 : template <typename T>
983 : inline
984 1021190 : numeric_index_type PetscVector<T>::local_size () const
985 : {
986 212356 : libmesh_assert (this->initialized());
987 :
988 1021190 : PetscInt petsc_size=0;
989 :
990 1021190 : LibmeshPetscCall(VecGetLocalSize(_vec, &petsc_size));
991 :
992 1021190 : return static_cast<numeric_index_type>(petsc_size);
993 : }
994 :
995 :
996 :
997 : template <typename T>
998 : inline
999 32403036 : numeric_index_type PetscVector<T>::first_local_index () const
1000 : {
1001 4316285 : libmesh_assert (this->initialized());
1002 :
1003 4316285 : numeric_index_type first = 0;
1004 :
1005 32403036 : if (_array_is_present) // Can we use cached values?
1006 1568641 : first = _first;
1007 : else
1008 : {
1009 30834395 : PetscInt petsc_first=0, petsc_last=0;
1010 30834395 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1011 30834395 : first = static_cast<numeric_index_type>(petsc_first);
1012 : }
1013 :
1014 32403036 : return first;
1015 : }
1016 :
1017 :
1018 :
1019 : template <typename T>
1020 : inline
1021 6315554 : numeric_index_type PetscVector<T>::last_local_index () const
1022 : {
1023 3247705 : libmesh_assert (this->initialized());
1024 :
1025 3247705 : numeric_index_type last = 0;
1026 :
1027 6315554 : if (_array_is_present) // Can we use cached values?
1028 1450354 : last = _last;
1029 : else
1030 : {
1031 4865200 : PetscInt petsc_first=0, petsc_last=0;
1032 4865200 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1033 4865200 : last = static_cast<numeric_index_type>(petsc_last);
1034 : }
1035 :
1036 6315554 : return last;
1037 : }
1038 :
1039 :
1040 :
1041 : template <typename T>
1042 : inline
1043 6965134839 : numeric_index_type PetscVector<T>::map_global_to_local_index (const numeric_index_type i) const
1044 : {
1045 551864446 : libmesh_assert (this->initialized());
1046 :
1047 551864446 : numeric_index_type first=0;
1048 551864446 : numeric_index_type last=0;
1049 :
1050 6965134839 : if (_array_is_present) // Can we use cached values?
1051 : {
1052 6965134839 : first = _first;
1053 6965134839 : last = _last;
1054 : }
1055 : else
1056 : {
1057 0 : PetscInt petsc_first=0, petsc_last=0;
1058 0 : LibmeshPetscCall(VecGetOwnershipRange (_vec, &petsc_first, &petsc_last));
1059 0 : first = static_cast<numeric_index_type>(petsc_first);
1060 0 : last = static_cast<numeric_index_type>(petsc_last);
1061 : }
1062 :
1063 :
1064 6965134839 : if ((i>=first) && (i<last))
1065 : {
1066 6797120963 : return i-first;
1067 : }
1068 :
1069 5089852 : GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1070 : #ifndef NDEBUG
1071 5089852 : const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1072 5089852 : if (it == end)
1073 : {
1074 0 : std::ostringstream error_message;
1075 0 : error_message << "No index " << i << " in ghosted vector.\n"
1076 0 : << "Vector contains [" << first << ',' << last << ")\n";
1077 0 : GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1078 0 : if (b == end)
1079 : {
1080 0 : error_message << "And empty ghost array.\n";
1081 : }
1082 : else
1083 : {
1084 0 : error_message << "And ghost array {" << b->first;
1085 0 : for (++b; b != end; ++b)
1086 0 : error_message << ',' << b->first;
1087 0 : error_message << "}\n";
1088 : }
1089 :
1090 0 : libmesh_error_msg(error_message.str());
1091 : }
1092 5089852 : libmesh_assert (it != _global_to_local_map.end());
1093 : #endif
1094 168013876 : return it->second+last-first;
1095 : }
1096 :
1097 :
1098 :
1099 : template <typename T>
1100 : inline
1101 6508448388 : T PetscVector<T>::operator() (const numeric_index_type i) const
1102 : {
1103 6508448388 : this->_get_array(true);
1104 :
1105 6508448388 : const numeric_index_type local_index = this->map_global_to_local_index(i);
1106 :
1107 : #ifndef NDEBUG
1108 509714296 : if (this->type() == GHOSTED)
1109 : {
1110 505941494 : libmesh_assert_less (local_index, _local_size);
1111 : }
1112 : #endif
1113 :
1114 6508448388 : return static_cast<T>(_read_only_values[local_index]);
1115 : }
1116 :
1117 :
1118 :
1119 : template <typename T>
1120 : inline
1121 88410034 : void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
1122 : T * values) const
1123 : {
1124 88410034 : this->_get_array(true);
1125 :
1126 16571813 : const std::size_t num = index.size();
1127 :
1128 545096485 : for (std::size_t i=0; i<num; i++)
1129 : {
1130 498838308 : const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1131 : #ifndef NDEBUG
1132 42150150 : if (this->type() == GHOSTED)
1133 : {
1134 42140562 : libmesh_assert_less (local_index, _local_size);
1135 : }
1136 : #endif
1137 456686451 : values[i] = static_cast<T>(_read_only_values[local_index]);
1138 : }
1139 88410034 : }
1140 :
1141 :
1142 : template <typename T>
1143 : inline
1144 0 : PetscScalar * PetscVector<T>::get_array()
1145 : {
1146 0 : _get_array(false);
1147 0 : _values_manually_retrieved = true;
1148 :
1149 0 : return _values;
1150 : }
1151 :
1152 :
1153 : template <typename T>
1154 : inline
1155 0 : const PetscScalar * PetscVector<T>::get_array_read() const
1156 : {
1157 0 : _get_array(true);
1158 0 : _values_manually_retrieved = true;
1159 :
1160 0 : return _read_only_values;
1161 : }
1162 :
1163 : template <typename T>
1164 : inline
1165 0 : void PetscVector<T>::restore_array()
1166 : {
1167 : // \note \p _values_manually_retrieved needs to be set to \p false
1168 : // \e before calling \p _restore_array()!
1169 0 : _values_manually_retrieved = false;
1170 0 : _restore_array();
1171 0 : }
1172 :
1173 : template <typename T>
1174 : inline
1175 0 : Real PetscVector<T>::min () const
1176 : {
1177 0 : parallel_object_only();
1178 :
1179 0 : this->_restore_array();
1180 :
1181 0 : PetscInt index=0;
1182 0 : PetscReal returnval=0.;
1183 :
1184 0 : LibmeshPetscCall(VecMin (_vec, &index, &returnval));
1185 :
1186 : // this return value is correct: VecMin returns a PetscReal
1187 0 : return static_cast<Real>(returnval);
1188 : }
1189 :
1190 :
1191 :
1192 : template <typename T>
1193 : inline
1194 0 : Real PetscVector<T>::max() const
1195 : {
1196 0 : parallel_object_only();
1197 :
1198 0 : this->_restore_array();
1199 :
1200 0 : PetscInt index=0;
1201 0 : PetscReal returnval=0.;
1202 :
1203 0 : LibmeshPetscCall(VecMax (_vec, &index, &returnval));
1204 :
1205 : // this return value is correct: VecMax returns a PetscReal
1206 0 : return static_cast<Real>(returnval);
1207 : }
1208 :
1209 :
1210 :
1211 : template <typename T>
1212 : inline
1213 1328188 : void PetscVector<T>::swap (NumericVector<T> & other)
1214 : {
1215 28720 : parallel_object_only();
1216 :
1217 28720 : NumericVector<T>::swap(other);
1218 :
1219 28720 : PetscVector<T> & v = cast_ref<PetscVector<T> &>(other);
1220 :
1221 28720 : std::swap(_vec, v._vec);
1222 28720 : std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1223 28720 : std::swap(_global_to_local_map, v._global_to_local_map);
1224 :
1225 : #ifdef LIBMESH_HAVE_CXX11_THREAD
1226 : // Only truly atomic for v... but swap() doesn't really need to be thread safe!
1227 1355276 : _array_is_present = v._array_is_present.exchange(_array_is_present);
1228 : #else
1229 : std::swap(_array_is_present, v._array_is_present);
1230 : #endif
1231 :
1232 28720 : std::swap(_local_form, v._local_form);
1233 28720 : std::swap(_values, v._values);
1234 1328188 : }
1235 :
1236 :
1237 :
1238 : template <typename T>
1239 : inline
1240 27455 : std::size_t PetscVector<T>::max_allowed_id () const
1241 : {
1242 : // The PetscInt type is used for indexing, it may be either a signed
1243 : // 4-byte or 8-byte integer depending on how PETSc is configured.
1244 27455 : return std::numeric_limits<PetscInt>::max();
1245 : }
1246 :
1247 :
1248 :
1249 : #ifdef LIBMESH_HAVE_CXX11
1250 : static_assert(sizeof(PetscInt) == sizeof(numeric_index_type),
1251 : "PETSc and libMesh integer sizes must match!");
1252 : #endif
1253 :
1254 :
1255 : inline
1256 6996138 : PetscInt * numeric_petsc_cast(const numeric_index_type * p)
1257 : {
1258 6996138 : return reinterpret_cast<PetscInt *>(const_cast<numeric_index_type *>(p));
1259 : }
1260 :
1261 : } // namespace libMesh
1262 :
1263 :
1264 : #endif // #ifdef LIBMESH_HAVE_PETSC
1265 : #endif // LIBMESH_PETSC_VECTOR_H
|
