libMesh
laspack_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 
56 template <typename T>
57 class LaspackVector final : public NumericVector<T>
58 {
59 public:
60 
64  explicit
65  LaspackVector (const Parallel::Communicator & comm,
66  const ParallelType = AUTOMATIC);
67 
71  explicit
72  LaspackVector (const Parallel::Communicator & comm,
73  const numeric_index_type n,
74  const ParallelType = AUTOMATIC);
75 
80  LaspackVector (const Parallel::Communicator & comm,
81  const numeric_index_type n,
82  const numeric_index_type n_local,
83  const ParallelType = AUTOMATIC);
84 
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 
101  LaspackVector<T> & operator= (const LaspackVector<T> & v);
102 
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 
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 
242  QVector _vec;
243 
247  friend class LaspackLinearSolver<T>;
248 };
249 
250 
251 
252 //----------------------------------------------------------
253 // LaspackVector inline methods
254 template <typename T>
255 inline
257  const ParallelType ptype)
258  : NumericVector<T>(comm, ptype)
259 {
260  this->_type = ptype;
261 }
262 
263 
264 
265 template <typename T>
266 inline
268  const numeric_index_type n,
269  const ParallelType ptype)
270  : NumericVector<T>(comm, ptype)
271 {
272  this->init(n, n, false, ptype);
273 }
274 
275 
276 
277 template <typename T>
278 inline
280  const numeric_index_type n,
281  const numeric_index_type n_local,
282  const ParallelType ptype)
283  : NumericVector<T>(comm, ptype)
284 {
285  this->init(n, n_local, false, ptype);
286 }
287 
288 
289 
290 template <typename T>
291 inline
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  : NumericVector<T>(comm, ptype)
298 {
299  this->init(N, n_local, ghost, false, ptype);
300 }
301 
302 
303 
304 template <typename T>
305 inline
307 {
308  this->clear ();
309 }
310 
311 
312 
313 template <typename T>
314 inline
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  libmesh_assert_equal_to (n, n_local);
323 
324  this->_type = SERIAL;
325 
326  // Clear initialized vectors
327  if (this->initialized())
328  this->clear();
329 
330  // create a sequential vector
331 
332  static int cnt = 0;
333  std::string foo{"Vec-"};
334  foo += std::to_string(cnt++);
335 
336  V_Constr(&_vec, const_cast<char *>(foo.c_str()), n, Normal, _LPTrue);
337 
338  this->_is_initialized = true;
339 #ifndef NDEBUG
340  this->_is_closed = true;
341 #endif
342 
343  // Optionally zero out all components
344  if (fast == false)
345  this->zero ();
346 
347  return;
348 }
349 
350 
351 
352 template <typename T>
353 inline
355  const bool fast,
356  const ParallelType ptype)
357 {
358  this->init(n,n,fast,ptype);
359 }
360 
361 
362 template <typename T>
363 inline
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  libmesh_assert(ghost.empty());
371  this->init(n,n_local,fast,ptype);
372 }
373 
374 
375 
376 /* Default implementation for solver packages for which ghosted
377  vectors are not yet implemented. */
378 template <class T>
380  const bool fast)
381 {
382  this->init(other.size(),other.local_size(),fast,other.type());
383 }
384 
385 
386 
387 template <typename T>
388 inline
390 {
391  libmesh_assert (this->initialized());
392 
393 #ifndef NDEBUG
394  this->_is_closed = true;
395 #endif
396 }
397 
398 
399 
400 template <typename T>
401 inline
403 {
404  if (this->initialized())
405  {
406  V_Destr (&_vec);
407  }
408 
409  this->_is_initialized = false;
410 #ifndef NDEBUG
411  this->_is_closed = false;
412 #endif
413 }
414 
415 
416 
417 template <typename T> inline
419 {
420  libmesh_assert (this->initialized());
421  libmesh_assert (this->closed());
422 
423  V_SetAllCmp (&_vec, 0.);
424 }
425 
426 
427 
428 template <typename T>
429 inline
430 std::unique_ptr<NumericVector<T>> LaspackVector<T>::zero_clone () const
431 {
432  NumericVector<T> * cloned_vector = new LaspackVector<T>(this->comm());
433 
434  cloned_vector->init(*this);
435 
436  return std::unique_ptr<NumericVector<T>>(cloned_vector);
437 }
438 
439 
440 
441 template <typename T>
442 inline
443 std::unique_ptr<NumericVector<T>> LaspackVector<T>::clone () const
444 {
445  NumericVector<T> * cloned_vector = new LaspackVector<T>(this->comm());
446 
447  cloned_vector->init(*this, true);
448 
449  *cloned_vector = *this;
450 
451  return std::unique_ptr<NumericVector<T>>(cloned_vector);
452 }
453 
454 
455 
456 template <typename T>
457 inline
459 {
460  libmesh_assert (this->initialized());
461 
462  return static_cast<numeric_index_type>(V_GetDim(const_cast<QVector*>(&_vec)));
463 }
464 
465 
466 
467 template <typename T>
468 inline
470 {
471  libmesh_assert (this->initialized());
472 
473  return this->size();
474 }
475 
476 
477 
478 template <typename T>
479 inline
481 {
482  libmesh_assert (this->initialized());
483 
484  return 0;
485 }
486 
487 
488 
489 template <typename T>
490 inline
492 {
493  libmesh_assert (this->initialized());
494 
495  return this->size();
496 }
497 
498 
499 
500 template <typename T>
501 inline
503 {
504  libmesh_assert (this->initialized());
505  libmesh_assert_less (i, this->size());
506 
507  std::scoped_lock lock(this->_numeric_vector_mutex);
508  V_SetCmp (&_vec, i+1, value);
509 
510 #ifndef NDEBUG
511  this->_is_closed = false;
512 #endif
513 }
514 
515 
516 
517 template <typename T>
518 inline
520 {
521  libmesh_assert (this->initialized());
522  libmesh_assert_less (i, this->size());
523 
524  std::scoped_lock lock(this->_numeric_vector_mutex);
525  V_AddCmp (&_vec, i+1, value);
526 
527 #ifndef NDEBUG
528  this->_is_closed = false;
529 #endif
530 }
531 
532 
533 
534 template <typename T>
535 inline
537 {
538  libmesh_assert (this->initialized());
539  libmesh_assert ( ((i >= this->first_local_index()) &&
540  (i < this->last_local_index())) );
541 
542 
543  return static_cast<T>(V_GetCmp(const_cast<QVector*>(&_vec), i+1));
544 }
545 
546 
547 
548 template <typename T>
549 inline
551 {
552  LaspackVector<T> & v = cast_ref<LaspackVector<T> &>(other);
553 
554  // This is all grossly dependent on Laspack version...
555 
556  std::swap(_vec.Name, v._vec.Name);
557  std::swap(_vec.Dim, v._vec.Dim);
558  std::swap(_vec.Instance, v._vec.Instance);
559  std::swap(_vec.LockLevel, v._vec.LockLevel);
560  std::swap(_vec.Multipl, v._vec.Multipl);
561  std::swap(_vec.OwnData, v._vec.OwnData);
562 
563  // This should still be O(1), since _vec.Cmp is just a pointer to
564  // data on the heap
565 
566  std::swap(_vec.Cmp, v._vec.Cmp);
567 }
568 
569 
570 
571 template <typename T>
572 inline
574 {
575  // The QVector type declares a "size_t Dim;"
576  return std::numeric_limits<std::size_t>::max();
577 }
578 
579 
580 } // namespace libMesh
581 
582 
583 #endif // #ifdef LIBMESH_HAVE_LASPACK
584 #endif // LIBMESH_LASPACK_VECTOR_H
virtual void add_vector(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:273
virtual void swap(NumericVector< T > &v) override
Swaps the contents of this with v.
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC) override
Change the dimension of the vector to n.
virtual void close() override
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual numeric_index_type size() const override
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
Subtracts v from *this, .
virtual void abs() override
Sets for each entry in the vector.
virtual void scale(const T factor) override
Scale each element of the vector by the given factor.
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
virtual T dot(const NumericVector< T > &v) const override
virtual void set(const numeric_index_type i, const T value) override
Sets v(i) = value.
virtual numeric_index_type size() const =0
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
LaspackVector< T > & operator=(const LaspackVector< T > &v)
Copy assignment operator.
LaspackVector(const Parallel::Communicator &comm, const ParallelType=AUTOMATIC)
Dummy-Constructor.
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:43
const Parallel::Communicator & comm() const
virtual T operator()(const numeric_index_type i) const override
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:280
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
Adds v to *this, .
virtual NumericVector< T > & operator*=(const NumericVector< T > &v) override
Computes the component-wise multiplication of this vector&#39;s entries by another&#39;s, ...
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
Change the dimension of the vector to n.
virtual std::size_t max_allowed_id() const override
This class provides an interface to Laspack iterative solvers that is compatible with the libMesh Lin...
uint8_t processor_id_type
virtual std::unique_ptr< NumericVector< T > > zero_clone() const override
virtual std::unique_ptr< NumericVector< T > > clone() const override
virtual Real max() const override
virtual Real l2_norm() const override
dof_id_type numeric_index_type
Definition: id_types.h:99
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:247
virtual void add(const numeric_index_type i, const T value) override
Adds value to the vector entry specified by i.
virtual Real min() const override
virtual void localize(std::vector< T > &v_local) const override
Creates a copy of the global vector in the local vector v_local.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libmesh_assert(ctx)
virtual numeric_index_type last_local_index() const override
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const override
Creates a local copy of the global vector in v_local only on processor proc_id.
ParallelType _type
Type of vector.
virtual numeric_index_type first_local_index() const override
This class provides a nice interface to the Laspack C-based data structures for serial vectors...
virtual void zero() override
Set all entries to zero.
ParallelType type() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void conjugate() override
Negates the imaginary component of each entry in the vector.
virtual Real linfty_norm() const override
virtual numeric_index_type local_size() const =0
virtual numeric_index_type local_size() const override
virtual void clear() override
Restores the NumericVector<T> to a pristine state.
static const bool value
Definition: xdr_io.C:54
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:266
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
Computes the component-wise division of this vector&#39;s entries by another&#39;s, .
virtual T sum() const override
QVector _vec
Actual Laspack vector datatype to hold vector entries.
virtual void reciprocal() override
Computes the component-wise reciprocal, .
virtual Real l1_norm() const override
virtual void pointwise_divide(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
ParallelType
Defines an enum for parallel data structure types.