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_DENSE_VECTOR_H
21 : #define LIBMESH_DENSE_VECTOR_H
22 :
23 : // Local Includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/dense_vector_base.h"
26 : #include "libmesh/compare_types.h"
27 : #include "libmesh/int_range.h"
28 : #include "libmesh/tensor_tools.h"
29 :
30 : #ifdef LIBMESH_HAVE_EIGEN
31 : #include "libmesh/ignore_warnings.h"
32 : #include <Eigen/Core>
33 : #include "libmesh/restore_warnings.h"
34 : #endif
35 :
36 : #ifdef LIBMESH_HAVE_METAPHYSICL
37 : #include "metaphysicl/raw_type.h"
38 : #endif
39 :
40 : // C++ includes
41 : #include <vector>
42 : #include <initializer_list>
43 :
44 : namespace libMesh
45 : {
46 :
47 : /**
48 : * Defines a dense vector for use in Finite Element-type computations.
49 : * This class is to basically compliment the \p DenseMatrix class. It
50 : * has additional capabilities over the \p std::vector that make it
51 : * useful for finite elements, particularly for systems of equations.
52 : * All overridden virtual functions are documented in dense_vector_base.h.
53 : *
54 : * \author Benjamin S. Kirk
55 : * \date 2003
56 : */
57 : template<typename T>
58 6502826 : class DenseVector : public DenseVectorBase<T>
59 : {
60 : public:
61 :
62 : /**
63 : * Constructor. Creates a dense vector of dimension \p n.
64 : */
65 : explicit
66 : DenseVector(const unsigned int n=0);
67 :
68 : /**
69 : * Constructor. Creates a dense vector of dimension \p n where all
70 : * entries have value \p val.
71 : */
72 : explicit
73 : DenseVector(const unsigned int n,
74 : const T & val);
75 :
76 : /**
77 : * Copy-constructor.
78 : */
79 : template <typename T2>
80 : DenseVector (const DenseVector<T2> & other_vector);
81 :
82 : /**
83 : * Copy-constructor, from a \p std::vector.
84 : */
85 : template <typename T2>
86 : DenseVector (const std::vector<T2> & other_vector);
87 :
88 : /**
89 : * Initializer list constructor.
90 : */
91 : template <typename T2>
92 : DenseVector (std::initializer_list<T2> init_list);
93 :
94 : /**
95 : * The 5 special functions can be defaulted for this class, as it
96 : * does not manage any memory itself.
97 : */
98 176256 : DenseVector (DenseVector &&) = default;
99 25896342 : DenseVector (const DenseVector &) = default;
100 28517995 : DenseVector & operator= (const DenseVector &) = default;
101 166464 : DenseVector & operator= (DenseVector &&) = default;
102 46124484 : virtual ~DenseVector() = default;
103 :
104 5425412 : virtual unsigned int size() const override final
105 : {
106 15550456 : return cast_int<unsigned int>(_val.size());
107 : }
108 :
109 19225 : virtual bool empty() const override final
110 19225 : { return _val.empty(); }
111 :
112 : virtual void zero() override final;
113 :
114 : /**
115 : * \returns Entry \p i of the vector as a const reference.
116 : */
117 : const T & operator() (const unsigned int i) const;
118 :
119 : /**
120 : * \returns Entry \p i of the vector as a writable reference.
121 : */
122 : T & operator() (const unsigned int i);
123 :
124 : /**
125 : * \returns Entry \p i of the vector as a const reference.
126 : */
127 0 : const T & operator[] (const unsigned int i) const { return (*this)(i); }
128 :
129 : /**
130 : * \returns Entry \p i of the vector as a writable reference.
131 : */
132 0 : T & operator[] (const unsigned int i) { return (*this)(i); }
133 :
134 0 : virtual T el(const unsigned int i) const override final
135 0 : { return (*this)(i); }
136 :
137 0 : virtual T & el(const unsigned int i) override final
138 0 : { return (*this)(i); }
139 :
140 : /**
141 : * Assignment operator.
142 : *
143 : * \returns A reference to *this.
144 : */
145 : template <typename T2>
146 : DenseVector<T> & operator = (const DenseVector<T2> & other_vector);
147 :
148 : /**
149 : * STL-like swap method
150 : */
151 : void swap(DenseVector<T> & other_vector);
152 :
153 : /**
154 : * Resize the vector. Sets all elements to 0.
155 : */
156 : void resize (const unsigned int n);
157 :
158 : /**
159 : * Append additional entries to (resizing, but unchanging) the
160 : * vector.
161 : */
162 : template <typename T2>
163 : void append (const DenseVector<T2> & other_vector);
164 :
165 : /**
166 : * Multiplies every element in the vector by \p factor.
167 : */
168 : void scale (const T factor);
169 :
170 : /**
171 : * Multiplies every element in the vector by \p factor.
172 : *
173 : * \returns A reference to *this.
174 : */
175 : DenseVector<T> & operator*= (const T factor);
176 :
177 : /**
178 : * Adds \p factor times \p vec to this vector.
179 : * This should only work if T += T2 * T3 is valid C++ and
180 : * if T2 is scalar. Return type is void
181 : *
182 : * \returns A reference to *this.
183 : */
184 : template <typename T2, typename T3>
185 : typename boostcopy::enable_if_c<
186 : ScalarTraits<T2>::value, void >::type
187 : add (const T2 factor,
188 : const DenseVector<T3> & vec);
189 :
190 : /**
191 : * \returns The dot product of *this with \p vec.
192 : *
193 : * In the complex-valued case, uses the complex conjugate of \p vec.
194 : */
195 : template <typename T2>
196 : typename CompareTypes<T, T2>::supertype dot (const DenseVector<T2> & vec) const;
197 :
198 : /**
199 : * \returns The dot product of *this with \p vec.
200 : *
201 : * In the complex-valued case, does not use the complex conjugate of
202 : * \p vec.
203 : */
204 : template <typename T2>
205 : typename CompareTypes<T, T2>::supertype indefinite_dot (const DenseVector<T2> & vec) const;
206 :
207 : /**
208 : * \returns \p true if \p vec is exactly equal to this vector, false otherwise.
209 : */
210 : template <typename T2>
211 : bool operator== (const DenseVector<T2> & vec) const;
212 :
213 : /**
214 : * \returns \p true if \p vec is not exactly equal to this vector, false otherwise.
215 : */
216 : template <typename T2>
217 : bool operator!= (const DenseVector<T2> & vec) const;
218 :
219 : /**
220 : * Adds \p vec to this vector.
221 : *
222 : * \returns A reference to *this.
223 : */
224 : template <typename T2>
225 : DenseVector<T> & operator+= (const DenseVector<T2> & vec);
226 :
227 : /**
228 : * Subtracts \p vec from this vector.
229 : *
230 : * \returns A reference to *this.
231 : */
232 : template <typename T2>
233 : DenseVector<T> & operator-= (const DenseVector<T2> & vec);
234 :
235 : /**
236 : * \returns The minimum entry of the vector, or the minimum real
237 : * part in the case of complex numbers.
238 : */
239 : Real min () const;
240 :
241 : /**
242 : * \returns The maximum entry of the vector, or the maximum real
243 : * part in the case of complex numbers.
244 : */
245 : Real max () const;
246 :
247 : /**
248 : * \returns The \f$l_1\f$-norm of the vector, i.e. the sum of the
249 : * absolute values of the entries.
250 : */
251 : Real l1_norm () const;
252 :
253 : /**
254 : * \returns The \f$l_2\f$-norm of the vector, i.e. the square root
255 : * of the sum of the squares of the entries.
256 : */
257 : Real l2_norm () const;
258 :
259 : /**
260 : * \returns The \f$l_\infty\f$-norm of the vector, i.e. the maximum
261 : * absolute value of the entries.
262 : */
263 : Real linfty_norm () const;
264 :
265 : /**
266 : * Puts the principal subvector of size \p sub_n (i.e. first sub_n
267 : * entries) into \p dest.
268 : */
269 : void get_principal_subvector (unsigned int sub_n, DenseVector<T> & dest) const;
270 :
271 : /**
272 : * \returns A reference to the underlying data storage vector.
273 : *
274 : * This should be used with caution (i.e. one should not change the
275 : * size of the vector, etc.) but is useful for interoperating with
276 : * low level BLAS routines which expect a simple array.
277 : */
278 75988403 : std::vector<T> & get_values() { return _val; }
279 :
280 : /**
281 : * \returns A constant reference to the underlying data storage vector.
282 : */
283 1164024 : const std::vector<T> & get_values() const { return _val; }
284 :
285 : /**
286 : * \returns An iterator pointing to the beginning of the vector
287 : */
288 0 : typename std::vector<T>::const_iterator begin() const { return _val.begin(); }
289 0 : typename std::vector<T>::iterator begin() { return _val.begin(); }
290 :
291 : /**
292 : * \returns An iterator pointing to the end of the vector
293 : */
294 0 : typename std::vector<T>::const_iterator end() const { return _val.end(); }
295 0 : typename std::vector<T>::iterator end() { return _val.end(); }
296 :
297 : private:
298 :
299 : /**
300 : * The actual data values, stored as a 1D array.
301 : */
302 : std::vector<T> _val;
303 : };
304 :
305 :
306 :
307 : // ------------------------------------------------------------
308 : // DenseVector member functions
309 : template<typename T>
310 : inline
311 40119737 : DenseVector<T>::DenseVector(const unsigned int n) :
312 39214343 : _val (n, T{})
313 : {
314 35669042 : }
315 :
316 :
317 : template<typename T>
318 : inline
319 0 : DenseVector<T>::DenseVector(const unsigned int n,
320 : const T & val) :
321 0 : _val (n, val)
322 : {
323 0 : }
324 :
325 :
326 :
327 : template<typename T>
328 : template<typename T2>
329 : inline
330 : DenseVector<T>::DenseVector (const DenseVector<T2> & other_vector) :
331 : DenseVectorBase<T>()
332 : {
333 : const std::vector<T2> & other_vals = other_vector.get_values();
334 :
335 : _val.clear();
336 :
337 : const int N = cast_int<int>(other_vals.size());
338 : _val.reserve(N);
339 :
340 : for (int i=0; i<N; i++)
341 : _val.push_back(other_vals[i]);
342 : }
343 :
344 :
345 :
346 : template<typename T>
347 : template<typename T2>
348 : inline
349 2409837 : DenseVector<T>::DenseVector (const std::vector<T2> & other_vector) :
350 2409837 : _val(other_vector)
351 : {
352 2134303 : }
353 :
354 :
355 : template<typename T>
356 : template <typename T2>
357 : inline
358 : DenseVector<T>::DenseVector (std::initializer_list<T2> init_list) :
359 : _val(init_list.begin(), init_list.end())
360 : {
361 : }
362 :
363 :
364 :
365 : template<typename T>
366 : template<typename T2>
367 : inline
368 : DenseVector<T> & DenseVector<T>::operator = (const DenseVector<T2> & other_vector)
369 : {
370 : const std::vector<T2> & other_vals = other_vector.get_values();
371 :
372 : _val.clear();
373 :
374 : const int N = cast_int<int>(other_vals.size());
375 : _val.reserve(N);
376 :
377 : for (int i=0; i<N; i++)
378 : _val.push_back(other_vals[i]);
379 :
380 : return *this;
381 : }
382 :
383 :
384 :
385 : template<typename T>
386 : inline
387 3906688 : void DenseVector<T>::swap(DenseVector<T> & other_vector)
388 : {
389 3906688 : _val.swap(other_vector._val);
390 3906688 : }
391 :
392 :
393 :
394 : template<typename T>
395 : inline
396 144375060 : void DenseVector<T>::resize(const unsigned int n)
397 : {
398 157777611 : _val.resize(n);
399 :
400 11581858 : zero();
401 144375060 : }
402 :
403 :
404 :
405 : template<typename T>
406 : template<typename T2>
407 : inline
408 0 : void DenseVector<T>::append (const DenseVector<T2> & other_vector)
409 : {
410 0 : const std::vector<T2> & other_vals = other_vector.get_values();
411 :
412 0 : _val.reserve(this->size() + other_vals.size());
413 0 : _val.insert(_val.end(), other_vals.begin(), other_vals.end());
414 0 : }
415 :
416 :
417 :
418 : template<typename T>
419 : inline
420 1531137 : void DenseVector<T>::zero()
421 : {
422 1531137 : std::fill (_val.begin(),
423 : _val.end(),
424 : T{});
425 1531137 : }
426 :
427 :
428 :
429 : template<typename T>
430 : inline
431 3222565 : const T & DenseVector<T>::operator () (const unsigned int i) const
432 : {
433 3222565 : libmesh_assert_less (i, _val.size());
434 :
435 513332817 : return _val[i];
436 : }
437 :
438 :
439 :
440 : template<typename T>
441 : inline
442 2700 : T & DenseVector<T>::operator () (const unsigned int i)
443 : {
444 2700 : libmesh_assert_less (i, _val.size());
445 :
446 6022638642 : return _val[i];
447 : }
448 :
449 :
450 :
451 : template<typename T>
452 : inline
453 40582384 : void DenseVector<T>::scale (const T factor)
454 : {
455 5986094 : const int N = cast_int<int>(_val.size());
456 442730123 : for (int i=0; i<N; i++)
457 432945952 : _val[i] *= factor;
458 40582384 : }
459 :
460 :
461 :
462 : template<typename T>
463 : inline
464 3997375 : DenseVector<T> & DenseVector<T>::operator*= (const T factor)
465 : {
466 40582384 : this->scale(factor);
467 3997375 : return *this;
468 : }
469 :
470 :
471 :
472 : template<typename T>
473 : template<typename T2, typename T3>
474 : inline
475 : typename boostcopy::enable_if_c<
476 : ScalarTraits<T2>::value, void >::type
477 24084850 : DenseVector<T>::add (const T2 factor,
478 : const DenseVector<T3> & vec)
479 : {
480 2355460 : libmesh_assert_equal_to (this->size(), vec.size());
481 :
482 2739688 : const int N = cast_int<int>(_val.size());
483 237529944 : for (int i=0; i<N; i++)
484 247786298 : (*this)(i) += static_cast<T>(factor)*vec(i);
485 24084850 : }
486 :
487 :
488 :
489 : template<typename T>
490 : template<typename T2>
491 : inline
492 13517648 : typename CompareTypes<T, T2>::supertype DenseVector<T>::dot (const DenseVector<T2> & vec) const
493 : {
494 14742872 : if (!_val.size())
495 61600 : return 0.;
496 :
497 1164024 : libmesh_assert_equal_to (this->size(), vec.size());
498 :
499 : #ifdef LIBMESH_HAVE_EIGEN
500 : // We reverse the order of the arguments to dot() here since
501 : // the convention in Eigen is to take the complex conjugate of the
502 : // *first* argument, while ours is to take the complex conjugate of
503 : // the second.
504 11679624 : return Eigen::Map<const typename Eigen::Matrix<T2, Eigen::Dynamic, 1>>(vec.get_values().data(), vec.size())
505 12843648 : .dot(Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()));
506 : #else
507 : typename CompareTypes<T, T2>::supertype val = 0.;
508 :
509 : const int N = cast_int<int>(_val.size());
510 : // The following pragma tells clang's vectorizer that it is safe to
511 : // reorder floating point operations for this loop.
512 : #ifdef __clang__
513 : #pragma clang loop vectorize(enable)
514 : #endif
515 : for (int i=0; i<N; i++)
516 : val += (*this)(i)*libmesh_conj(vec(i));
517 :
518 : return val;
519 : #endif
520 : }
521 :
522 : template<typename T>
523 : template<typename T2>
524 : inline
525 : typename CompareTypes<T, T2>::supertype DenseVector<T>::indefinite_dot (const DenseVector<T2> & vec) const
526 : {
527 : libmesh_assert_equal_to (this->size(), vec.size());
528 :
529 : typename CompareTypes<T, T2>::supertype val = 0.;
530 :
531 : const int N = cast_int<int>(_val.size());
532 : for (int i=0; i<N; i++)
533 : val += (*this)(i)*(vec(i));
534 :
535 : return val;
536 : }
537 :
538 : template<typename T>
539 : template<typename T2>
540 : inline
541 : bool DenseVector<T>::operator== (const DenseVector<T2> & vec) const
542 : {
543 : libmesh_assert_equal_to (this->size(), vec.size());
544 :
545 : const int N = cast_int<int>(_val.size());
546 : for (int i=0; i<N; i++)
547 : if ((*this)(i) != vec(i))
548 : return false;
549 :
550 : return true;
551 : }
552 :
553 :
554 :
555 : template<typename T>
556 : template<typename T2>
557 : inline
558 : bool DenseVector<T>::operator!= (const DenseVector<T2> & vec) const
559 : {
560 : libmesh_assert_equal_to (this->size(), vec.size());
561 :
562 : const int N = cast_int<int>(_val.size());
563 : for (int i=0; i<N; i++)
564 : if ((*this)(i) != vec(i))
565 : return true;
566 :
567 : return false;
568 : }
569 :
570 :
571 :
572 : template<typename T>
573 : template<typename T2>
574 : inline
575 933728 : DenseVector<T> & DenseVector<T>::operator+= (const DenseVector<T2> & vec)
576 : {
577 92736 : libmesh_assert_equal_to (this->size(), vec.size());
578 :
579 139104 : const int N = cast_int<int>(_val.size());
580 43113344 : for (int i=0; i<N; i++)
581 47488616 : (*this)(i) += vec(i);
582 :
583 933728 : return *this;
584 : }
585 :
586 :
587 :
588 : template<typename T>
589 : template<typename T2>
590 : inline
591 20540847 : DenseVector<T> & DenseVector<T>::operator-= (const DenseVector<T2> & vec)
592 : {
593 2023299 : libmesh_assert_equal_to (this->size(), vec.size());
594 :
595 4048379 : const int N = cast_int<int>(_val.size());
596 229289807 : for (int i=0; i<N; i++)
597 234344342 : (*this)(i) -= vec(i);
598 :
599 20540847 : return *this;
600 : }
601 :
602 :
603 :
604 : template<typename T>
605 : inline
606 0 : Real DenseVector<T>::min () const
607 : {
608 0 : libmesh_assert (this->size());
609 0 : Real my_min = libmesh_real((*this)(0));
610 :
611 0 : const int N = cast_int<int>(_val.size());
612 0 : for (int i=1; i!=N; i++)
613 : {
614 0 : Real current = libmesh_real((*this)(i));
615 0 : my_min = (my_min < current? my_min : current);
616 : }
617 0 : return my_min;
618 : }
619 :
620 :
621 :
622 : template<typename T>
623 : inline
624 0 : Real DenseVector<T>::max () const
625 : {
626 0 : libmesh_assert (this->size());
627 0 : Real my_max = libmesh_real((*this)(0));
628 :
629 0 : const int N = cast_int<int>(_val.size());
630 0 : for (int i=1; i!=N; i++)
631 : {
632 0 : Real current = libmesh_real((*this)(i));
633 0 : my_max = (my_max > current? my_max : current);
634 : }
635 0 : return my_max;
636 : }
637 :
638 :
639 :
640 : template<typename T>
641 : inline
642 0 : Real DenseVector<T>::l1_norm () const
643 : {
644 0 : if (!_val.size())
645 0 : return 0.;
646 :
647 : #ifdef LIBMESH_HAVE_EIGEN
648 0 : return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<1>();
649 : #else
650 : Real my_norm = 0.;
651 : const int N = cast_int<int>(_val.size());
652 : for (int i=0; i!=N; i++)
653 : my_norm += std::abs((*this)(i));
654 :
655 : return my_norm;
656 : #endif
657 : }
658 :
659 :
660 :
661 : template<typename T>
662 : inline
663 0 : Real DenseVector<T>::l2_norm () const
664 : {
665 0 : if (!_val.size())
666 0 : return 0.;
667 :
668 : #ifdef LIBMESH_HAVE_EIGEN
669 0 : return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).norm();
670 : #else
671 : Real my_norm = 0.;
672 : const int N = cast_int<int>(_val.size());
673 : // The following pragma tells clang's vectorizer that it is safe to
674 : // reorder floating point operations for this loop.
675 : #ifdef __clang__
676 : #pragma clang loop vectorize(enable)
677 : #endif
678 : for (int i=0; i!=N; i++)
679 : my_norm += TensorTools::norm_sq((*this)(i));
680 :
681 : return sqrt(my_norm);
682 : #endif
683 : }
684 :
685 :
686 :
687 : template<typename T>
688 : inline
689 0 : Real DenseVector<T>::linfty_norm () const
690 : {
691 0 : if (!_val.size())
692 0 : return 0.;
693 :
694 : #ifdef LIBMESH_HAVE_EIGEN
695 0 : return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<Eigen::Infinity>();
696 : #else
697 : Real my_norm = TensorTools::norm_sq((*this)(0));
698 :
699 : const int N = cast_int<int>(_val.size());
700 : for (int i=1; i!=N; i++)
701 : {
702 : Real current = TensorTools::norm_sq((*this)(i));
703 : my_norm = (my_norm > current? my_norm : current);
704 : }
705 : return sqrt(my_norm);
706 : #endif
707 : }
708 :
709 :
710 :
711 : template<typename T>
712 : inline
713 11312742 : void DenseVector<T>::get_principal_subvector (unsigned int sub_n,
714 : DenseVector<T> & dest) const
715 : {
716 1008892 : libmesh_assert_less_equal ( sub_n, this->size() );
717 :
718 10303850 : dest.resize(sub_n);
719 1008892 : const int N = cast_int<int>(sub_n);
720 113671452 : for (int i=0; i<N; i++)
721 120645390 : dest(i) = _val[i];
722 11312742 : }
723 :
724 : } // namespace libMesh
725 :
726 : #ifdef LIBMESH_HAVE_METAPHYSICL
727 : namespace MetaPhysicL
728 : {
729 : template <typename T>
730 : struct RawType<libMesh::DenseVector<T>>
731 : {
732 : typedef libMesh::DenseVector<typename RawType<T>::value_type> value_type;
733 :
734 : static value_type value (const libMesh::DenseVector<T> & in)
735 : {
736 : const auto s = in.size();
737 : value_type ret(s);
738 : for (unsigned int i = 0; i < s; ++i)
739 : ret(i) = raw_value(in(i));
740 :
741 : return ret;
742 : }
743 : };
744 : }
745 : #endif
746 :
747 : #endif // LIBMESH_DENSE_VECTOR_H
|
