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