laspack_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 // C++ includes
34 #include <cstdio> // for std::sprintf
35 
36 // Laspack includes
37 #include <operats.h>
38 #include <qvector.h>
39 
40 namespace libMesh
41 {
42 
43 // Forward declarations
44 template <typename T> class LaspackLinearSolver;
45 template <typename T> class SparseMatrix;
46 
55 template <typename T>
56 class LaspackVector final : public NumericVector<T>
57 {
58 public:
59 
63  explicit
64  LaspackVector (const Parallel::Communicator & comm,
65  const ParallelType = AUTOMATIC);
66 
70  explicit
71  LaspackVector (const Parallel::Communicator & comm,
72  const numeric_index_type n,
73  const ParallelType = AUTOMATIC);
74 
79  LaspackVector (const Parallel::Communicator & comm,
80  const numeric_index_type n,
81  const numeric_index_type n_local,
82  const ParallelType = AUTOMATIC);
83 
89  LaspackVector (const Parallel::Communicator & comm,
90  const numeric_index_type N,
91  const numeric_index_type n_local,
92  const std::vector<numeric_index_type> & ghost,
93  const ParallelType = AUTOMATIC);
94 
100  LaspackVector<T> & operator= (const LaspackVector<T> & v);
101 
107  LaspackVector (LaspackVector &&) = delete;
108  LaspackVector (const LaspackVector &) = delete;
109  LaspackVector & operator= (LaspackVector &&) = delete;
110  virtual ~LaspackVector ();
111 
112  virtual void close () override;
113 
114  virtual void clear () override;
115 
116  virtual void zero () override;
117 
118  virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
119 
120  virtual std::unique_ptr<NumericVector<T>> clone () const override;
121 
122  virtual void init (const numeric_index_type N,
123  const numeric_index_type n_local,
124  const bool fast=false,
125  const ParallelType ptype=AUTOMATIC) override;
126 
127  virtual void init (const numeric_index_type N,
128  const bool fast=false,
129  const ParallelType ptype=AUTOMATIC) override;
130 
131  virtual void init (const numeric_index_type N,
132  const numeric_index_type n_local,
133  const std::vector<numeric_index_type> & ghost,
134  const bool fast = false,
135  const ParallelType = AUTOMATIC) override;
136 
137  virtual void init (const NumericVector<T> & other,
138  const bool fast = false) override;
139 
140  virtual NumericVector<T> & operator= (const T s) override;
141 
142  virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
143 
144  virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
145 
146  virtual Real min () const override;
147 
148  virtual Real max () const override;
149 
150  virtual T sum () const override;
151 
152  virtual Real l1_norm () const override;
153 
154  virtual Real l2_norm () const override;
155 
156  virtual Real linfty_norm () const override;
157 
158  virtual numeric_index_type size () const override;
159 
160  virtual numeric_index_type local_size() const override;
161 
162  virtual numeric_index_type first_local_index() const override;
163 
164  virtual numeric_index_type last_local_index() const override;
165 
166  virtual T operator() (const numeric_index_type i) const override;
167 
168  virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
169 
170  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
171 
172  virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
173 
174  virtual void reciprocal() override;
175 
176  virtual void conjugate() override;
177 
178  virtual void set (const numeric_index_type i, const T value) override;
179 
180  virtual void add (const numeric_index_type i, const T value) override;
181 
182  virtual void add (const T s) override;
183 
184  virtual void add (const NumericVector<T> & v) override;
185 
186  virtual void add (const T a, const NumericVector<T> & v) override;
187 
193 
194  virtual void add_vector (const NumericVector<T> & v,
195  const SparseMatrix<T> & A) override;
196 
197  virtual void add_vector_transpose (const NumericVector<T> & v,
198  const SparseMatrix<T> & A) override;
199 
200  virtual void scale (const T factor) override;
201 
202  virtual void abs() override;
203 
204  virtual T dot(const NumericVector<T> & v) const override;
205 
206  virtual void localize (std::vector<T> & v_local) const override;
207 
208  virtual void localize (NumericVector<T> & v_local) const override;
209 
210  virtual void localize (NumericVector<T> & v_local,
211  const std::vector<numeric_index_type> & send_list) const override;
212 
213  virtual void localize (std::vector<T> & v_local,
214  const std::vector<numeric_index_type> & indices) const override;
215 
216  virtual void localize (const numeric_index_type first_local_idx,
217  const numeric_index_type last_local_idx,
218  const std::vector<numeric_index_type> & send_list) override;
219 
220  virtual void localize_to_one (std::vector<T> & v_local,
221  const processor_id_type proc_id=0) const override;
222 
223  virtual void pointwise_mult (const NumericVector<T> & vec1,
224  const NumericVector<T> & vec2) override;
225 
226  virtual void swap (NumericVector<T> & v) override;
227 
228 private:
229 
234  QVector _vec;
235 
239  friend class LaspackLinearSolver<T>;
240 };
241 
242 
243 
244 //----------------------------------------------------------
245 // LaspackVector inline methods
246 template <typename T>
247 inline
249  const ParallelType ptype)
250  : NumericVector<T>(comm, ptype)
251 {
252  this->_type = ptype;
253 }
254 
255 
256 
257 template <typename T>
258 inline
260  const numeric_index_type n,
261  const ParallelType ptype)
262  : NumericVector<T>(comm, ptype)
263 {
264  this->init(n, n, false, ptype);
265 }
266 
267 
268 
269 template <typename T>
270 inline
272  const numeric_index_type n,
273  const numeric_index_type n_local,
274  const ParallelType ptype)
275  : NumericVector<T>(comm, ptype)
276 {
277  this->init(n, n_local, false, ptype);
278 }
279 
280 
281 
282 template <typename T>
283 inline
285  const numeric_index_type N,
286  const numeric_index_type n_local,
287  const std::vector<numeric_index_type> & ghost,
288  const ParallelType ptype)
289  : NumericVector<T>(comm, ptype)
290 {
291  this->init(N, n_local, ghost, false, ptype);
292 }
293 
294 
295 
296 template <typename T>
297 inline
299 {
300  this->clear ();
301 }
302 
303 
304 
305 template <typename T>
306 inline
308  const numeric_index_type libmesh_dbg_var(n_local),
309  const bool fast,
310  const ParallelType)
311 {
312  // Laspack vectors only for serial cases,
313  // but can provide a "parallel" vector on one processor.
314  libmesh_assert_equal_to (n, n_local);
315 
316  this->_type = SERIAL;
317 
318  // Clear initialized vectors
319  if (this->initialized())
320  this->clear();
321 
322  // create a sequential vector
323 
324  static int cnt = 0;
325  char foo[80];
326  std::sprintf(foo, "Vec-%d", cnt++);
327 
328  V_Constr(&_vec, const_cast<char *>(foo), n, Normal, _LPTrue);
329 
330  this->_is_initialized = true;
331 #ifndef NDEBUG
332  this->_is_closed = true;
333 #endif
334 
335  // Optionally zero out all components
336  if (fast == false)
337  this->zero ();
338 
339  return;
340 }
341 
342 
343 
344 template <typename T>
345 inline
347  const bool fast,
348  const ParallelType ptype)
349 {
350  this->init(n,n,fast,ptype);
351 }
352 
353 
354 template <typename T>
355 inline
357  const numeric_index_type n_local,
358  const std::vector<numeric_index_type> & libmesh_dbg_var(ghost),
359  const bool fast,
360  const ParallelType ptype)
361 {
362  libmesh_assert(ghost.empty());
363  this->init(n,n_local,fast,ptype);
364 }
365 
366 
367 
368 /* Default implementation for solver packages for which ghosted
369  vectors are not yet implemented. */
370 template <class T>
372  const bool fast)
373 {
374  this->init(other.size(),other.local_size(),fast,other.type());
375 }
376 
377 
378 
379 template <typename T>
380 inline
382 {
383  libmesh_assert (this->initialized());
384 
385 #ifndef NDEBUG
386  this->_is_closed = true;
387 #endif
388 }
389 
390 
391 
392 template <typename T>
393 inline
395 {
396  if (this->initialized())
397  {
398  V_Destr (&_vec);
399  }
400 
401  this->_is_initialized = false;
402 #ifndef NDEBUG
403  this->_is_closed = false;
404 #endif
405 }
406 
407 
408 
409 template <typename T> inline
411 {
412  libmesh_assert (this->initialized());
413  libmesh_assert (this->closed());
414 
415  V_SetAllCmp (&_vec, 0.);
416 }
417 
418 
419 
420 template <typename T>
421 inline
422 std::unique_ptr<NumericVector<T>> LaspackVector<T>::zero_clone () const
423 {
424  NumericVector<T> * cloned_vector = new LaspackVector<T>(this->comm());
425 
426  cloned_vector->init(*this);
427 
428  return std::unique_ptr<NumericVector<T>>(cloned_vector);
429 }
430 
431 
432 
433 template <typename T>
434 inline
435 std::unique_ptr<NumericVector<T>> LaspackVector<T>::clone () const
436 {
437  NumericVector<T> * cloned_vector = new LaspackVector<T>(this->comm());
438 
439  cloned_vector->init(*this, true);
440 
441  *cloned_vector = *this;
442 
443  return std::unique_ptr<NumericVector<T>>(cloned_vector);
444 }
445 
446 
447 
448 template <typename T>
449 inline
451 {
452  libmesh_assert (this->initialized());
453 
454  return static_cast<numeric_index_type>(V_GetDim(const_cast<QVector*>(&_vec)));
455 }
456 
457 
458 
459 template <typename T>
460 inline
462 {
463  libmesh_assert (this->initialized());
464 
465  return this->size();
466 }
467 
468 
469 
470 template <typename T>
471 inline
473 {
474  libmesh_assert (this->initialized());
475 
476  return 0;
477 }
478 
479 
480 
481 template <typename T>
482 inline
484 {
485  libmesh_assert (this->initialized());
486 
487  return this->size();
488 }
489 
490 
491 
492 template <typename T>
493 inline
495 {
496  libmesh_assert (this->initialized());
497  libmesh_assert_less (i, this->size());
498 
499  V_SetCmp (&_vec, i+1, value);
500 
501 #ifndef NDEBUG
502  this->_is_closed = false;
503 #endif
504 }
505 
506 
507 
508 template <typename T>
509 inline
511 {
512  libmesh_assert (this->initialized());
513  libmesh_assert_less (i, this->size());
514 
515  V_AddCmp (&_vec, i+1, value);
516 
517 #ifndef NDEBUG
518  this->_is_closed = false;
519 #endif
520 }
521 
522 
523 
524 template <typename T>
525 inline
527 {
528  libmesh_assert (this->initialized());
529  libmesh_assert ( ((i >= this->first_local_index()) &&
530  (i < this->last_local_index())) );
531 
532 
533  return static_cast<T>(V_GetCmp(const_cast<QVector*>(&_vec), i+1));
534 }
535 
536 
537 
538 template <typename T>
539 inline
541 {
542  LaspackVector<T> & v = cast_ref<LaspackVector<T> &>(other);
543 
544  // This is all grossly dependent on Laspack version...
545 
546  std::swap(_vec.Name, v._vec.Name);
547  std::swap(_vec.Dim, v._vec.Dim);
548  std::swap(_vec.Instance, v._vec.Instance);
549  std::swap(_vec.LockLevel, v._vec.LockLevel);
550  std::swap(_vec.Multipl, v._vec.Multipl);
551  std::swap(_vec.OwnData, v._vec.OwnData);
552 
553  // This should still be O(1), since _vec.Cmp is just a pointer to
554  // data on the heap
555 
556  std::swap(_vec.Cmp, v._vec.Cmp);
557 }
558 
559 
560 } // namespace libMesh
561 
562 
563 #endif // #ifdef LIBMESH_HAVE_LASPACK
564 #endif // LIBMESH_LASPACK_VECTOR_H
virtual void add_vector(const NumericVector< T > &v, const SparseMatrix< T > &A) override
bool closed()
Definition: libmesh.C:265
virtual void swap(NumericVector< T > &v) override
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC) override
virtual void close() override
virtual numeric_index_type size() const override
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
virtual void abs() override
virtual void scale(const T factor) override
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
virtual T dot(const NumericVector< T > &v) const override
virtual void set(const numeric_index_type i, const T value) override
virtual numeric_index_type size() const =0
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
LaspackVector< T > & operator=(const LaspackVector< T > &v)
LaspackVector(const Parallel::Communicator &comm, const ParallelType=AUTOMATIC)
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
uint8_t processor_id_type
Definition: id_types.h:99
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
const Parallel::Communicator & comm() const
virtual T operator()(const numeric_index_type i) const override
const Number zero
Definition: libmesh.h:239
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
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:92
virtual void add(const numeric_index_type i, const T value) override
virtual Real min() const override
virtual void localize(std::vector< T > &v_local) const override
void init(triangulateio &t)
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
virtual numeric_index_type first_local_index() const override
virtual void zero() override
ParallelType type() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void conjugate() override
virtual Real linfty_norm() const override
static PetscErrorCode Mat * A
void swap(Iterator &lhs, Iterator &rhs)
virtual numeric_index_type local_size() const =0
virtual numeric_index_type local_size() const override
virtual void clear() override
static const bool value
Definition: xdr_io.C:109
bool initialized()
Definition: libmesh.C:258
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
virtual T sum() const override
virtual void reciprocal() override
virtual Real l1_norm() const override