Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 :
21 : #ifndef LIBMESH_LASPACK_VECTOR_H
22 : #define LIBMESH_LASPACK_VECTOR_H
23 :
24 :
25 :
26 : #include "libmesh/libmesh_common.h"
27 :
28 : #ifdef LIBMESH_HAVE_LASPACK
29 :
30 : // Local includes
31 : #include "libmesh/numeric_vector.h"
32 :
33 : // Laspack includes
34 : #include <operats.h>
35 : #include <qvector.h>
36 :
37 : // C++ includes
38 : #include <limits>
39 : #include <mutex>
40 :
41 : namespace libMesh
42 : {
43 :
44 : // Forward declarations
45 : template <typename T> class LaspackLinearSolver;
46 : template <typename T> class SparseMatrix;
47 :
48 : /**
49 : * This class provides a nice interface to the Laspack C-based data
50 : * structures for serial vectors. All overridden virtual functions
51 : * are documented in numeric_vector.h.
52 : *
53 : * \author Benjamin S. Kirk
54 : * \date 2002
55 : */
56 : template <typename T>
57 : class LaspackVector final : public NumericVector<T>
58 : {
59 : public:
60 :
61 : /**
62 : * Dummy-Constructor. Dimension=0
63 : */
64 : explicit
65 : LaspackVector (const Parallel::Communicator & comm,
66 : const ParallelType = AUTOMATIC);
67 :
68 : /**
69 : * Constructor. Set dimension to \p n and initialize all elements with zero.
70 : */
71 : explicit
72 : LaspackVector (const Parallel::Communicator & comm,
73 : const numeric_index_type n,
74 : const ParallelType = AUTOMATIC);
75 :
76 : /**
77 : * Constructor. Set local dimension to \p n_local, the global dimension
78 : * to \p n, and initialize all elements with zero.
79 : */
80 : LaspackVector (const Parallel::Communicator & comm,
81 : const numeric_index_type n,
82 : const numeric_index_type n_local,
83 : const ParallelType = AUTOMATIC);
84 :
85 : /**
86 : * Constructor. Set local dimension to \p n_local, the global
87 : * dimension to \p n, but additionally reserve memory for the
88 : * indices specified by the \p ghost argument.
89 : */
90 : LaspackVector (const Parallel::Communicator & comm,
91 : const numeric_index_type N,
92 : const numeric_index_type n_local,
93 : const std::vector<numeric_index_type> & ghost,
94 : const ParallelType = AUTOMATIC);
95 :
96 : /**
97 : * Copy assignment operator.
98 : * Calls Asgn_VV() to assign the contents of one vector to another.
99 : * \returns A reference to *this as the derived type.
100 : */
101 : LaspackVector<T> & operator= (const LaspackVector<T> & v);
102 :
103 : /**
104 : * This class manages a C-style struct (QVector) manually, so we
105 : * don't want to allow any automatic copy/move functions to be
106 : * generated, and we can't default the destructor.
107 : */
108 : LaspackVector (LaspackVector &&) = delete;
109 : LaspackVector (const LaspackVector &) = delete;
110 : LaspackVector & operator= (LaspackVector &&) = delete;
111 : virtual ~LaspackVector ();
112 :
113 : virtual void close () override;
114 :
115 : virtual void clear () override;
116 :
117 : virtual void zero () override;
118 :
119 : virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
120 :
121 : virtual std::unique_ptr<NumericVector<T>> clone () const override;
122 :
123 : virtual void init (const numeric_index_type N,
124 : const numeric_index_type n_local,
125 : const bool fast=false,
126 : const ParallelType ptype=AUTOMATIC) override;
127 :
128 : virtual void init (const numeric_index_type N,
129 : const bool fast=false,
130 : const ParallelType ptype=AUTOMATIC) override;
131 :
132 : virtual void init (const numeric_index_type N,
133 : const numeric_index_type n_local,
134 : const std::vector<numeric_index_type> & ghost,
135 : const bool fast = false,
136 : const ParallelType = AUTOMATIC) override;
137 :
138 : virtual void init (const NumericVector<T> & other,
139 : const bool fast = false) override;
140 :
141 : virtual NumericVector<T> & operator= (const T s) override;
142 :
143 : virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
144 :
145 : virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
146 :
147 : virtual Real min () const override;
148 :
149 : virtual Real max () const override;
150 :
151 : virtual T sum () const override;
152 :
153 : virtual Real l1_norm () const override;
154 :
155 : virtual Real l2_norm () const override;
156 :
157 : virtual Real linfty_norm () const override;
158 :
159 : virtual numeric_index_type size () const override;
160 :
161 : virtual numeric_index_type local_size() const override;
162 :
163 : virtual numeric_index_type first_local_index() const override;
164 :
165 : virtual numeric_index_type last_local_index() const override;
166 :
167 : virtual T operator() (const numeric_index_type i) const override;
168 :
169 : virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
170 :
171 : virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
172 :
173 : virtual NumericVector<T> & operator *= (const NumericVector<T> & v) override;
174 :
175 : virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
176 :
177 : virtual void reciprocal() override;
178 :
179 : virtual void conjugate() override;
180 :
181 : virtual void set (const numeric_index_type i, const T value) override;
182 :
183 : virtual void add (const numeric_index_type i, const T value) override;
184 :
185 : virtual void add (const T s) override;
186 :
187 : virtual void add (const NumericVector<T> & v) override;
188 :
189 : virtual void add (const T a, const NumericVector<T> & v) override;
190 :
191 : /**
192 : * We override one NumericVector<T>::add_vector() method but don't
193 : * want to hide the other defaults.
194 : */
195 : using NumericVector<T>::add_vector;
196 :
197 : virtual void add_vector (const NumericVector<T> & v,
198 : const SparseMatrix<T> & A) override;
199 :
200 : virtual void add_vector_transpose (const NumericVector<T> & v,
201 : const SparseMatrix<T> & A) override;
202 :
203 : virtual void scale (const T factor) override;
204 :
205 : virtual void abs() override;
206 :
207 : virtual T dot(const NumericVector<T> & v) const override;
208 :
209 : virtual void localize (std::vector<T> & v_local) const override;
210 :
211 : virtual void localize (NumericVector<T> & v_local) const override;
212 :
213 : virtual void localize (NumericVector<T> & v_local,
214 : const std::vector<numeric_index_type> & send_list) const override;
215 :
216 : virtual void localize (std::vector<T> & v_local,
217 : const std::vector<numeric_index_type> & indices) const override;
218 :
219 : virtual void localize (const numeric_index_type first_local_idx,
220 : const numeric_index_type last_local_idx,
221 : const std::vector<numeric_index_type> & send_list) override;
222 :
223 : virtual void localize_to_one (std::vector<T> & v_local,
224 : const processor_id_type proc_id=0) const override;
225 :
226 : virtual void pointwise_mult (const NumericVector<T> & vec1,
227 : const NumericVector<T> & vec2) override;
228 :
229 : virtual void pointwise_divide (const NumericVector<T> & vec1,
230 : const NumericVector<T> & vec2) override;
231 :
232 : virtual void swap (NumericVector<T> & v) override;
233 :
234 : virtual std::size_t max_allowed_id() const override;
235 :
236 : private:
237 :
238 : /**
239 : * Actual Laspack vector datatype
240 : * to hold vector entries
241 : */
242 : QVector _vec;
243 :
244 : /**
245 : * Make other Laspack datatypes friends
246 : */
247 : friend class LaspackLinearSolver<T>;
248 : };
249 :
250 :
251 :
252 : //----------------------------------------------------------
253 : // LaspackVector inline methods
254 : template <typename T>
255 : inline
256 0 : LaspackVector<T>::LaspackVector (const Parallel::Communicator & comm,
257 : const ParallelType ptype)
258 0 : : NumericVector<T>(comm, ptype)
259 : {
260 0 : this->_type = ptype;
261 0 : }
262 :
263 :
264 :
265 : template <typename T>
266 : inline
267 0 : LaspackVector<T>::LaspackVector (const Parallel::Communicator & comm,
268 : const numeric_index_type n,
269 : const ParallelType ptype)
270 0 : : NumericVector<T>(comm, ptype)
271 : {
272 0 : this->init(n, n, false, ptype);
273 0 : }
274 :
275 :
276 :
277 : template <typename T>
278 : inline
279 0 : LaspackVector<T>::LaspackVector (const Parallel::Communicator & comm,
280 : const numeric_index_type n,
281 : const numeric_index_type n_local,
282 : const ParallelType ptype)
283 0 : : NumericVector<T>(comm, ptype)
284 : {
285 0 : this->init(n, n_local, false, ptype);
286 0 : }
287 :
288 :
289 :
290 : template <typename T>
291 : inline
292 0 : LaspackVector<T>::LaspackVector (const Parallel::Communicator & comm,
293 : const numeric_index_type N,
294 : const numeric_index_type n_local,
295 : const std::vector<numeric_index_type> & ghost,
296 : const ParallelType ptype)
297 0 : : NumericVector<T>(comm, ptype)
298 : {
299 0 : this->init(N, n_local, ghost, false, ptype);
300 0 : }
301 :
302 :
303 :
304 : template <typename T>
305 : inline
306 0 : LaspackVector<T>::~LaspackVector ()
307 : {
308 0 : this->clear ();
309 0 : }
310 :
311 :
312 :
313 : template <typename T>
314 : inline
315 0 : void LaspackVector<T>::init (const numeric_index_type n,
316 : const numeric_index_type libmesh_dbg_var(n_local),
317 : const bool fast,
318 : const ParallelType)
319 : {
320 : // Laspack vectors only for serial cases,
321 : // but can provide a "parallel" vector on one processor.
322 0 : libmesh_assert_equal_to (n, n_local);
323 :
324 0 : this->_type = SERIAL;
325 :
326 : // Clear initialized vectors
327 0 : if (this->initialized())
328 0 : this->clear();
329 :
330 : // create a sequential vector
331 :
332 : static int cnt = 0;
333 0 : std::string foo{"Vec-"};
334 0 : foo += std::to_string(cnt++);
335 :
336 0 : V_Constr(&_vec, const_cast<char *>(foo.c_str()), n, Normal, _LPTrue);
337 :
338 0 : this->_is_initialized = true;
339 : #ifndef NDEBUG
340 0 : this->_is_closed = true;
341 : #endif
342 :
343 : // Optionally zero out all components
344 0 : if (fast == false)
345 0 : this->zero ();
346 :
347 0 : return;
348 : }
349 :
350 :
351 :
352 : template <typename T>
353 : inline
354 0 : void LaspackVector<T>::init (const numeric_index_type n,
355 : const bool fast,
356 : const ParallelType ptype)
357 : {
358 0 : this->init(n,n,fast,ptype);
359 0 : }
360 :
361 :
362 : template <typename T>
363 : inline
364 0 : void LaspackVector<T>::init (const numeric_index_type n,
365 : const numeric_index_type n_local,
366 : const std::vector<numeric_index_type> & libmesh_dbg_var(ghost),
367 : const bool fast,
368 : const ParallelType ptype)
369 : {
370 0 : libmesh_assert(ghost.empty());
371 0 : this->init(n,n_local,fast,ptype);
372 0 : }
373 :
374 :
375 :
376 : /* Default implementation for solver packages for which ghosted
377 : vectors are not yet implemented. */
378 : template <class T>
379 0 : void LaspackVector<T>::init (const NumericVector<T> & other,
380 : const bool fast)
381 : {
382 0 : this->init(other.size(),other.local_size(),fast,other.type());
383 0 : }
384 :
385 :
386 :
387 : template <typename T>
388 : inline
389 0 : void LaspackVector<T>::close ()
390 : {
391 0 : libmesh_assert (this->initialized());
392 :
393 : #ifndef NDEBUG
394 0 : this->_is_closed = true;
395 : #endif
396 0 : }
397 :
398 :
399 :
400 : template <typename T>
401 : inline
402 0 : void LaspackVector<T>::clear ()
403 : {
404 0 : if (this->initialized())
405 : {
406 0 : V_Destr (&_vec);
407 : }
408 :
409 0 : this->_is_initialized = false;
410 : #ifndef NDEBUG
411 0 : this->_is_closed = false;
412 : #endif
413 0 : }
414 :
415 :
416 :
417 : template <typename T> inline
418 0 : void LaspackVector<T>::zero ()
419 : {
420 0 : libmesh_assert (this->initialized());
421 0 : libmesh_assert (this->closed());
422 :
423 0 : V_SetAllCmp (&_vec, 0.);
424 0 : }
425 :
426 :
427 :
428 : template <typename T>
429 : inline
430 0 : std::unique_ptr<NumericVector<T>> LaspackVector<T>::zero_clone () const
431 : {
432 0 : std::unique_ptr<NumericVector<T>> cloned_vector =
433 : std::make_unique<LaspackVector<T>>(this->comm());
434 :
435 0 : cloned_vector->init(*this);
436 :
437 0 : return cloned_vector;
438 : }
439 :
440 :
441 :
442 : template <typename T>
443 : inline
444 0 : std::unique_ptr<NumericVector<T>> LaspackVector<T>::clone () const
445 : {
446 0 : std::unique_ptr<NumericVector<T>> cloned_vector =
447 : std::make_unique<LaspackVector<T>>(this->comm());
448 :
449 0 : cloned_vector->init(*this, true);
450 :
451 0 : *cloned_vector = *this;
452 :
453 0 : return cloned_vector;
454 : }
455 :
456 :
457 :
458 : template <typename T>
459 : inline
460 0 : numeric_index_type LaspackVector<T>::size () const
461 : {
462 0 : libmesh_assert (this->initialized());
463 :
464 64 : return static_cast<numeric_index_type>(V_GetDim(const_cast<QVector*>(&_vec)));
465 : }
466 :
467 :
468 :
469 : template <typename T>
470 : inline
471 0 : numeric_index_type LaspackVector<T>::local_size () const
472 : {
473 0 : libmesh_assert (this->initialized());
474 :
475 0 : return this->size();
476 : }
477 :
478 :
479 :
480 : template <typename T>
481 : inline
482 0 : numeric_index_type LaspackVector<T>::first_local_index () const
483 : {
484 0 : libmesh_assert (this->initialized());
485 :
486 0 : return 0;
487 : }
488 :
489 :
490 :
491 : template <typename T>
492 : inline
493 0 : numeric_index_type LaspackVector<T>::last_local_index () const
494 : {
495 0 : libmesh_assert (this->initialized());
496 :
497 0 : return this->size();
498 : }
499 :
500 :
501 :
502 : template <typename T>
503 : inline
504 120 : void LaspackVector<T>::set (const numeric_index_type i, const T value)
505 : {
506 0 : libmesh_assert (this->initialized());
507 0 : libmesh_assert_less (i, this->size());
508 :
509 120 : std::scoped_lock lock(this->_numeric_vector_mutex);
510 120 : V_SetCmp (&_vec, i+1, value);
511 :
512 : #ifndef NDEBUG
513 0 : this->_is_closed = false;
514 : #endif
515 120 : }
516 :
517 :
518 :
519 : template <typename T>
520 : inline
521 200 : void LaspackVector<T>::add (const numeric_index_type i, const T value)
522 : {
523 0 : libmesh_assert (this->initialized());
524 0 : libmesh_assert_less (i, this->size());
525 :
526 200 : std::scoped_lock lock(this->_numeric_vector_mutex);
527 200 : V_AddCmp (&_vec, i+1, value);
528 :
529 : #ifndef NDEBUG
530 0 : this->_is_closed = false;
531 : #endif
532 200 : }
533 :
534 :
535 :
536 : template <typename T>
537 : inline
538 80 : T LaspackVector<T>::operator() (const numeric_index_type i) const
539 : {
540 0 : libmesh_assert (this->initialized());
541 0 : libmesh_assert ( ((i >= this->first_local_index()) &&
542 : (i < this->last_local_index())) );
543 :
544 :
545 600 : return static_cast<T>(V_GetCmp(const_cast<QVector*>(&_vec), i+1));
546 : }
547 :
548 :
549 :
550 : template <typename T>
551 : inline
552 0 : void LaspackVector<T>::swap (NumericVector<T> & other)
553 : {
554 0 : LaspackVector<T> & v = cast_ref<LaspackVector<T> &>(other);
555 :
556 : // This is all grossly dependent on Laspack version...
557 :
558 0 : std::swap(_vec.Name, v._vec.Name);
559 0 : std::swap(_vec.Dim, v._vec.Dim);
560 0 : std::swap(_vec.Instance, v._vec.Instance);
561 0 : std::swap(_vec.LockLevel, v._vec.LockLevel);
562 0 : std::swap(_vec.Multipl, v._vec.Multipl);
563 0 : std::swap(_vec.OwnData, v._vec.OwnData);
564 :
565 : // This should still be O(1), since _vec.Cmp is just a pointer to
566 : // data on the heap
567 :
568 0 : std::swap(_vec.Cmp, v._vec.Cmp);
569 0 : }
570 :
571 :
572 :
573 : template <typename T>
574 : inline
575 0 : std::size_t LaspackVector<T>::max_allowed_id () const
576 : {
577 : // The QVector type declares a "size_t Dim;"
578 0 : return std::numeric_limits<std::size_t>::max();
579 : }
580 :
581 :
582 : } // namespace libMesh
583 :
584 :
585 : #endif // #ifdef LIBMESH_HAVE_LASPACK
586 : #endif // LIBMESH_LASPACK_VECTOR_H
|
