libMesh
laspack_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 libmesh_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 
99  ~LaspackVector ();
100 
101  virtual void close () libmesh_override;
102 
103  virtual void clear () libmesh_override;
104 
105  virtual void zero () libmesh_override;
106 
107  virtual UniquePtr<NumericVector<T>> zero_clone () const libmesh_override;
108 
109  virtual UniquePtr<NumericVector<T>> clone () const libmesh_override;
110 
111  virtual void init (const numeric_index_type N,
112  const numeric_index_type n_local,
113  const bool fast=false,
114  const ParallelType ptype=AUTOMATIC) libmesh_override;
115 
116  virtual void init (const numeric_index_type N,
117  const bool fast=false,
118  const ParallelType ptype=AUTOMATIC) libmesh_override;
119 
120  virtual void init (const numeric_index_type N,
121  const numeric_index_type n_local,
122  const std::vector<numeric_index_type> & ghost,
123  const bool fast = false,
124  const ParallelType = AUTOMATIC) libmesh_override;
125 
126  virtual void init (const NumericVector<T> & other,
127  const bool fast = false) libmesh_override;
128 
129  virtual NumericVector<T> & operator= (const T s) libmesh_override;
130 
131  virtual NumericVector<T> & operator= (const NumericVector<T> & v) libmesh_override;
132 
138  LaspackVector<T> & operator= (const LaspackVector<T> & v);
139 
140  virtual NumericVector<T> & operator= (const std::vector<T> & v) libmesh_override;
141 
142  virtual Real min () const libmesh_override;
143 
144  virtual Real max () const libmesh_override;
145 
146  virtual T sum () const libmesh_override;
147 
148  virtual Real l1_norm () const libmesh_override;
149 
150  virtual Real l2_norm () const libmesh_override;
151 
152  virtual Real linfty_norm () const libmesh_override;
153 
154  virtual numeric_index_type size () const libmesh_override;
155 
156  virtual numeric_index_type local_size() const libmesh_override;
157 
158  virtual numeric_index_type first_local_index() const libmesh_override;
159 
160  virtual numeric_index_type last_local_index() const libmesh_override;
161 
162  virtual T operator() (const numeric_index_type i) const libmesh_override;
163 
164  virtual NumericVector<T> & operator += (const NumericVector<T> & v) libmesh_override;
165 
166  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) libmesh_override;
167 
168  virtual NumericVector<T> & operator /= (NumericVector<T> & v) libmesh_override;
169 
170  virtual void reciprocal() libmesh_override;
171 
172  virtual void conjugate() libmesh_override;
173 
174  virtual void set (const numeric_index_type i, const T value) libmesh_override;
175 
176  virtual void add (const numeric_index_type i, const T value) libmesh_override;
177 
178  virtual void add (const T s) libmesh_override;
179 
180  virtual void add (const NumericVector<T> & v) libmesh_override;
181 
182  virtual void add (const T a, const NumericVector<T> & v) libmesh_override;
183 
188  using NumericVector<T>::add_vector;
189 
190  virtual void add_vector (const NumericVector<T> & v,
191  const SparseMatrix<T> & A) libmesh_override;
192 
193  virtual void add_vector_transpose (const NumericVector<T> & v,
194  const SparseMatrix<T> & A) libmesh_override;
195 
196  virtual void scale (const T factor) libmesh_override;
197 
198  virtual void abs() libmesh_override;
199 
200  virtual T dot(const NumericVector<T> & v) const libmesh_override;
201 
202  virtual void localize (std::vector<T> & v_local) const libmesh_override;
203 
204  virtual void localize (NumericVector<T> & v_local) const libmesh_override;
205 
206  virtual void localize (NumericVector<T> & v_local,
207  const std::vector<numeric_index_type> & send_list) const libmesh_override;
208 
209  virtual void localize (std::vector<T> & v_local,
210  const std::vector<numeric_index_type> & indices) const libmesh_override;
211 
212  virtual void localize (const numeric_index_type first_local_idx,
213  const numeric_index_type last_local_idx,
214  const std::vector<numeric_index_type> & send_list) libmesh_override;
215 
216  virtual void localize_to_one (std::vector<T> & v_local,
217  const processor_id_type proc_id=0) const libmesh_override;
218 
219  virtual void pointwise_mult (const NumericVector<T> & vec1,
220  const NumericVector<T> & vec2) libmesh_override;
221 
222  virtual void swap (NumericVector<T> & v) libmesh_override;
223 
224 private:
225 
230  QVector _vec;
231 
235  friend class LaspackLinearSolver<T>;
236 };
237 
238 
239 
240 //----------------------------------------------------------
241 // LaspackVector inline methods
242 template <typename T>
243 inline
244 LaspackVector<T>::LaspackVector (const Parallel::Communicator & comm,
245  const ParallelType ptype)
246  : NumericVector<T>(comm, ptype)
247 {
248  this->_type = ptype;
249 }
250 
251 
252 
253 template <typename T>
254 inline
256  const numeric_index_type n,
257  const ParallelType ptype)
258  : NumericVector<T>(comm, ptype)
259 {
260  this->init(n, n, false, ptype);
261 }
262 
263 
264 
265 template <typename T>
266 inline
268  const numeric_index_type n,
269  const numeric_index_type n_local,
270  const ParallelType ptype)
271  : NumericVector<T>(comm, ptype)
272 {
273  this->init(n, n_local, false, ptype);
274 }
275 
276 
277 
278 template <typename T>
279 inline
281  const numeric_index_type N,
282  const numeric_index_type n_local,
283  const std::vector<numeric_index_type> & ghost,
284  const ParallelType ptype)
285  : NumericVector<T>(comm, ptype)
286 {
287  this->init(N, n_local, ghost, false, ptype);
288 }
289 
290 
291 
292 template <typename T>
293 inline
295 {
296  this->clear ();
297 }
298 
299 
300 
301 template <typename T>
302 inline
304  const numeric_index_type libmesh_dbg_var(n_local),
305  const bool fast,
306  const ParallelType)
307 {
308  // Laspack vectors only for serial cases,
309  // but can provide a "parallel" vector on one processor.
310  libmesh_assert_equal_to (n, n_local);
311 
312  this->_type = SERIAL;
313 
314  // Clear initialized vectors
315  if (this->initialized())
316  this->clear();
317 
318  // create a sequential vector
319 
320  static int cnt = 0;
321  char foo[80];
322  std::sprintf(foo, "Vec-%d", cnt++);
323 
324  V_Constr(&_vec, const_cast<char *>(foo), n, Normal, _LPTrue);
325 
326  this->_is_initialized = true;
327 #ifndef NDEBUG
328  this->_is_closed = true;
329 #endif
330 
331  // Optionally zero out all components
332  if (fast == false)
333  this->zero ();
334 
335  return;
336 }
337 
338 
339 
340 template <typename T>
341 inline
343  const bool fast,
344  const ParallelType ptype)
345 {
346  this->init(n,n,fast,ptype);
347 }
348 
349 
350 template <typename T>
351 inline
353  const numeric_index_type n_local,
354  const std::vector<numeric_index_type> & libmesh_dbg_var(ghost),
355  const bool fast,
356  const ParallelType ptype)
357 {
358  libmesh_assert(ghost.empty());
359  this->init(n,n_local,fast,ptype);
360 }
361 
362 
363 
364 /* Default implementation for solver packages for which ghosted
365  vectors are not yet implemented. */
366 template <class T>
367 void LaspackVector<T>::init (const NumericVector<T> & other,
368  const bool fast)
369 {
370  this->init(other.size(),other.local_size(),fast,other.type());
371 }
372 
373 
374 
375 template <typename T>
376 inline
378 {
379  libmesh_assert (this->initialized());
380 
381 #ifndef NDEBUG
382  this->_is_closed = true;
383 #endif
384 }
385 
386 
387 
388 template <typename T>
389 inline
391 {
392  if (this->initialized())
393  {
394  V_Destr (&_vec);
395  }
396 
397  this->_is_initialized = false;
398 #ifndef NDEBUG
399  this->_is_closed = false;
400 #endif
401 }
402 
403 
404 
405 template <typename T> inline
407 {
408  libmesh_assert (this->initialized());
409  libmesh_assert (this->closed());
410 
411  V_SetAllCmp (&_vec, 0.);
412 }
413 
414 
415 
416 template <typename T>
417 inline
419 {
420  NumericVector<T> * cloned_vector = new LaspackVector<T>(this->comm());
421 
422  cloned_vector->init(*this);
423 
424  return UniquePtr<NumericVector<T>>(cloned_vector);
425 }
426 
427 
428 
429 template <typename T>
430 inline
432 {
433  NumericVector<T> * cloned_vector = new LaspackVector<T>(this->comm());
434 
435  cloned_vector->init(*this, true);
436 
437  *cloned_vector = *this;
438 
439  return UniquePtr<NumericVector<T>>(cloned_vector);
440 }
441 
442 
443 
444 template <typename T>
445 inline
447 {
448  libmesh_assert (this->initialized());
449 
450  return static_cast<numeric_index_type>(V_GetDim(const_cast<QVector*>(&_vec)));
451 }
452 
453 
454 
455 template <typename T>
456 inline
458 {
459  libmesh_assert (this->initialized());
460 
461  return this->size();
462 }
463 
464 
465 
466 template <typename T>
467 inline
469 {
470  libmesh_assert (this->initialized());
471 
472  return 0;
473 }
474 
475 
476 
477 template <typename T>
478 inline
480 {
481  libmesh_assert (this->initialized());
482 
483  return this->size();
484 }
485 
486 
487 
488 template <typename T>
489 inline
490 void LaspackVector<T>::set (const numeric_index_type i, const T value)
491 {
492  libmesh_assert (this->initialized());
493  libmesh_assert_less (i, this->size());
494 
495  V_SetCmp (&_vec, i+1, value);
496 
497 #ifndef NDEBUG
498  this->_is_closed = false;
499 #endif
500 }
501 
502 
503 
504 template <typename T>
505 inline
506 void LaspackVector<T>::add (const numeric_index_type i, const T value)
507 {
508  libmesh_assert (this->initialized());
509  libmesh_assert_less (i, this->size());
510 
511  V_AddCmp (&_vec, i+1, value);
512 
513 #ifndef NDEBUG
514  this->_is_closed = false;
515 #endif
516 }
517 
518 
519 
520 template <typename T>
521 inline
523 {
524  libmesh_assert (this->initialized());
525  libmesh_assert ( ((i >= this->first_local_index()) &&
526  (i < this->last_local_index())) );
527 
528 
529  return static_cast<T>(V_GetCmp(const_cast<QVector*>(&_vec), i+1));
530 }
531 
532 
533 
534 template <typename T>
535 inline
537 {
538  LaspackVector<T> & v = cast_ref<LaspackVector<T> &>(other);
539 
540  // This is all grossly dependent on Laspack version...
541 
542  std::swap(_vec.Name, v._vec.Name);
543  std::swap(_vec.Dim, v._vec.Dim);
544  std::swap(_vec.Instance, v._vec.Instance);
545  std::swap(_vec.LockLevel, v._vec.LockLevel);
546  std::swap(_vec.Multipl, v._vec.Multipl);
547  std::swap(_vec.OwnData, v._vec.OwnData);
548 
549  // This should still be O(1), since _vec.Cmp is just a pointer to
550  // data on the heap
551 
552  std::swap(_vec.Cmp, v._vec.Cmp);
553 }
554 
555 
556 } // namespace libMesh
557 
558 
559 #endif // #ifdef LIBMESH_HAVE_LASPACK
560 #endif // LIBMESH_LASPACK_VECTOR_H
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:281
double abs(double a)
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
virtual numeric_index_type size() const =0
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
uint8_t processor_id_type
Definition: id_types.h:99
Numeric vector.
Definition: dof_map.h:66
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
const Number zero
.
Definition: libmesh.h:178
long double max(long double a, double b)
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
This class provides an interface to Laspack iterative solvers that is compatible with the libMesh Lin...
ParallelType
Defines an enum for parallel data structure types.
dof_id_type numeric_index_type
Definition: id_types.h:92
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:255
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual numeric_index_type local_size() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
X_input swap(X_system)
static PetscErrorCode Mat * A
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
static const bool value
Definition: xdr_io.C:108
ParallelType type() const
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:274
long double min(long double a, double b)