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 : #ifndef LIBMESH_NUMERIC_VECTOR_H
21 : #define LIBMESH_NUMERIC_VECTOR_H
22 :
23 : // Local includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/enum_parallel_type.h" // AUTOMATIC
26 : #include "libmesh/id_types.h"
27 : #include "libmesh/int_range.h"
28 : #include "libmesh/reference_counted_object.h"
29 : #include "libmesh/libmesh.h"
30 : #include "libmesh/parallel_object.h"
31 : #include "libmesh/dense_subvector.h"
32 : #include "libmesh/dense_vector.h"
33 :
34 : // C++ includes
35 : #include <cstddef>
36 : #include <set>
37 : #include <vector>
38 : #include <memory>
39 : #include <mutex>
40 :
41 : namespace libMesh
42 : {
43 :
44 :
45 : // forward declarations
46 : template <typename T> class NumericVector;
47 : template <typename T> class DenseVector;
48 : template <typename T> class DenseSubVector;
49 : template <typename T> class SparseMatrix;
50 : template <typename T> class ShellMatrix;
51 : enum SolverPackage : int;
52 :
53 : /**
54 : * \brief Provides a uniform interface to vector storage schemes for different
55 : * linear algebra libraries.
56 : *
57 : * \note This class is the abstract base class for different implementations
58 : * of numeric vectors. Most of the time you should use a System object to
59 : * create numeric vectors. If this is not desired, you can instantiate one of
60 : * the derived classes (PetscVector, EigenSparseVector, etc.) or use the
61 : * NumericVector::build method. When creating the vector yourself, make sure
62 : * that you initialize the vector properly (NumericVector::init).
63 : *
64 : * \author Benjamin S. Kirk
65 : * \date 2003
66 : */
67 : template <typename T>
68 : class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
69 : public ParallelObject
70 : {
71 : public:
72 :
73 : /**
74 : * Dummy-Constructor. Dimension=0
75 : */
76 : explicit
77 : NumericVector (const Parallel::Communicator & comm_in,
78 : const ParallelType ptype = AUTOMATIC);
79 :
80 : /**
81 : * Constructor. Set dimension to \p n and initialize all elements with zero.
82 : */
83 : explicit
84 : NumericVector (const Parallel::Communicator & comm_in,
85 : const numeric_index_type n,
86 : const ParallelType ptype = AUTOMATIC);
87 :
88 : /**
89 : * Constructor. Set local dimension to \p n_local, the global dimension
90 : * to \p n, and initialize all elements with zero.
91 : */
92 : NumericVector (const Parallel::Communicator & comm_in,
93 : const numeric_index_type n,
94 : const numeric_index_type n_local,
95 : const ParallelType ptype = AUTOMATIC);
96 :
97 : /**
98 : * Constructor. Set local dimension to \p n_local, the global
99 : * dimension to \p n, but additionally reserve memory for the
100 : * indices specified by the \p ghost argument.
101 : */
102 : NumericVector (const Parallel::Communicator & comm_in,
103 : const numeric_index_type N,
104 : const numeric_index_type n_local,
105 : const std::vector<numeric_index_type> & ghost,
106 : const ParallelType ptype = AUTOMATIC);
107 :
108 : /**
109 : * This _looks_ like a copy assignment operator, but note that,
110 : * unlike normal copy assignment operators, it is pure virtual. This
111 : * function should be overridden in derived classes so that they can
112 : * be copied correctly via references to the base class. This design
113 : * usually isn't a good idea in general, but in this context it
114 : * works because we usually don't have a mix of different kinds of
115 : * NumericVectors active in the library at a single time.
116 : *
117 : * \returns A reference to *this as the base type.
118 : */
119 : virtual NumericVector<T> & operator= (const NumericVector<T> & v) = 0;
120 :
121 : /**
122 : * These 3 special functions can be defaulted for this class, as it
123 : * does not manage any memory itself.
124 : */
125 : NumericVector (NumericVector &&) = default;
126 : NumericVector (const NumericVector &) = default;
127 : NumericVector & operator= (NumericVector &&) = default;
128 :
129 : /**
130 : * While this class doesn't manage any memory, the derived class might and
131 : * users may be deleting through a pointer to this base class
132 : */
133 76280 : virtual ~NumericVector() = default;
134 :
135 : /**
136 : * Builds a \p NumericVector on the processors in communicator
137 : * \p comm using the linear solver package specified by
138 : * \p solver_package
139 : */
140 : static std::unique_ptr<NumericVector<T>>
141 : build(const Parallel::Communicator & comm,
142 : SolverPackage solver_package = libMesh::default_solver_package(),
143 : ParallelType parallel_type = AUTOMATIC);
144 :
145 : /**
146 : * \returns \p true if the vector has been initialized,
147 : * false otherwise.
148 : */
149 1532130191 : virtual bool initialized() const { return _is_initialized; }
150 :
151 : /**
152 : * \returns The type (SERIAL, PARALLEL, GHOSTED) of the vector.
153 : *
154 : * This is metadata reflecting what type of vector was requested,
155 : * which should reflect what level of replication and ghosting is
156 : * necessary in parallel.
157 : */
158 2809389 : ParallelType type() const { return _type; }
159 :
160 : /**
161 : * \returns Whether the vector is effectively serial, i.e. whether
162 : * all vector data is accessible on every processor.
163 : *
164 : * This is true for explicitly SERIAL vectors, but also for vectors
165 : * on serial (1-processor) communicators.
166 : */
167 :
168 329670 : bool is_effectively_serial() const
169 4252501 : { return (this->n_processors()==1) || (_type == SERIAL); }
170 :
171 : /**
172 : * \returns Whether the vector is effectively ghosted, i.e. whether
173 : * vector operations require handling of ghosted coefficients.
174 : *
175 : * This is true for GHOSTED vectors in parallel, but in serial
176 : * subclasses may return false here for GHOSTED vectors since every
177 : * ghosted vector is effectively serial and thus may be able to take
178 : * advantage of serial-specific code paths.
179 : */
180 :
181 573492878 : bool is_effectively_ghosted() const
182 595574439 : { return (this->n_processors()!=1) && (_type == GHOSTED); }
183 :
184 : /**
185 : * \returns The type (SERIAL, PARALLEL, GHOSTED) of the vector.
186 : *
187 : * \deprecated because it is dangerous to change the ParallelType
188 : * of an already-initialized NumericVector. See NumericVector::set_type()
189 : * for a safer, non-deprecated setter.
190 : */
191 : #ifdef LIBMESH_ENABLE_DEPRECATED
192 460348 : ParallelType & type() { return _type; }
193 : #endif
194 :
195 : /**
196 : * Allow the user to change the ParallelType of the NumericVector
197 : * under some circumstances. If the NumericVector has not been
198 : * initialized yet, then it is generally safe to change the
199 : * ParallelType. otherwise, if the NumericVector has already been
200 : * initialized with a specific type, we cannot change it without
201 : * doing some extra copying/reinitialization work, and we currently
202 : * throw an error if this is requested.
203 : */
204 : void set_type(ParallelType t);
205 :
206 : /**
207 : * \returns \p true if the vector is closed and ready for
208 : * computation, false otherwise.
209 : */
210 974572 : virtual bool closed() const { return _is_closed; }
211 :
212 : /**
213 : * Calls the NumericVector's internal assembly routines, ensuring
214 : * that the values are consistent across processors.
215 : */
216 : virtual void close () = 0;
217 :
218 : /**
219 : * Restores the \p NumericVector<T> to a pristine state.
220 : */
221 : virtual void clear ();
222 :
223 : /**
224 : * Set all entries to zero. Equivalent to \p v = 0, but more obvious and
225 : * faster.
226 : */
227 : virtual void zero () = 0;
228 :
229 : /**
230 : * \returns A smart pointer to a copy of this vector with the same
231 : * type, size, and partitioning, but with all zero entries.
232 : *
233 : * \note This must be overridden in the derived classes.
234 : */
235 : virtual std::unique_ptr<NumericVector<T>> zero_clone () const = 0;
236 :
237 : /**
238 : * \returns A copy of this vector wrapped in a smart pointer.
239 : *
240 : * \note This must be overridden in the derived classes.
241 : */
242 : virtual std::unique_ptr<NumericVector<T>> clone () const = 0;
243 :
244 : /**
245 : * Change the dimension of the vector to \p n. The reserved memory
246 : * for this vector remains unchanged if possible. If \p n==0, all
247 : * memory is freed. Therefore, if you want to resize the vector and
248 : * release the memory not needed, you have to first call \p init(0)
249 : * and then \p init(n). This behaviour is analogous to that of the
250 : * STL containers.
251 : *
252 : * On \p fast==false, the vector is filled by zeros.
253 : */
254 : virtual void init (const numeric_index_type n,
255 : const numeric_index_type n_local,
256 : const bool fast = false,
257 : const ParallelType ptype = AUTOMATIC) = 0;
258 :
259 : /**
260 : * Call \p init() with n_local = N.
261 : */
262 : virtual void init (const numeric_index_type n,
263 : const bool fast = false,
264 : const ParallelType ptype = AUTOMATIC) = 0;
265 :
266 : /**
267 : * Create a vector that holds the local indices plus those specified
268 : * in the \p ghost argument.
269 : */
270 : virtual void init (const numeric_index_type n,
271 : const numeric_index_type n_local,
272 : const std::vector<numeric_index_type> & ghost,
273 : const bool fast = false,
274 : const ParallelType ptype = AUTOMATIC) = 0;
275 :
276 : /**
277 : * Creates a vector that has the same dimension and storage type as
278 : * \p other, including ghost dofs.
279 : */
280 : virtual void init (const NumericVector<T> & other,
281 : const bool fast = false) = 0;
282 :
283 : /**
284 : * Sets all entries of the vector to the value \p s.
285 : *
286 : * \returns A reference to *this.
287 : */
288 : virtual NumericVector<T> & operator= (const T s) = 0;
289 :
290 : /**
291 : * Sets (*this)(i) = v(i) for each entry of the vector.
292 : *
293 : * \returns A reference to *this as the base type.
294 : */
295 : virtual NumericVector<T> & operator= (const std::vector<T> & v) = 0;
296 :
297 : /**
298 : * \returns The minimum entry in the vector, or the minimum real
299 : * part in the case of complex numbers.
300 : */
301 : virtual Real min () const = 0;
302 :
303 : /**
304 : * \returns The maximum entry in the vector, or the maximum real
305 : * part in the case of complex numbers.
306 : */
307 : virtual Real max () const = 0;
308 :
309 : /**
310 : * \returns The sum of all values in the vector.
311 : */
312 : virtual T sum() const = 0;
313 :
314 : /**
315 : * \returns The \f$ \ell_1 \f$-norm of the vector, i.e. the sum of the
316 : * absolute values of the entries.
317 : */
318 : virtual Real l1_norm () const = 0;
319 :
320 : /**
321 : * \returns The \f$ \ell_2 \f$-norm of the vector, i.e. the square root
322 : * of the sum of the squares of the entries.
323 : */
324 : virtual Real l2_norm () const = 0;
325 :
326 : /**
327 : * \returns The \f$ \ell_{\infty} \f$-norm of the vector, i.e. the maximum
328 : * absolute value of the entries of the vector.
329 : */
330 : virtual Real linfty_norm () const = 0;
331 :
332 : /**
333 : * \returns The \f$ \ell_1 \f$-norm of the vector, i.e. the sum of the
334 : * absolute values for the specified entries in the vector.
335 : *
336 : * \note The indices must necessarily live on this processor.
337 : */
338 : virtual Real subset_l1_norm (const std::set<numeric_index_type> & indices) const;
339 :
340 : /**
341 : * \returns The \f$ \ell_2 \f$-norm of the vector, i.e. the square root
342 : * of the sum of the squares of the elements for the specified
343 : * entries in the vector.
344 : *
345 : * \note The indices must necessarily live on this processor.
346 : */
347 : virtual Real subset_l2_norm (const std::set<numeric_index_type> & indices) const;
348 :
349 : /**
350 : * \returns The maximum absolute value of the specified entries of
351 : * this vector, which is the \f$ \ell_{\infty} \f$-norm of a vector.
352 : *
353 : * \note The indices must necessarily live on this processor.
354 : */
355 : virtual Real subset_linfty_norm (const std::set<numeric_index_type> & indices) const;
356 :
357 : /**
358 : * \returns The \f$ \ell_2 \f$-norm of \f$ \vec{u} - \vec{v} \f$, where
359 : * \f$ \vec{u} \f$ is \p this.
360 : */
361 : Real l2_norm_diff (const NumericVector<T> & other_vec) const;
362 :
363 : /**
364 : * \returns The \f$ \ell_1 \f$-norm of \f$ \vec{u} - \vec{v} \f$, where
365 : * \f$ \vec{u} \f$ is \p this.
366 : */
367 : Real l1_norm_diff (const NumericVector<T> & other_vec) const;
368 :
369 : /**
370 : * \returns The size of the vector.
371 : */
372 : virtual numeric_index_type size () const = 0;
373 :
374 : /**
375 : * \returns The local size of the vector, i.e. \p index_stop - \p index_start.
376 : */
377 : virtual numeric_index_type local_size() const = 0;
378 :
379 : /**
380 : * \returns The index of the first vector element actually stored on
381 : * this processor.
382 : *
383 : * \note The minimum for this index is \p 0.
384 : */
385 : virtual numeric_index_type first_local_index() const = 0;
386 :
387 : /**
388 : * \returns The index+1 of the last vector element actually stored
389 : * on this processor.
390 : *
391 : * \note The maximum for this index is \p size().
392 : */
393 : virtual numeric_index_type last_local_index() const = 0;
394 :
395 : /**
396 : * \returns A copy of the ith entry of the vector.
397 : */
398 : virtual T operator() (const numeric_index_type i) const = 0;
399 :
400 : /**
401 : * \returns (*this)(i).
402 : */
403 0 : virtual T el(const numeric_index_type i) const { return (*this)(i); }
404 :
405 : /**
406 : * Access multiple components at once. \p values will *not* be
407 : * reallocated; it should already have enough space. The default
408 : * implementation calls \p operator() for each index, but some
409 : * implementations may supply faster methods here.
410 : */
411 : virtual void get(const std::vector<numeric_index_type> & index,
412 : T * values) const;
413 :
414 : /**
415 : * Access multiple components at once. \p values will be resized,
416 : * if necessary, and filled. The default implementation calls \p
417 : * operator() for each index, but some implementations may supply
418 : * faster methods here.
419 : */
420 : void get(const std::vector<numeric_index_type> & index,
421 : std::vector<T> & values) const;
422 :
423 : /**
424 : * Adds \p v to *this,
425 : * \f$ \vec{u} \leftarrow \vec{u} + \vec{v} \f$.
426 : * Equivalent to \p u.add(1, v).
427 : *
428 : * \returns A reference to *this.
429 : */
430 : virtual NumericVector<T> & operator += (const NumericVector<T> & v) = 0;
431 :
432 : /**
433 : * Subtracts \p v from *this,
434 : * \f$ \vec{u} \leftarrow \vec{u} - \vec{v} \f$.
435 : * Equivalent to \p u.add(-1, v).
436 : *
437 : * \returns A reference to *this.
438 : */
439 : virtual NumericVector<T> & operator -= (const NumericVector<T> & v) = 0;
440 :
441 : /**
442 : * Scales the vector by \p a,
443 : * \f$ \vec{u} \leftarrow a\vec{u} \f$.
444 : * Equivalent to \p u.scale(a)
445 : *
446 : * \returns A reference to *this.
447 : */
448 11572 : NumericVector<T> & operator *= (const T a) { this->scale(a); return *this; }
449 :
450 : /**
451 : * Computes the component-wise multiplication of this vector's entries by
452 : * another's,
453 : * \f$ u_i \leftarrow u_i v_i \, \forall i\f$
454 : *
455 : * \returns A reference to *this.
456 : */
457 : virtual NumericVector<T> & operator *= (const NumericVector<T> & v) = 0;
458 :
459 : /**
460 : * Scales the vector by \p 1/a,
461 : * \f$ \vec{u} \leftarrow \frac{1}{a}\vec{u} \f$.
462 : * Equivalent to \p u.scale(1./a)
463 : *
464 : * \returns A reference to *this.
465 : */
466 7430 : NumericVector<T> & operator /= (const T a) { this->scale(1./a); return *this; }
467 :
468 : /**
469 : * Computes the component-wise division of this vector's entries by another's,
470 : * \f$ u_i \leftarrow \frac{u_i}{v_i} \, \forall i\f$
471 : *
472 : * \returns A reference to *this.
473 : */
474 : virtual NumericVector<T> & operator /= (const NumericVector<T> & v) = 0;
475 :
476 : /**
477 : * Computes the component-wise reciprocal,
478 : * \f$ u_i \leftarrow \frac{1}{u_i} \, \forall i\f$
479 : */
480 : virtual void reciprocal() = 0;
481 :
482 : /**
483 : * Negates the imaginary component of each entry in the vector.
484 : */
485 : virtual void conjugate() = 0;
486 :
487 : /**
488 : * Sets v(i) = \p value. Note that library implementations of this method are
489 : * thread safe, e.g. we will lock \p _numeric_vector_mutex before writing to the vector
490 : */
491 : virtual void set (const numeric_index_type i, const T value) = 0;
492 :
493 : /**
494 : * Adds \p value to the vector entry specified by \p i. Note that library implementations of this
495 : * method are thread safe, e.g. we will lock \p _numeric_vector_mutex before writing to the vector
496 : */
497 : virtual void add (const numeric_index_type i, const T value) = 0;
498 :
499 : /**
500 : * Adds \p s to each entry of the vector,
501 : * \f$ u_i \leftarrow u_i + s \f$
502 : */
503 : virtual void add (const T s) = 0;
504 :
505 : /**
506 : * Adds \p v to \p this,
507 : * \f$ \vec{u} \leftarrow \vec{u} + \vec{v} \f$.
508 : * Equivalent to calling \p operator+=().
509 : */
510 : virtual void add (const NumericVector<T> & v) = 0;
511 :
512 : /**
513 : * Vector addition with a scalar multiple,
514 : * \f$ \vec{u} \leftarrow \vec{u} + a\vec{v} \f$.
515 : * Equivalent to calling \p operator+=().
516 : */
517 : virtual void add (const T a, const NumericVector<T> & v) = 0;
518 :
519 : /**
520 : * Computes \f$ \vec{u} \leftarrow \vec{u} + \vec{v} \f$,
521 : * where \p v is a pointer and each \p dof_indices[i] specifies where
522 : * to add value \p v[i]. This should be overridden in subclasses
523 : * for efficiency. Note that library implementations of this
524 : * method are thread safe
525 : */
526 : virtual void add_vector (const T * v,
527 : const std::vector<numeric_index_type> & dof_indices);
528 :
529 : /**
530 : * Computes \f$ \vec{u} \leftarrow \vec{u} + \vec{v} \f$,
531 : * where \p v is a std::vector and each \p dof_indices[i] specifies
532 : * where to add value \p v[i]. This method is guaranteed to be thread safe as long as library
533 : * implementations of \p NumericVector are used
534 : */
535 : void add_vector (const std::vector<T> & v,
536 : const std::vector<numeric_index_type> & dof_indices);
537 :
538 : /**
539 : * Computes \f$ \vec{u} \leftarrow \vec{u} + \vec{v} \f$,
540 : * where \p v is a NumericVector and each \p dof_indices[i]
541 : * specifies where to add value \p v(i). This method is
542 : * guaranteed to be thread-safe as long as library implementations of \p NumericVector are used
543 : */
544 : void add_vector (const NumericVector<T> & v,
545 : const std::vector<numeric_index_type> & dof_indices);
546 :
547 : /**
548 : * Computes \f$ \vec{u} \leftarrow \vec{u} + \vec{v} \f$,
549 : * where \p v is a DenseVector and each \p dof_indices[i] specifies
550 : * where to add value \p v(i). This method is guaranteed to be thread safe as long as library
551 : * implementations of \p NumericVector are used
552 : */
553 : void add_vector (const DenseVector<T> & v,
554 : const std::vector<numeric_index_type> & dof_indices);
555 :
556 : /**
557 : * Computes \f$ \vec{u} \leftarrow \vec{u} + A \vec{v} \f$,
558 : * i.e. adds the product of a \p SparseMatrix \p A and a \p
559 : * NumericVector \p v to \p this.
560 : */
561 : virtual void add_vector (const NumericVector<T> & v,
562 : const SparseMatrix<T> & A) = 0;
563 :
564 : /**
565 : * Computes \f$ \vec{u} \leftarrow \vec{u} + A \vec{v} \f$,
566 : * i.e. adds the product of a \p ShellMatrix \p A and a \p
567 : * NumericVector \p v to \p this.
568 : */
569 : void add_vector (const NumericVector<T> & v,
570 : const ShellMatrix<T> & A);
571 :
572 : /**
573 : * Computes \f$ \vec{u} \leftarrow \vec{u} + A^T \vec{v} \f$,
574 : * i.e. adds the product of the transpose of a \p SparseMatrix \p A
575 : * and a \p NumericVector \p v to \p this.
576 : */
577 : virtual void add_vector_transpose (const NumericVector<T> & v,
578 : const SparseMatrix<T> & A) = 0;
579 :
580 : /**
581 : * Inserts the entries of \p v in *this at the locations specified by \p v. Note that library
582 : * implementations of this method are thread safe
583 : */
584 : virtual void insert (const T * v,
585 : const std::vector<numeric_index_type> & dof_indices);
586 :
587 : /**
588 : * Inserts the entries of \p v in *this at the locations specified by \p v. This method is
589 : * guaranteed to be thread-safe as long as library implementations of \p NumericVector are used
590 : */
591 : void insert (const std::vector<T> & v,
592 : const std::vector<numeric_index_type> & dof_indices);
593 :
594 : /**
595 : * Inserts the entries of \p v in *this at the locations specified by \p v. This method is
596 : * guaranteed to be thread-safe as long as library implementations of \p NumericVector are used
597 : */
598 : void insert (const NumericVector<T> & v,
599 : const std::vector<numeric_index_type> & dof_indices);
600 :
601 : /**
602 : * Inserts the entries of \p v in *this at the locations specified by \p v. This method is
603 : * guaranteed to be thread-safe as long as library implementations of \p NumericVector are used
604 : */
605 : void insert (const DenseVector<T> & v,
606 : const std::vector<numeric_index_type> & dof_indices);
607 :
608 : /**
609 : * Inserts the entries of \p v in *this at the locations specified by \p v. This method is
610 : * guaranteed to be thread-safe as long as library implementations of \p NumericVector are used
611 : */
612 : void insert (const DenseSubVector<T> & v,
613 : const std::vector<numeric_index_type> & dof_indices);
614 :
615 : /**
616 : * Scale each element of the vector by the given \p factor.
617 : */
618 : virtual void scale (const T factor) = 0;
619 :
620 : /**
621 : * Sets \f$ u_i \leftarrow |u_i| \f$ for each entry in the vector.
622 : */
623 : virtual void abs() = 0;
624 :
625 : /**
626 : * \returns \f$ \vec{u} \cdot \vec{v} \f$, the dot product of
627 : * (*this) with the vector \p v.
628 : *
629 : * Uses the complex-conjugate of \p v in the complex-valued case.
630 : */
631 : virtual T dot(const NumericVector<T> & v) const = 0;
632 :
633 : /**
634 : * Creates a copy of the global vector in the local vector \p
635 : * v_local.
636 : */
637 : virtual void localize (std::vector<T> & v_local) const = 0;
638 :
639 : /**
640 : * Same, but fills a \p NumericVector<T> instead of a \p
641 : * std::vector.
642 : */
643 : virtual void localize (NumericVector<T> & v_local) const = 0;
644 :
645 : /**
646 : * Creates a local vector \p v_local containing only information
647 : * relevant to this processor, as defined by the \p send_list.
648 : */
649 : virtual void localize (NumericVector<T> & v_local,
650 : const std::vector<numeric_index_type> & send_list) const = 0;
651 :
652 : /**
653 : * Fill in the local std::vector "v_local" with the global indices
654 : * given in "indices".
655 : *
656 : * \note The indices can be different on every processor, and the
657 : * same index can be localized to more than one processor. The
658 : * resulting v_local can be shorter than the original, and the
659 : * entries will be in the order specified by indices.
660 : *
661 : * Example:
662 : * \verbatim
663 : * On 4 procs *this = {a, b, c, d, e, f, g, h, i} is a parallel vector.
664 : * On each proc, the indices arrays are set up as:
665 : * proc0, indices = {1,2,4,5}
666 : * proc1, indices = {2,5,6,8}
667 : * proc2, indices = {2,3,6,7}
668 : * proc3, indices = {0,1,2,3}
669 : *
670 : * After calling this version of localize, the v_local vectors are:
671 : * proc0, v_local = {b,c,e,f}
672 : * proc1, v_local = {c,f,g,i}
673 : * proc2, v_local = {c,d,g,h}
674 : * proc3, v_local = {a,b,c,d}
675 : * \endverbatim
676 : *
677 : * This function is useful in parallel I/O routines, when you have a
678 : * parallel vector of solution values which you want to write a
679 : * subset of.
680 : */
681 : virtual void localize (std::vector<T> & v_local,
682 : const std::vector<numeric_index_type> & indices) const = 0;
683 :
684 : /**
685 : * Updates a local vector with selected values from neighboring
686 : * processors, as defined by \p send_list.
687 : */
688 : virtual void localize (const numeric_index_type first_local_idx,
689 : const numeric_index_type last_local_idx,
690 : const std::vector<numeric_index_type> & send_list) = 0;
691 :
692 : /**
693 : * Creates a local copy of the global vector in \p v_local only on
694 : * processor \p proc_id. By default the data is sent to processor
695 : * 0. This method is useful for outputting data from one processor.
696 : */
697 : virtual void localize_to_one (std::vector<T> & v_local,
698 : const processor_id_type proc_id=0) const = 0;
699 :
700 : /**
701 : * \returns \p -1 when \p this is equivalent to \p other_vector
702 : * (up to the given \p threshold), or the first index where
703 : * \p abs(a[i]-b[i]) exceeds the threshold.
704 : */
705 : virtual int compare (const NumericVector<T> & other_vector,
706 : const Real threshold = TOLERANCE) const;
707 :
708 : /**
709 : * \returns \p -1 when \p this is equivalent to \p other_vector, (up
710 : * to the given local relative \p threshold), or the first index
711 : * where \p abs(a[i]-b[i])/max(a[i],b[i]) exceeds the threshold.
712 : */
713 : virtual int local_relative_compare (const NumericVector<T> & other_vector,
714 : const Real threshold = TOLERANCE) const;
715 :
716 : /**
717 : * \returns \p -1 when \p this is equivalent to \p other_vector (up
718 : * to the given global relative \p threshold), or the first index
719 : * where \p abs(a[i]-b[i])/max_j(a[j],b[j]) exceeds the threshold.
720 : */
721 : virtual int global_relative_compare (const NumericVector<T> & other_vector,
722 : const Real threshold = TOLERANCE) const;
723 :
724 : /**
725 : * Computes \f$ u_i \leftarrow u_i v_i \f$ (summation not implied)
726 : * i.e. the pointwise (component-wise) product of \p vec1 and
727 : * \p vec2, and stores the result in \p *this.
728 : */
729 : virtual void pointwise_mult (const NumericVector<T> & vec1,
730 : const NumericVector<T> & vec2) = 0;
731 :
732 : /**
733 : * Computes \f$ u_i \leftarrow \frac{v_{1,i}}{v_{2,i}} \f$ (summation not implied)
734 : * i.e. the pointwise (component-wise) division of \p vec1 and
735 : * \p vec2, and stores the result in \p *this.
736 : */
737 : virtual void pointwise_divide (const NumericVector<T> & vec1,
738 : const NumericVector<T> & vec2) = 0;
739 :
740 : /**
741 : * Prints the local contents of the vector, by default to
742 : * libMesh::out
743 : */
744 : virtual void print(std::ostream & os=libMesh::out) const;
745 :
746 : /**
747 : * Prints the global contents of the vector, by default to
748 : * libMesh::out
749 : */
750 : virtual void print_global(std::ostream & os=libMesh::out) const;
751 :
752 : /**
753 : * Same as above but allows you to use stream syntax.
754 : */
755 0 : friend std::ostream & operator << (std::ostream & os, const NumericVector<T> & v)
756 : {
757 0 : v.print_global(os);
758 0 : return os;
759 : }
760 :
761 : /**
762 : * Print the contents of the vector in Matlab's sparse matrix
763 : * format. Optionally prints the vector to the file named \p name.
764 : * If \p name is not specified it is dumped to the screen.
765 : */
766 : virtual void print_matlab(const std::string & filename = "") const;
767 :
768 : /**
769 : * Read the contents of the vector from the Matlab-script format
770 : * used by PETSc.
771 : *
772 : * If the size of the matrix in \p filename appears consistent with
773 : * the existing partitioning of \p this then the existing parallel
774 : * decomposition will be retained.
775 : * If not, then \p this will be initialized with the size from
776 : * the file, linearly partitioned onto the number of processors
777 : * available.
778 : */
779 : virtual void read_matlab(const std::string & filename);
780 :
781 : /**
782 : * Fills in \p subvector from this vector using the indices in \p rows.
783 : *
784 : * \p supplying_global_rows communicates whether
785 : * the supplied vector of rows corresponds to all the rows that
786 : * should be used in the subvector's index set, e.g. whether the
787 : * rows correspond to the global collective. If the rows supplied
788 : * are only the local indices, then \p supplying_global_rows should
789 : * be set to false
790 : *
791 : * This is currently only implemented for \p PetscVector and \p
792 : * EigenSparseVector.
793 : */
794 0 : virtual void create_subvector(NumericVector<T> & /* subvector */,
795 : const std::vector<numeric_index_type> & /* rows */,
796 : bool /* supplying_global_rows */ = true) const
797 : {
798 0 : libmesh_not_implemented();
799 : }
800 :
801 : /**
802 : * Creates a view into this vector using the indices in \p
803 : * rows. Calls to this API should always be paired with a call to \p restore_subvector,
804 : * else any changes made in the subvector will not be reflected in (\p this) original vector.
805 : *
806 : * This is currently only implemented for \p PetscVector and \p
807 : * EigenSparseVector.
808 : */
809 : virtual std::unique_ptr<NumericVector<T>>
810 0 : get_subvector(const std::vector<numeric_index_type> & /* rows */)
811 : {
812 0 : libmesh_not_implemented();
813 : }
814 :
815 : /**
816 : * Restores a view into this vector using the indices in \p
817 : * rows. The \p subvector should have been acquired from \p
818 : * get_subvector() with the same rows.
819 : *
820 : * This is currently only implemented for \p PetscVector and \p
821 : * EigenSparseVector.
822 : */
823 0 : virtual void restore_subvector(std::unique_ptr<NumericVector<T>> /* subvector */,
824 : const std::vector<numeric_index_type> & /* rows */)
825 : {
826 0 : libmesh_not_implemented();
827 : }
828 :
829 : /**
830 : * Swaps the contents of this with \p v. There should be enough
831 : * indirection in subclasses to make this an O(1) header-swap
832 : * operation.
833 : */
834 : virtual void swap (NumericVector<T> & v);
835 :
836 : /**
837 : * \returns the max number of entries (across all procs) that the
838 : * NumericVector can have.
839 : *
840 : * Although we define numeric_index_type, and it is typically
841 : * required that the NumericVector's underlying implementation's
842 : * indices have _size_ equal to sizeof(numeric_index_type), these
843 : * two types can still have different _signedness_, which affects
844 : * the maximum number of values which can be stored in the
845 : * NumericVector.
846 : */
847 : virtual std::size_t max_allowed_id() const = 0;
848 :
849 : /**
850 : * \returns whether or not this vector is able to be read for
851 : * global operations
852 : */
853 : bool readable() const;
854 :
855 : /**
856 : * \returns whether or not this vector and the vector \p v are
857 : * able to be read for global operations with one-another
858 : */
859 : bool compatible(const NumericVector<T> & v) const;
860 :
861 : protected:
862 :
863 : /**
864 : * Flag which tracks whether the vector's values are consistent on
865 : * all processors after insertion or addition of values has occurred
866 : * on some or all processors.
867 : */
868 : bool _is_closed;
869 :
870 : /**
871 : * \p true once init() has been called.
872 : */
873 : bool _is_initialized;
874 :
875 : /**
876 : * Type of vector.
877 : */
878 : ParallelType _type;
879 :
880 : /**
881 : * Mutex for performing thread-safe operations
882 : */
883 : std::mutex _numeric_vector_mutex;
884 : };
885 :
886 :
887 : /*----------------------- Inline functions ----------------------------------*/
888 :
889 :
890 :
891 : template <typename T>
892 : inline
893 3320820 : NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
894 : const ParallelType ptype) :
895 : ParallelObject(comm_in),
896 3155920 : _is_closed(false),
897 3155920 : _is_initialized(false),
898 3320148 : _type(ptype)
899 : {
900 3232200 : }
901 :
902 :
903 :
904 : template <typename T>
905 : inline
906 0 : NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
907 : const numeric_index_type /*n*/,
908 : const ParallelType ptype) :
909 : ParallelObject(comm_in),
910 0 : _is_closed(false),
911 0 : _is_initialized(false),
912 0 : _type(ptype)
913 : {
914 0 : libmesh_not_implemented(); // Abstract base class!
915 : // init(n, n, false, ptype);
916 : }
917 :
918 :
919 :
920 : template <typename T>
921 : inline
922 0 : NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
923 : const numeric_index_type /*n*/,
924 : const numeric_index_type /*n_local*/,
925 : const ParallelType ptype) :
926 : ParallelObject(comm_in),
927 0 : _is_closed(false),
928 0 : _is_initialized(false),
929 0 : _type(ptype)
930 : {
931 0 : libmesh_not_implemented(); // Abstract base class!
932 : // init(n, n_local, false, ptype);
933 : }
934 :
935 :
936 :
937 : template <typename T>
938 : inline
939 0 : NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
940 : const numeric_index_type /*n*/,
941 : const numeric_index_type /*n_local*/,
942 : const std::vector<numeric_index_type> & /*ghost*/,
943 : const ParallelType ptype) :
944 : ParallelObject(comm_in),
945 0 : _is_closed(false),
946 0 : _is_initialized(false),
947 0 : _type(ptype)
948 : {
949 0 : libmesh_not_implemented(); // Abstract base class!
950 : // init(n, n_local, ghost, false, ptype);
951 : }
952 :
953 :
954 :
955 : template <typename T>
956 : inline
957 0 : void NumericVector<T>::clear ()
958 : {
959 0 : _is_closed = false;
960 0 : _is_initialized = false;
961 0 : }
962 :
963 :
964 :
965 : template <typename T>
966 : inline
967 2054333 : void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
968 : T * values) const
969 : {
970 0 : const std::size_t num = index.size();
971 20422419 : for (std::size_t i=0; i<num; i++)
972 : {
973 18368086 : values[i] = (*this)(index[i]);
974 : }
975 2054333 : }
976 :
977 :
978 :
979 : template <typename T>
980 : inline
981 95559144 : void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
982 : std::vector<T> & values) const
983 : {
984 17213819 : const std::size_t num = index.size();
985 95559144 : values.resize(num);
986 95559144 : if (!num)
987 24428 : return;
988 :
989 95276907 : this->get(index, values.data());
990 : }
991 :
992 :
993 :
994 : template <typename T>
995 : inline
996 0 : void NumericVector<T>::add_vector(const std::vector<T> & v,
997 : const std::vector<numeric_index_type> & dof_indices)
998 : {
999 0 : libmesh_assert(v.size() == dof_indices.size());
1000 0 : if (!v.empty())
1001 0 : this->add_vector(v.data(), dof_indices);
1002 0 : }
1003 :
1004 :
1005 :
1006 : template <typename T>
1007 : inline
1008 19557562 : void NumericVector<T>::add_vector(const DenseVector<T> & v,
1009 : const std::vector<numeric_index_type> & dof_indices)
1010 : {
1011 1886303 : libmesh_assert(v.size() == dof_indices.size());
1012 21476893 : if (!v.empty())
1013 21476893 : this->add_vector(&v(0), dof_indices);
1014 19557562 : }
1015 :
1016 :
1017 :
1018 : template <typename T>
1019 : inline
1020 0 : void NumericVector<T>::insert(const std::vector<T> & v,
1021 : const std::vector<numeric_index_type> & dof_indices)
1022 : {
1023 0 : libmesh_assert(v.size() == dof_indices.size());
1024 0 : if (!v.empty())
1025 0 : this->insert(v.data(), dof_indices);
1026 0 : }
1027 :
1028 :
1029 :
1030 : template <typename T>
1031 : inline
1032 0 : void NumericVector<T>::insert(const DenseVector<T> & v,
1033 : const std::vector<numeric_index_type> & dof_indices)
1034 : {
1035 0 : libmesh_assert(v.size() == dof_indices.size());
1036 0 : if (!v.empty())
1037 0 : this->insert(&v(0), dof_indices);
1038 0 : }
1039 :
1040 :
1041 :
1042 : template <typename T>
1043 : inline
1044 1728 : void NumericVector<T>::insert(const DenseSubVector<T> & v,
1045 : const std::vector<numeric_index_type> & dof_indices)
1046 : {
1047 0 : libmesh_assert(v.size() == dof_indices.size());
1048 1728 : if (!v.empty())
1049 1728 : this->insert(&v(0), dof_indices);
1050 1728 : }
1051 :
1052 :
1053 :
1054 : // Full specialization of the print() member for complex
1055 : // variables. This must precede the non-specialized
1056 : // version, at least according to icc v7.1
1057 : template <>
1058 : inline
1059 0 : void NumericVector<Complex>::print(std::ostream & os) const
1060 : {
1061 : libmesh_assert (this->initialized());
1062 0 : os << "Size\tglobal = " << this->size()
1063 0 : << "\t\tlocal = " << this->local_size() << std::endl;
1064 :
1065 : // std::complex<>::operator<<() is defined, but use this form
1066 0 : os << "#\tReal part\t\tImaginary part" << std::endl;
1067 0 : for (auto i : index_range(*this))
1068 0 : os << i << "\t"
1069 0 : << (*this)(i).real() << "\t\t"
1070 0 : << (*this)(i).imag() << std::endl;
1071 0 : }
1072 :
1073 :
1074 :
1075 : template <typename T>
1076 : inline
1077 0 : void NumericVector<T>::print(std::ostream & os) const
1078 : {
1079 0 : libmesh_assert (this->initialized());
1080 0 : os << "Size\tglobal = " << this->size()
1081 0 : << "\t\tlocal = " << this->local_size() << std::endl;
1082 :
1083 0 : os << "#\tValue" << std::endl;
1084 0 : for (auto i : index_range(*this))
1085 0 : os << i << "\t" << (*this)(i) << std::endl;
1086 0 : }
1087 :
1088 :
1089 :
1090 : template <>
1091 : inline
1092 0 : void NumericVector<Complex>::print_global(std::ostream & os) const
1093 : {
1094 : libmesh_assert (this->initialized());
1095 :
1096 0 : std::vector<Complex> v(this->size());
1097 0 : this->localize(v);
1098 :
1099 : // Right now we only want one copy of the output
1100 0 : if (this->processor_id())
1101 : return;
1102 :
1103 0 : os << "Size\tglobal = " << this->size() << std::endl;
1104 0 : os << "#\tReal part\t\tImaginary part" << std::endl;
1105 0 : for (auto i : make_range(v.size()))
1106 0 : os << i << "\t"
1107 0 : << v[i].real() << "\t\t"
1108 : << v[i].imag() << std::endl;
1109 : }
1110 :
1111 :
1112 : template <typename T>
1113 : inline
1114 0 : void NumericVector<T>::print_global(std::ostream & os) const
1115 : {
1116 0 : libmesh_assert (this->initialized());
1117 :
1118 0 : std::vector<T> v(this->size());
1119 0 : this->localize(v);
1120 :
1121 : // Right now we only want one copy of the output
1122 0 : if (this->processor_id())
1123 0 : return;
1124 :
1125 0 : os << "Size\tglobal = " << this->size() << std::endl;
1126 0 : os << "#\tValue" << std::endl;
1127 0 : for (auto i : make_range(v.size()))
1128 0 : os << i << "\t" << v[i] << std::endl;
1129 : }
1130 :
1131 :
1132 :
1133 : template <typename T>
1134 : inline
1135 28720 : void NumericVector<T>::swap (NumericVector<T> & v)
1136 : {
1137 28720 : std::swap(_is_closed, v._is_closed);
1138 28720 : std::swap(_is_initialized, v._is_initialized);
1139 28720 : std::swap(_type, v._type);
1140 28720 : }
1141 :
1142 : template <typename T>
1143 : auto
1144 : l1_norm(const NumericVector<T> & vec)
1145 : {
1146 : return vec.l1_norm();
1147 : }
1148 :
1149 : template <typename T>
1150 : auto
1151 : l1_norm_diff(const NumericVector<T> & vec1, const NumericVector<T> & vec2)
1152 : {
1153 : return vec1.l1_norm_diff(vec2);
1154 : }
1155 :
1156 : } // namespace libMesh
1157 :
1158 :
1159 : // Workaround for weird boost/NumericVector interaction bug
1160 : #ifdef LIBMESH_DEFAULT_QUADRUPLE_PRECISION
1161 : #include <boost/mpl/bool.hpp>
1162 :
1163 : namespace boost { namespace multiprecision { namespace detail {
1164 : template <typename T, typename To>
1165 : struct is_lossy_conversion<libMesh::NumericVector<T>, To> {
1166 : typedef boost::mpl::true_ type;
1167 : static const bool value = type::value;
1168 : };
1169 : }}}
1170 : #endif
1171 :
1172 :
1173 : #endif // LIBMESH_NUMERIC_VECTOR_H
|
