laspack_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
44 // Forward declarations
45 template <typename T> class LaspackLinearSolver;
46 template <typename T> class SparseMatrix;
47 
48 
56 template <typename T>
57 class LaspackVector : 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,
83  const ParallelType = AUTOMATIC);
84 
90  LaspackVector (const Parallel::Communicator &comm,
91  const numeric_index_type N,
93  const std::vector<numeric_index_type>& ghost,
94  const ParallelType = AUTOMATIC);
95 
100  ~LaspackVector ();
101 
105  void close ();
106 
110  void clear ();
111 
116  void zero ();
117 
123  virtual AutoPtr<NumericVector<T> > zero_clone () const;
124 
128  AutoPtr<NumericVector<T> > clone () const;
129 
143  void init (const numeric_index_type N,
145  const bool fast=false,
146  const ParallelType ptype=AUTOMATIC);
147 
151  void init (const numeric_index_type N,
152  const bool fast=false,
153  const ParallelType ptype=AUTOMATIC);
154 
159  void init (const numeric_index_type /*N*/,
160  const numeric_index_type /*n_local*/,
161  const std::vector<numeric_index_type>& /*ghost*/,
162  const bool /*fast*/ = false,
163  const ParallelType = AUTOMATIC);
164 
169  virtual void init (const NumericVector<T>& other,
170  const bool fast = false);
171 
175  NumericVector<T> & operator= (const T s);
176 
180  NumericVector<T> & operator= (const NumericVector<T> &V);
181 
185  LaspackVector<T> & operator= (const LaspackVector<T> &V);
186 
190  NumericVector<T> & operator= (const std::vector<T> &v);
191 
197  Real min () const;
198 
204  Real max () const;
205 
209  T sum () const;
210 
215  Real l1_norm () const;
216 
222  Real l2_norm () const;
223 
229  Real linfty_norm () const;
230 
238  numeric_index_type size () const;
239 
245 
251 
257 
261  T operator() (const numeric_index_type i) const;
262 
267  NumericVector<T> & operator += (const NumericVector<T> &V);
268 
273  NumericVector<T> & operator -= (const NumericVector<T> &V);
274 
278  virtual NumericVector<T> & operator /= (NumericVector<T> & v);
279 
283  virtual void reciprocal();
284 
289  virtual void conjugate();
290 
294  void set (const numeric_index_type i, const T value);
295 
299  void add (const numeric_index_type i, const T value);
300 
306  void add (const T s);
307 
313  void add (const NumericVector<T>& V);
314 
320  void add (const T a, const NumericVector<T>& v);
321 
327  void add_vector (const std::vector<T>& v,
328  const std::vector<numeric_index_type>& dof_indices);
329 
336  void add_vector (const NumericVector<T>& V,
337  const std::vector<numeric_index_type>& dof_indices);
338 
343  void add_vector (const NumericVector<T> &,
344  const SparseMatrix<T> &);
345 
352  void add_vector (const DenseVector<T>& V,
353  const std::vector<numeric_index_type>& dof_indices);
354 
359  void add_vector_transpose (const NumericVector<T> &,
360  const SparseMatrix<T> &);
361 
366  virtual void insert (const std::vector<T>& v,
367  const std::vector<numeric_index_type>& dof_indices);
368 
375  virtual void insert (const NumericVector<T>& V,
376  const std::vector<numeric_index_type>& dof_indices);
377 
383  virtual void insert (const DenseVector<T>& V,
384  const std::vector<numeric_index_type>& dof_indices);
385 
391  virtual void insert (const DenseSubVector<T>& V,
392  const std::vector<numeric_index_type>& dof_indices);
393 
398  void scale (const T factor);
399 
404  virtual void abs();
405 
409  virtual T dot(const NumericVector<T>& V) const;
410 
415  void localize (std::vector<T>& v_local) const;
416 
421  void localize (NumericVector<T>& v_local) const;
422 
428  void localize (NumericVector<T>& v_local,
429  const std::vector<numeric_index_type>& send_list) const;
430 
435  void localize (const numeric_index_type first_local_idx,
436  const numeric_index_type last_local_idx,
437  const std::vector<numeric_index_type>& send_list);
438 
439 
446  void localize_to_one (std::vector<T>& v_local,
447  const processor_id_type proc_id=0) const;
448 
453  virtual void pointwise_mult (const NumericVector<T>& vec1,
454  const NumericVector<T>& vec2);
455 
459  virtual void swap (NumericVector<T> &v);
460 
461  private:
462 
467  QVector _vec;
468 
472  friend class LaspackLinearSolver<T>;
473 };
474 
475 
476 
477 //----------------------- ----------------------------------
478 // LaspackVector inline methods
479 template <typename T>
480 inline
482  const ParallelType ptype)
483  : NumericVector<T>(comm, ptype)
484 {
485  this->_type = ptype;
486 }
487 
488 
489 
490 template <typename T>
491 inline
493  const numeric_index_type n,
494  const ParallelType ptype)
495  : NumericVector<T>(comm, ptype)
496 {
497  this->init(n, n, false, ptype);
498 }
499 
500 
501 
502 template <typename T>
503 inline
505  const numeric_index_type n,
507  const ParallelType ptype)
508  : NumericVector<T>(comm, ptype)
509 {
510  this->init(n, n_local, false, ptype);
511 }
512 
513 
514 
515 template <typename T>
516 inline
518  const numeric_index_type N,
520  const std::vector<numeric_index_type>& ghost,
521  const ParallelType ptype)
522  : NumericVector<T>(comm, ptype)
523 {
524  this->init(N, n_local, ghost, false, ptype);
525 }
526 
527 
528 
529 template <typename T>
530 inline
532 {
533  this->clear ();
534 }
535 
536 
537 
538 template <typename T>
539 inline
541  const numeric_index_type libmesh_dbg_var(n_local),
542  const bool fast,
543  const ParallelType)
544 {
545  // Laspack vectors only for serial cases,
546  // but can provide a "parallel" vector on one processor.
547  libmesh_assert_equal_to (n, n_local);
548 
549  this->_type = SERIAL;
550 
551  // Clear initialized vectors
552  if (this->initialized())
553  this->clear();
554 
555  // create a sequential vector
556 
557  static int cnt = 0;
558  char foo[80];
559  std::sprintf(foo, "Vec-%d", cnt++);
560 
561  V_Constr(&_vec, const_cast<char*>(foo), n, Normal, _LPTrue);
562 
563  this->_is_initialized = true;
564 #ifndef NDEBUG
565  this->_is_closed = true;
566 #endif
567 
568  // Optionally zero out all components
569  if (fast == false)
570  this->zero ();
571 
572  return;
573 }
574 
575 
576 
577 template <typename T>
578 inline
580  const bool fast,
581  const ParallelType ptype)
582 {
583  this->init(n,n,fast,ptype);
584 }
585 
586 
587 template <typename T>
588 inline
591  const std::vector<numeric_index_type>& libmesh_dbg_var(ghost),
592  const bool fast,
593  const ParallelType ptype)
594 {
595  libmesh_assert(ghost.empty());
596  this->init(n,n_local,fast,ptype);
597 }
598 
599 
600 
601 /* Default implementation for solver packages for which ghosted
602  vectors are not yet implemented. */
603 template <class T>
605  const bool fast)
606 {
607  this->init(other.size(),other.local_size(),fast,other.type());
608 }
609 
610 
611 
612 template <typename T>
613 inline
615 {
616  libmesh_assert (this->initialized());
617 
618 #ifndef NDEBUG
619  this->_is_closed = true;
620 #endif
621 }
622 
623 
624 
625 template <typename T>
626 inline
628 {
629  if (this->initialized())
630  {
631  V_Destr (&_vec);
632  }
633 
634  this->_is_initialized = false;
635 #ifndef NDEBUG
636  this->_is_closed = false;
637 #endif
638 }
639 
640 
641 
642 template <typename T> inline
644 {
645  libmesh_assert (this->initialized());
646  libmesh_assert (this->closed());
647 
648  V_SetAllCmp (&_vec, 0.);
649 }
650 
651 
652 
653 template <typename T>
654 inline
656 {
657  AutoPtr<NumericVector<T> > cloned_vector
658  (new LaspackVector<T>(this->comm()));
659 
660  cloned_vector->init(*this);
661 
662  return cloned_vector;
663 }
664 
665 
666 
667 template <typename T>
668 inline
670 {
671  AutoPtr<NumericVector<T> > cloned_vector
672  (new LaspackVector<T>(this->comm()));
673 
674  cloned_vector->init(*this, true);
675 
676  *cloned_vector = *this;
677 
678  return cloned_vector;
679 }
680 
681 
682 
683 template <typename T>
684 inline
686 {
687  libmesh_assert (this->initialized());
688 
689  return static_cast<numeric_index_type>(V_GetDim(const_cast<QVector*>(&_vec)));
690 }
691 
692 
693 
694 template <typename T>
695 inline
697 {
698  libmesh_assert (this->initialized());
699 
700  return this->size();
701 }
702 
703 
704 
705 template <typename T>
706 inline
708 {
709  libmesh_assert (this->initialized());
710 
711  return 0;
712 }
713 
714 
715 
716 template <typename T>
717 inline
719 {
720  libmesh_assert (this->initialized());
721 
722  return this->size();
723 }
724 
725 
726 
727 template <typename T>
728 inline
729 void LaspackVector<T>::set (const numeric_index_type i, const T value)
730 {
731  libmesh_assert (this->initialized());
732  libmesh_assert_less (i, this->size());
733 
734  V_SetCmp (&_vec, i+1, value);
735 
736 #ifndef NDEBUG
737  this->_is_closed = false;
738 #endif
739 }
740 
741 
742 
743 template <typename T>
744 inline
745 void LaspackVector<T>::add (const numeric_index_type i, const T value)
746 {
747  libmesh_assert (this->initialized());
748  libmesh_assert_less (i, this->size());
749 
750  V_AddCmp (&_vec, i+1, value);
751 
752 #ifndef NDEBUG
753  this->_is_closed = false;
754 #endif
755 }
756 
757 
758 
759 template <typename T>
760 inline
762 {
763  libmesh_assert (this->initialized());
764  libmesh_assert ( ((i >= this->first_local_index()) &&
765  (i < this->last_local_index())) );
766 
767 
768  return static_cast<T>(V_GetCmp(const_cast<QVector*>(&_vec), i+1));
769 }
770 
771 
772 
773 template <typename T>
774 inline
776 {
777  LaspackVector<T>& v = libmesh_cast_ref<LaspackVector<T>&>(other);
778 
779  // This is all grossly dependent on Laspack version...
780 
781  std::swap(_vec.Name, v._vec.Name);
782  std::swap(_vec.Dim, v._vec.Dim);
783  std::swap(_vec.Instance, v._vec.Instance);
784  std::swap(_vec.LockLevel, v._vec.LockLevel);
785  std::swap(_vec.Multipl, v._vec.Multipl);
786  std::swap(_vec.OwnData, v._vec.OwnData);
787 
788  // This should still be O(1), since _vec.Cmp is just a pointer to
789  // data on the heap
790 
791  std::swap(_vec.Cmp, v._vec.Cmp);
792 }
793 
794 
795 } // namespace libMesh
796 
797 
798 #endif // #ifdef LIBMESH_HAVE_LASPACK
799 #endif // LIBMESH_LASPACK_VECTOR_H

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo