• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

MLAPI_MultiVector.h

Go to the documentation of this file.
00001 #ifndef ML_MULTIVECTOR_H
00002 #define ML_MULTIVECTOR_H
00003 
00013 /* ******************************************************************** */
00014 /* See the file COPYRIGHT for a complete copyright notice, contact      */
00015 /* person and disclaimer.                                               */        
00016 /* ******************************************************************** */
00017 
00018 //#include "ml_lapack.h"
00019 #include "MLAPI_Error.h"
00020 #include "MLAPI_BaseObject.h"
00021 #include "MLAPI_Space.h"
00022 #include "MLAPI_CompObject.h"
00023 #include "MLAPI_TimeObject.h"
00024 #include "MLAPI_BaseLinearCombination.h"
00025 #include "Teuchos_RefCountPtr.hpp"
00026 #include "Teuchos_BLAS_wrappers.hpp"
00027 #include "Teuchos_LAPACK_wrappers.hpp"
00028 #include <iomanip>
00029 #include <vector>
00030 #include <algorithm>
00031 #include "ml_utils.h"
00032 
00033 namespace MLAPI {
00034 
00035 class DoubleVector {
00036 
00037 public:
00038   DoubleVector(const int size)
00039   {
00040     ownership_ = true;
00041     ptr_ = new double[size];
00042   }
00043 
00044   DoubleVector(double* ptr)
00045   {
00046     ownership_ = false;
00047     ptr_ = ptr;
00048   }
00049 
00050   inline double& operator[](const int i)
00051   {
00052     return(ptr_[i]);
00053   }
00054 
00055   inline const double& operator[](const int i) const
00056   {
00057     return(ptr_[i]);
00058   }
00059 
00060   ~DoubleVector()
00061   {
00062     if (ownership_)
00063       delete[] ptr_;
00064   }
00065 
00066   double* Values()
00067   {
00068     return(ptr_);
00069   }
00070 
00071   const double* Values() const
00072   {
00073     return(ptr_);
00074   }
00075 
00076 private:
00077   DoubleVector(const DoubleVector& rhs)
00078   {}
00079 
00080   DoubleVector& operator=(const DoubleVector& rhs)
00081   {
00082     return(*this);
00083   }
00084   
00085   double* ptr_;
00086   bool    ownership_;
00087 };
00088 
00089 class MultiVector;
00090 
00091     
00102 class MultiVector : public BaseObject, public CompObject, public TimeObject {
00103 
00104 public:
00105 
00107 
00109   MultiVector() 
00110   { 
00111     NumVectors_ = 0;
00112   }
00113 
00115   MultiVector(const Space& VectorSpace, const int NumVectors = 1,
00116                bool SetToZero = true)
00117   {
00118     StackPush();
00119 
00120     NumVectors_  = NumVectors;
00121     VectorSpace_ = VectorSpace;
00122     SetRCPLength(NumVectors);
00123     if (GetMyLength()) {
00124       for (int v = 0 ; v < NumVectors ; ++v)
00125         SetRCPValues(Teuchos::rcp(new DoubleVector(GetMyLength())), v);
00126 
00127       if (SetToZero)
00128         *this = 0.0;
00129     }
00130 
00131     StackPop();
00132   }
00133 
00135   MultiVector(const Space& VectorSpace, double** Values,
00136               const int NumVectors = 1)
00137   {
00138     StackPush();
00139 
00140     NumVectors_  = NumVectors;
00141     VectorSpace_ = VectorSpace;
00142     SetRCPLength(GetNumVectors());
00143     for (int v = 0 ; v < GetNumVectors() ; ++v)
00144       SetRCPValues(Teuchos::rcp(new DoubleVector(Values[v])), v);
00145 
00146     StackPop();
00147   }
00148 
00150   MultiVector(const Space& VectorSpace, Teuchos::RefCountPtr<DoubleVector> RCPValues)
00151   {
00152     StackPush();
00153 
00154     NumVectors_  = 1;
00155     VectorSpace_ = VectorSpace;
00156     SetRCPLength(GetNumVectors());
00157     SetRCPValues(RCPValues, 0);
00158 
00159     StackPop();
00160   }
00161 
00163   MultiVector(const Space& VectorSpace, 
00164               std::vector<Teuchos::RefCountPtr<DoubleVector> > RCPValues)
00165   {
00166     StackPush();
00167 
00168     NumVectors_  = (int)RCPValues.size();
00169     VectorSpace_ = VectorSpace;
00170     SetRCPLength(GetNumVectors());
00171     for (int v = 0 ; v < GetNumVectors() ; ++v)
00172       SetRCPValues(RCPValues[v], v);
00173 
00174     StackPop();
00175   }
00176 
00178   MultiVector(const MultiVector& rhs)
00179   {
00180     NumVectors_  = rhs.GetNumVectors();
00181     VectorSpace_ = rhs.GetVectorSpace();
00182     SetRCPLength(GetNumVectors());
00183     for (int v = 0 ; v < GetNumVectors() ; ++v)
00184       SetRCPValues(rhs.GetRCPValues(v), v);
00185   }
00186 
00188   ~MultiVector() 
00189   {
00190     for (int v = 0 ; v < GetNumVectors() ; ++v)
00191       SetRCPValues(Teuchos::null, v);
00192   }
00193 
00194   // @}
00195   // @{ \name Reshape methods
00196   
00198   void Reshape()
00199   {
00200     StackPush();
00201 
00202     for (int v = 0 ; v < GetNumVectors() ; ++v)
00203       SetRCPValues(Teuchos::null, v);
00204     SetRCPLength(0);
00205     GetVectorSpace().Reshape();
00206     NumVectors_ = 0;
00207 
00208     StackPop();
00209   }
00210 
00212   void Reshape(const Space& S, const int NumVectors = 1,
00213                const bool SetToZero = true)
00214   {
00215     StackPush();
00216 
00217     NumVectors_ = NumVectors;
00218     VectorSpace_ = S;
00219     SetRCPLength(GetNumVectors());
00220     for (int v = 0 ; v < GetNumVectors() ; ++v) {
00221       if (GetMyLength())
00222         SetRCPValues(Teuchos::rcp(new DoubleVector(GetMyLength())), v);
00223       else
00224         SetRCPValues(Teuchos::null, v);
00225     }
00226 
00227     if (SetToZero) *this = 0.0;
00228 
00229     StackPop();
00230   }
00231 
00233   void Append(const int NumVectors = 1, const bool SetToZero = true)
00234   {
00235     int n = GetMyLength();
00236 
00237     for (int v = 0 ; v < NumVectors ; ++v) {
00238       if (GetMyLength()) {
00239         RCPValues_.push_back(Teuchos::rcp(new DoubleVector(n)));
00240       }
00241       else
00242         RCPValues_.push_back(Teuchos::null);
00243 
00244       ++NumVectors_;
00245 
00246       if (SetToZero) {
00247         Update(0.0, GetNumVectors() - 1);
00248       }
00249     }
00250   }
00251    
00253   void Append(MultiVector rhs)
00254   {
00255     StackPush();
00256 
00257     CheckSpaces(rhs);
00258 
00259     for (int v = 0 ; v < rhs.GetNumVectors() ; ++v) {
00260       RCPValues_.push_back(rhs.GetRCPValues(v));
00261       ++NumVectors_;
00262     }
00263 
00264     StackPop();
00265   }
00266 
00268   void Delete(const int v)
00269   {
00270     StackPush();
00271 
00272     CheckVector(v);
00273 
00274     std::vector<Teuchos::RefCountPtr<DoubleVector> > NewList;
00275 
00276     for (int i = 0 ; i < GetNumVectors() ; ++i) {
00277       if (i != v)
00278         NewList.push_back(GetRCPValues(i));
00279     }
00280 
00281     RCPValues_ = NewList;
00282 
00283     NumVectors_--;
00284 
00285     StackPop();
00286   }
00287 
00288   // @}
00289   // @{ \name Overloaded operators
00290 
00292   MultiVector& operator=(double rhs) 
00293   {
00294     StackPush();
00295 
00296     for (int v = 0 ; v < GetNumVectors() ; ++v)
00297       for (int i = 0 ; i < GetMyLength() ; ++i)
00298         GetValues(v)[i] = rhs;
00299 
00300     StackPop();
00301 
00302     return(*this);
00303   }
00304 
00306   MultiVector& operator=(const MultiVector& rhs) 
00307   {
00308     StackPush();
00309 
00310     if (this != &rhs) {
00311       NumVectors_  = rhs.GetNumVectors();
00312       VectorSpace_ = rhs.GetVectorSpace();
00313       SetRCPLength(GetNumVectors());
00314       for (int v = 0 ; v < GetNumVectors() ; ++v)
00315         SetRCPValues(rhs.GetRCPValues(v), v);
00316       SetLabel(rhs.GetLabel());
00317     }
00318 
00319     StackPop();
00320 
00321     return(*this);
00322   }
00323 
00325   MultiVector& operator=(const BaseLinearCombination& rhs)
00326   {
00327     StackPush();
00328 
00329     rhs.Set(*this);
00330 
00331     StackPop();
00332 
00333     return(*this);
00334   }
00335 
00336   bool operator==(const MultiVector& rhs) const
00337   {
00338     return(false);
00339   }
00340 
00342   MultiVector& operator=(const string& Label)
00343   {
00344     SetLabel(Label);
00345     return(*this);
00346   }
00347 
00349   inline const double& operator() (const int i) const
00350   {
00351     CheckSingleVector();
00352     CheckEntry(i);
00353     CheckVector(0);
00354 
00355     return(GetValues(0)[i]);
00356   }
00357 
00359   inline double& operator() (const int i) 
00360   {
00361     CheckSingleVector();
00362     CheckEntry(i);
00363     CheckVector(0);
00364 
00365     return(GetValues(0)[i]);
00366   }
00367 
00369   inline const double& operator()(const int i, const int v) const 
00370   {
00371     CheckEntry(i);
00372     CheckVector(v);
00373 
00374     return(GetValues(v)[i]);
00375   }
00376 
00378   inline double& operator()(const int i, const int v) 
00379   {
00380     CheckEntry(i);
00381     CheckVector(v);
00382 
00383     return(GetValues(v)[i]);
00384   }
00385 
00386   // @}
00387   // @{ \name Set and Get methods
00388   
00390   inline const Space& GetVectorSpace() const 
00391   {
00392     return(VectorSpace_);
00393   }
00394 
00396   inline Space& GetVectorSpace() 
00397   {
00398     return(VectorSpace_);
00399   }
00401   inline int GetNumVectors() const
00402   {
00403     return(NumVectors_);
00404   }
00405 
00407   inline int GetMyLength() const
00408   {
00409     return(VectorSpace_.GetNumMyElements());
00410   }
00411 
00413   inline int GetGlobalLength() const
00414   {
00415     return(VectorSpace_.GetNumGlobalElements());
00416   }
00417 
00419   inline double* GetValues(const int v)
00420   {
00421     return(RCPValues_[v].get()->Values());
00422   }
00423 
00425   inline const double* GetValues(const int v) const
00426   {
00427     return(RCPValues_[v].get()->Values());
00428   }
00429   
00431   inline Teuchos::RefCountPtr<DoubleVector>& GetRCPValues(const int v) 
00432   {
00433     CheckVector(v);
00434 
00435     return(RCPValues_[v]);
00436   }
00437 
00439   inline const Teuchos::RefCountPtr<DoubleVector>& GetRCPValues(const int v) const
00440   {
00441     CheckVector(v);
00442 
00443     return(RCPValues_[v]);
00444   }
00445 
00447   inline void SetRCPValues(const Teuchos::RefCountPtr<DoubleVector>& RCPValues,
00448                            const int v)
00449   {
00450     CheckVector(v);
00451 
00452     RCPValues_[v] = RCPValues;
00453   }
00454 
00455   // @}
00456   // @{ \name Mathematical methods
00457   
00459   void Update(const double alpha, int v = -1)
00460   {
00461     if (v == -1) {
00462       CheckSingleVector();
00463       v = 0;
00464     }
00465 
00466     if (v >= GetNumVectors())
00467       ML_THROW("Requested wrong vector, " + GetString(v) +
00468                " while NumVectors = " + GetString(GetNumVectors()), -1);
00469 
00470     for (int i = 0 ; i < GetMyLength() ; ++i)
00471       GetValues(v)[i] = alpha;
00472   }
00473 
00475   void Update(const MultiVector& rhs)
00476   {
00477     ResetTimer();
00478     StackPush();
00479 
00480     int n = GetMyLength();
00481     if (n == 0) return;
00482 
00483     CheckSpaces(rhs);
00484     CheckNumVectors(rhs.GetNumVectors());
00485 
00486     int incr = 1;
00487     for (int v = 0 ; v < GetNumVectors() ; ++v)
00488       // copy rhs into this
00489       DCOPY_F77(&n, (double*)rhs.GetValues(v), &incr, GetValues(v), &incr);
00490 
00491     StackPop();
00492     UpdateTime();
00493   }
00494   
00496   void Update(double alpha, const MultiVector& rhs)
00497   {
00498     ResetTimer();
00499     StackPush();
00500 
00501     int n = GetMyLength();
00502     if (n == 0) return;
00503 
00504     CheckSpaces(rhs);
00505     CheckNumVectors(rhs.GetNumVectors());
00506 
00507     int incr = 1;
00508     for (int v = 0 ; v < GetNumVectors() ; ++v) {
00509       // copy rhs into this
00510       DCOPY_F77(&n, (double*)rhs.GetValues(v), &incr, GetValues(v), &incr);
00511       // scale this
00512       DSCAL_F77(&n, &alpha, GetValues(v), &incr);
00513     }
00514 
00515     StackPop();
00516     UpdateFlops(1.0 * GetNumVectors() * GetGlobalLength());
00517     UpdateTime();
00518 
00519   }
00520 
00522   void Update(double alpha, const MultiVector& x,
00523               double beta,  const MultiVector& y)
00524   {
00525     ResetTimer();
00526     StackPush();
00527 
00528     int n = GetMyLength();
00529     if (n == 0) return;
00530 
00531     CheckSpaces(x);
00532     CheckSpaces(y);
00533     CheckNumVectors(x.GetNumVectors());
00534     CheckNumVectors(y.GetNumVectors());
00535 
00536     int incr = 1;
00537     for (int v = 0 ; v < GetNumVectors() ; ++v)
00538       // copy rhs into this
00539       DCOPY_F77(&n, (double*)x.GetValues(v), &incr, GetValues(v), &incr);
00540 
00541     StackPop();
00542     Update(beta,y,alpha);
00543     UpdateTime();
00544 
00545   }
00546 
00548   void Update(double alpha, const MultiVector& rhs, double beta)
00549   {
00550     ResetTimer();
00551     StackPush();
00552 
00553     int n = GetMyLength();
00554     if (n == 0) return;
00555 
00556     CheckSpaces(rhs);
00557     CheckNumVectors(rhs.GetNumVectors());
00558 
00559     for (int v = 0 ; v < GetNumVectors() ; ++v) 
00560     {
00561       double* ptr_this = GetValues(v);
00562       double* ptr_rhs  = (double*)(rhs.GetValues(v));
00563 
00564       if (alpha == 1.0 && beta == 1.0)
00565       {
00566         for (int i = 0 ; i < GetMyLength() ; ++i)
00567           ptr_this[i] += ptr_rhs[i];
00568         UpdateFlops(GetGlobalLength());
00569       }
00570       else if (alpha == 1.0 && beta == 0.0)
00571       {
00572         for (int i = 0 ; i < GetMyLength() ; ++i)
00573           ptr_this[i] = ptr_rhs[i];
00574       }
00575       else if (alpha == 0.0 && beta == 1.0)
00576       {
00577         // do nothing here
00578         if (false)
00579           cout << "blablablaaaaa" << endl;
00580       }
00581       else if (alpha == 1.0 && beta == -1.0)
00582       {
00583         for (int i = 0 ; i < GetMyLength() ; ++i)
00584           ptr_this[i] = ptr_rhs[i] - ptr_this[i];
00585         UpdateFlops(GetGlobalLength()); 
00586       }
00587       else if (alpha == -1.0 && beta == 1.0)
00588       {
00589         for (int i = 0 ; i < GetMyLength() ; ++i)
00590           ptr_this[i] -= ptr_rhs[i];
00591         UpdateFlops(GetGlobalLength()); 
00592       }
00593       else
00594       {
00595         for (int i = 0 ; i < GetMyLength() ; ++i)
00596           ptr_this[i] = ptr_rhs[i] * alpha + ptr_this[i] * beta;
00597         UpdateFlops(3.0 * GetGlobalLength()); 
00598       }
00599     }
00600 
00601     StackPop();
00602     UpdateTime();
00603   }
00604 
00605 #if 0
00606   void Add(double alpha)
00607   {
00608     for (int v = 0 ; v < GetNumVectors() ; ++v)
00609       for (int i = 0 ; i < GetMyLength() ; ++i)
00610         GetValues(v)[i] += alpha;
00611   }
00612 #endif
00613 
00615   inline double DotProduct(const MultiVector& rhs, int v = -1) const 
00616   {
00617     ResetTimer();
00618     StackPush();
00619 
00620     if (rhs.GetVectorSpace() != GetVectorSpace()) {
00621       ML_THROW("rhs.GetVectorSpace() is not equal to this->GetVectorSpace()", -1);
00622     }
00623 
00624     CheckNumVectors(rhs.GetNumVectors());
00625     
00626     if (v == -1) {
00627       CheckSingleVector();
00628       v = 0;
00629     }
00630 
00631     double MyResult = 0.0;
00632     double Result   = 0.0;
00633     int n           = GetMyLength();
00634     int incr        = 1;
00635     double* ptr     = (double*)GetValues(v);
00636     double* rhs_ptr = (double*)rhs.GetValues(v);
00637     MyResult        = DDOT_F77(&n, ptr, &incr, rhs_ptr, &incr);
00638     Result          = ML_Comm_GsumDouble(GetML_Comm(),MyResult);
00639 
00640     StackPop();
00641     UpdateFlops(2.0 * GetGlobalLength()); // DDOT
00642     UpdateTime();
00643 
00644     return(Result);
00645   }
00646 
00648   inline double Norm2(int v = -1) const 
00649   {
00650     ResetTimer();
00651     StackPush();
00652 
00653     if (v == -1) {
00654       CheckSingleVector();
00655       v = 0;
00656     }
00657 
00658     double MyResult = 0.0;
00659     double Result   = 0.0;
00660     int n           = GetMyLength();
00661     int incr        = 1;
00662     double* ptr     = (double*)GetValues(v);
00663     MyResult        = DDOT_F77(&n, ptr, &incr, ptr, &incr);
00664     Result          = ML_Comm_GsumDouble(GetML_Comm(),MyResult);
00665     
00666     StackPop();
00667     UpdateFlops(2.0 * GetGlobalLength()); // DDOT
00668     UpdateTime();
00669 
00670     return(sqrt(Result));
00671   }
00672 
00674   inline double NormInf(int v = -1) const 
00675   {
00676     ResetTimer();
00677     StackPush();
00678 
00679     if (v == -1) {
00680       CheckSingleVector();
00681       v = 0;
00682     }
00683 
00684     double MyResult = 0.0;
00685     double Result   = 0.0;
00686     int n           = GetMyLength();
00687     double* ptr     = (double*)GetValues(v);
00688     int incr        = 1;
00689     int i           = IDAMAX_F77(&n, ptr, &incr);
00690     MyResult        = fabs(ptr[i - 1]);
00691     // FIXME: delete below
00692     /*
00693     for (int i = 0 ; i < n ; ++i)
00694       if (MyResult < fabs(ptr[i])) 
00695         MyResult = fabs(ptr[i]);
00696         */
00697 
00698     Result          = ML_Comm_GmaxDouble(GetML_Comm(),MyResult);
00699 
00700     StackPop();
00701     UpdateTime();
00702     return(Result);
00703   }
00704 
00706   inline double NormOne(int v = -1) const 
00707   {
00708     ResetTimer();
00709     StackPush();
00710 
00711     if (v == -1) {
00712       CheckSingleVector();
00713       v = 0;
00714     }
00715 
00716     double MyResult = 0.0;
00717     double Result   = 0.0;
00718     double* ptr     = (double*)GetValues(v);
00719     for (int i = 0 ; i < GetMyLength() ; ++i)
00720       MyResult += fabs(ptr[i]);
00721 
00722     Result          = ML_Comm_GsumDouble(GetML_Comm(),MyResult);
00723 
00724     StackPop();
00725     UpdateTime();
00726     return(Result);
00727   }
00728 
00730   inline void Reciprocal(int v = -1) 
00731   {
00732     ResetTimer();
00733     StackPush();
00734     
00735     if (v == -1) {
00736       CheckSingleVector();
00737       v = 0;
00738     }
00739 
00740     for (int i = 0 ; i < GetMyLength() ; ++i) {
00741       if (GetValues(v)[i] != 0.0)
00742         GetValues(v)[i] = 1.0 / GetValues(v)[i];
00743     }
00744 
00745     StackPop();
00746 
00747     UpdateFlops(1.0 * GetGlobalLength());
00748     UpdateTime();
00749   }
00750 
00752   inline void Scale(const double Factor, int v = -1) 
00753   {
00754     ResetTimer();
00755     StackPush();
00756 
00757     if (v == -1) {
00758       CheckSingleVector();
00759       v = 0;
00760     }
00761 
00762     int n = GetMyLength();
00763     if (n == 0) return;
00764 
00765     int incr = 1;
00766     DSCAL_F77(&n, (double*)&Factor, GetValues(v), &incr);
00767 
00768     StackPop();
00769 
00770     UpdateFlops(1.0 * GetGlobalLength()); 
00771     UpdateTime();
00772   }
00773 
00774   // @}
00775   // @{ \name Miscellanous methods
00776 
00778   inline void Random(int v = -1) 
00779   {
00780     ResetTimer();
00781     StackPush();
00782 
00783     if (v == -1) {
00784       CheckSingleVector();
00785       v = 0;
00786     }
00787 
00788     ML_random_vec(GetValues(v),GetMyLength(),MLAPI::GetML_Comm());
00789 
00790     StackPop();
00791     UpdateTime();
00792     return;
00793   }
00794 
00796   void Sort(int v = -1, const bool IsIncreasing = false)
00797   {
00798     ResetTimer();
00799     StackPush();
00800 
00801     if (v == -1) {
00802       CheckSingleVector();
00803       v = 0;
00804     }
00805 
00806     CheckVector(v);
00807 
00808     std::vector<double> tmp(GetMyLength());
00809     for (int i = 0 ; i < GetMyLength() ; ++i)
00810       tmp[i] = (*this)(i, v);
00811 
00812     if (IsIncreasing)
00813       std::sort(tmp.begin(), tmp.end(), std::greater<double>());
00814     else
00815       std::sort(tmp.begin(), tmp.end());
00816 
00817     for (int i = 0 ; i < GetMyLength() ; ++i)
00818       (*this)(i,v) = tmp[i];
00819 
00820     StackPop();
00821     UpdateTime();
00822   }
00823 
00825   virtual std::ostream& Print(std::ostream& os,
00826                               const bool verbose = true) const
00827   {
00828     ResetTimer();
00829     StackPush();
00830 
00831     if (GetMyPID() == 0) {
00832       os << endl;
00833       os << "*** MLAPI::MultiVector ***" << endl;
00834       os << "Label             = " << GetLabel() << endl;
00835       os << "Local length      = " << GetMyLength() << endl;
00836       os << "Global length     = " << GetGlobalLength() << endl;
00837       os << "Number of vectors = " << GetNumVectors() << endl;
00838       os << "Flop count        = " << GetFlops() << endl;
00839       os << "Cumulative time   = " << GetTime() << endl;
00840       if (GetTime() != 0.0)
00841         os << "MFlops rate       = " << 1.0e-6 * GetFlops() / GetTime() << endl;
00842       else
00843         os << "MFlops rate       = 0.0" << endl;
00844       os << endl << endl;
00845     }
00846 
00847     if (verbose) {
00848 
00849       if (GetMyPID() == 0) {
00850         os.width(10);
00851         os << "ProcID";
00852         os.width(20);
00853         os << "LID";
00854         os.width(20);
00855         os << "GID";
00856         for (int v = 0 ; v < GetNumVectors() ; ++v) {
00857           os.width(20);
00858           os << "value " << v;
00859         }
00860         os << endl << endl;
00861       }
00862       
00863       for (int iproc = 0 ; iproc < GetNumProcs() ; ++iproc) {
00864 
00865         if (GetMyPID() == iproc) {
00866 
00867           for (int i = 0 ; i < GetMyLength() ; ++i) {
00868             os.width(10);
00869             os << GetMyPID();
00870             os.width(20);
00871             os << i;
00872             os.width(20);
00873             os << GetVectorSpace()(i);
00874             for (int v = 0 ; v < GetNumVectors() ; ++v) {
00875               os.width(20);
00876               os << (*this)(i,v);
00877             }
00878             os << endl;
00879           }
00880         }
00881 
00882         Barrier();
00883       }
00884       if (GetMyPID() == 0)
00885         os << endl;
00886     }
00887 
00888     StackPop();
00889     UpdateTime();
00890 
00891     return(os);
00892   }
00893   // @}
00894 
00895   bool IsAlias(const MultiVector& rhs) const
00896   {
00897     if (rhs.GetValues(0) == GetValues(0))
00898       return(true);
00899     else
00900       return(false);
00901   }
00902 
00903 private:
00904 
00906   void SetRCPLength(const int NumVectors)
00907   {
00908     RCPValues_.resize(NumVectors);
00909   }
00910 
00912   inline void Initialize()
00913   {
00914     for (int v = 0 ; v < GetNumVectors() ; ++v)
00915       RCPValues_[v] = Teuchos::null;
00916   }
00917 
00919   void CheckSpaces(const MultiVector rhs)  const
00920   {
00921     if (rhs.GetVectorSpace() != GetVectorSpace()) {
00922       ML_THROW("rhs.GetVectorSpace() is not equal to this->GetVectorSpace()", -1);
00923     }
00924 
00925     if (IsAlias(rhs))
00926       ML_THROW("updating a vector with its alias...", -1);
00927   }
00928 
00930   inline void CheckEntry(const int i) const
00931   {
00932 #ifdef MLAPI_CHECK
00933     if ((i < 0) || (i >= GetMyLength()))
00934       ML_THROW("Requested component " + GetString(i) +
00935                ", while MyLength() = " + GetString(GetMyLength()), -1);
00936 #endif
00937   }
00938 
00940   inline void CheckVector(const int v) const
00941   {
00942 #ifdef MLAPI_CHECK
00943     if (v < 0 || v >= GetNumVectors())
00944       ML_THROW("Requested vector " + GetString(v) +
00945                ", while NumVectors() = " + GetString(GetNumVectors()), -1);
00946 #endif
00947   }
00948 
00950   inline void CheckNumVectors(const int NumVectors) const
00951   {
00952     if (GetNumVectors() != NumVectors)
00953       ML_THROW("Incompatible number of vectors, " +
00954                GetString(GetNumVectors()) + " vs. " +
00955                GetString(NumVectors), -1);
00956   }
00957 
00959   inline void CheckSingleVector() const
00960   {
00961     if (GetNumVectors() != 1)
00962       ML_THROW("Implicitly requested vector 0, while NumVectors = "
00963                + GetString(GetNumVectors()), -1);
00964   }
00965 
00967   std::vector<Teuchos::RefCountPtr<DoubleVector> > RCPValues_;
00969   Space VectorSpace_;
00971   int NumVectors_;
00972 
00973 }; // MultiVector
00974 
00975 } // namespace MLAPI
00976 
00977 #endif // if ML_MULTIVECTOR_H