Anasazi  Version of the Day
AnasaziSpecializedEpetraAdapter.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
30 #include "Teuchos_ScalarTraits.hpp"
31 
36 namespace Anasazi {
37 
39  //
40  //--------Anasazi::EpetraOpMultiVec Implementation-------------------------------
41  //
43 
44  // Construction/Destruction
45 
47  : Epetra_OP( Op )
48  {
49  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector(Map_in, numvecs) );
50  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector(Map_in, numvecs) );
51  }
52 
54  const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride)
55  : Epetra_OP( Op )
56  {
57  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::Copy, Map_in, array, stride, numvecs) );
58  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector(Map_in, numvecs) );
59  }
60 
62  Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index)
63  : Epetra_OP( Op )
64  {
65  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector(CV, P_vec, &(const_cast<std::vector<int> &>(index))[0], index.size()) );
66  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector( P_vec.Map(), index.size()) );
67  }
68 
70  : Epetra_OP( P_vec.Epetra_OP )
71  {
72  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector( *(P_vec.Epetra_MV) ) );
73  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector( *(P_vec.Epetra_MV_Temp) ) );
74  }
75 
76  //
77  // member functions inherited from Anasazi::MultiVec
78  //
79  //
80  // Simulating a virtual copy constructor. If we could rely on the co-variance
81  // of virtual functions, we could return a pointer to EpetraOpMultiVec
82  // (the derived type) instead of a pointer to the pure virtual base class.
83  //
84 
85  MultiVec<double>* EpetraOpMultiVec::Clone ( const int numvecs ) const
86  {
87  EpetraOpMultiVec *ptr_apv = new EpetraOpMultiVec( Epetra_OP, Epetra_MV->Map(), numvecs );
88  return ptr_apv; // safe upcast.
89  }
90  //
91  // the following is a virtual copy constructor returning
92  // a pointer to the pure virtual class. vector values are
93  // copied.
94  //
95 
97  {
98  EpetraOpMultiVec *ptr_apv = new EpetraOpMultiVec(*this);
99  return ptr_apv; // safe upcast
100  }
101 
102 
103  MultiVec<double>* EpetraOpMultiVec::CloneCopy ( const std::vector<int>& index ) const
104  {
105  EpetraOpMultiVec * ptr_apv = new EpetraOpMultiVec( Epetra_OP, Copy, *Epetra_MV, index);
106  return ptr_apv; // safe upcast.
107  }
108 
109 
110  MultiVec<double>* EpetraOpMultiVec::CloneViewNonConst ( const std::vector<int>& index )
111  {
112  EpetraOpMultiVec * ptr_apv = new EpetraOpMultiVec( Epetra_OP, Epetra_DataAccess::View, *Epetra_MV, index);
113  return ptr_apv; // safe upcast.
114  }
115 
116  const MultiVec<double>* EpetraOpMultiVec::CloneView ( const std::vector<int>& index ) const
117  {
118  EpetraOpMultiVec * ptr_apv = new EpetraOpMultiVec( Epetra_OP, Epetra_DataAccess::View, *Epetra_MV, index);
119  return ptr_apv; // safe upcast.
120  }
121 
122 
123  void EpetraOpMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index )
124  {
125  // this should be revisited to e
126  EpetraOpMultiVec temp_vec(Epetra_OP, Epetra_DataAccess::View, *Epetra_MV, index);
127 
128  int numvecs = index.size();
129  if ( A.GetNumberVecs() != numvecs ) {
130  std::vector<int> index2( numvecs );
131  for(int i=0; i<numvecs; i++)
132  index2[i] = i;
133  EpetraOpMultiVec *tmp_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
134  TEUCHOS_TEST_FOR_EXCEPTION( tmp_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::SetBlocks() cast of MultiVec<double> to EpetraOpMultiVec failed.");
135  EpetraOpMultiVec A_vec(Epetra_OP, Epetra_DataAccess::View, *(tmp_vec->GetEpetraMultiVector()), index2);
136  temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
137  }
138  else {
139  temp_vec.MvAddMv( 1.0, A, 0.0, A );
140  }
141  }
142 
143  //-------------------------------------------------------------
144  //
145  // *this <- alpha * A * B + beta * (*this)
146  //
147  //-------------------------------------------------------------
148 
149  void EpetraOpMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A,
150  const Teuchos::SerialDenseMatrix<int,double>& B, double beta )
151  {
152  Epetra_LocalMap LocalMap(B.numRows(), 0, Epetra_MV->Map().Comm());
153  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
154 
155  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
156  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::SetBlocks() cast of MultiVec<double> to EpetraOpMultiVec failed.");
157 
159  Epetra_MV->Multiply( 'N', 'N', alpha, *(A_vec->GetEpetraMultiVector()), B_Pvec, beta ) != 0,
160  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
161  }
162 
163  //-------------------------------------------------------------
164  //
165  // *this <- alpha * A + beta * B
166  //
167  //-------------------------------------------------------------
168 
169  void EpetraOpMultiVec::MvAddMv ( double alpha , const MultiVec<double>& A,
170  double beta, const MultiVec<double>& B)
171  {
172  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
173  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::MvAddMv() cast of MultiVec<double> to EpetraOpMultiVec failed.");
174  EpetraOpMultiVec *B_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(B));
175  TEUCHOS_TEST_FOR_EXCEPTION( B_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::MvAddMv() cast of MultiVec<double> to EpetraOpMultiVec failed.");
176 
178  Epetra_MV->Update( alpha, *(A_vec->GetEpetraMultiVector()), beta, *(B_vec->GetEpetraMultiVector()), 0.0 ) != 0,
179  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvAddMv() call to Epetra_MultiVec::Update() returned a nonzero value.");
180  }
181 
182  //-------------------------------------------------------------
183  //
184  // dense B <- alpha * A^T * OP * (*this)
185  //
186  //-------------------------------------------------------------
187 
188  void EpetraOpMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
190 #ifdef HAVE_ANASAZI_EXPERIMENTAL
191  , ConjType conj
192 #endif
193  ) const
194  {
195  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
196 
197  if (A_vec) {
198  Epetra_LocalMap LocalMap(B.numRows(), 0, Epetra_MV->Map().Comm());
199  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
200 
201  int info = Epetra_OP->Apply( *Epetra_MV, *Epetra_MV_Temp );
203  "Anasazi::EpetraOpMultiVec::MvTransMv(): Error returned from Epetra_Operator::Apply()" );
204 
206  B_Pvec.Multiply( 'T', 'N', alpha, *(A_vec->GetEpetraMultiVector()), *Epetra_MV_Temp, 0.0 ) != 0,
207  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvTransMv() call to Epetra_MultiVector::Multiply() returned a nonzero value.");
208  }
209  }
210 
211  //-------------------------------------------------------------
212  //
213  // b[i] = A[i]^T * OP * this[i]
214  //
215  //-------------------------------------------------------------
216 
217  void EpetraOpMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double> & b
218 #ifdef HAVE_ANASAZI_EXPERIMENTAL
219  , ConjType conj
220 #endif
221  ) const
222  {
223  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
224  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::MvDot() cast of MultiVec<double> to EpetraOpMultiVec failed.");
225 
226  int info = Epetra_OP->Apply( *Epetra_MV, *Epetra_MV_Temp );
228  "Anasazi::EpetraOpMultiVec::MvDot(): Error returned from Epetra_Operator::Apply()" );
229 
230  if (( (int)b.size() >= A_vec->GetNumberVecs() ) ) {
232  Epetra_MV_Temp->Dot( *(A_vec->GetEpetraMultiVector()), &b[0] ) != 0,
233  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvDot() call to Epetra_MultiVector::Dot() returned an error.");
234  }
235  }
236 
237  //-------------------------------------------------------------
238  //
239  // normvec[i] = || this[i] ||_OP
240  //
241  //-------------------------------------------------------------
242 
243  void EpetraOpMultiVec::MvNorm ( std::vector<double> & normvec ) const
244  {
245  int info = Epetra_OP->Apply( *Epetra_MV, *Epetra_MV_Temp );
247  "Anasazi::EpetraOpMultiVec::MvNorm(): Error returned from Epetra_Operator::Apply()" );
248 
249  if (( (int)normvec.size() >= Epetra_MV->NumVectors() ) ) {
251  Epetra_MV_Temp->Dot( *Epetra_MV, &normvec[0] ) != 0,
252  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvNorm() call to Epetra_MultiVector::Dot() returned an error.");
253  }
254 
255  for (int i=0; i<Epetra_MV->NumVectors(); ++i)
256  normvec[i] = Teuchos::ScalarTraits<double>::squareroot( normvec[i] );
257  }
258 
259  //-------------------------------------------------------------
260  //
261  // this[i] = alpha[i] * this[i]
262  //
263  //-------------------------------------------------------------
264  void EpetraOpMultiVec::MvScale ( const std::vector<double>& alpha )
265  {
266  // Check to make sure the vector is as long as the multivector has columns.
267  int numvecs = this->GetNumberVecs();
268  TEUCHOS_TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument,
269  "Anasazi::EpetraOpMultiVec::MvScale() alpha argument size was inconsistent with number of vectors in mv.");
270 
271  std::vector<int> tmp_index( 1, 0 );
272  for (int i=0; i<numvecs; i++) {
273  Epetra_MultiVector temp_vec(Epetra_DataAccess::View, *Epetra_MV, &tmp_index[0], 1);
275  temp_vec.Scale( alpha[i] ) != 0,
276  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvScale() call to Epetra_MultiVector::Scale() returned a nonzero value.");
277  tmp_index[0]++;
278  }
279  }
280 
281 } // namespace Anasazi
const Epetra_BlockMap & Map() const
Declarations of specialized Anasazi multi-vector and operator classes using Epetra_MultiVector and Ep...
int Dot(const Epetra_MultiVector &A, double *Result) const
ScalarType * values() const
void MvNorm(std::vector< double > &normvec) const
Compute the 2-norm of each individual vector of *this. Upon return, normvec[i] holds the 2-norm of t...
Copy
EpetraSpecializedMultiVecFailure is thrown when a return value from an Epetra call on an Epetra_Multi...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MultiVec< double > * Clone(const int numvecs) const
Creates a new empty EpetraOpMultiVec containing numvecs columns.
EpetraOpMultiVec(const Teuchos::RCP< Epetra_Operator > &Op, const Epetra_BlockMap &Map_in, const int numvecs)
Basic EpetraOpMultiVec constructor.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
int Scale(double ScalarValue)
ConjType
Enumerated types used to specify conjugation arguments.
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
OrdinalType numRows() const
int NumVectors() const
void MvDot(const MultiVec< double > &A, std::vector< double > &b) const
Compute a vector b where the components are the individual dot-products, i.e. where A[i] is the i-th...
void MvScale(double alpha)
Scale each element of the vectors in *this with alpha.
Specialized adapter class for Anasazi::MultiVec that uses Epetra_MultiVector and Epetra_Operator to d...
const Epetra_Comm & Comm() const
MultiVec< double > * CloneViewNonConst(const std::vector< int > &index)
Creates a new EpetraOpMultiVec that shares the selected contents of *this.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
void MvTransMv(double alpha, const MultiVec< double > &A, Teuchos::SerialDenseMatrix< int, double > &B) const
Compute a dense matrix B through the matrix-matrix multiply .
void SetBlock(const MultiVec< double > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
void MvTimesMatAddMv(double alpha, const MultiVec< double > &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta)
Update *this with .
OrdinalType stride() const
const MultiVec< double > * CloneView(const std::vector< int > &index) const
Creates a new EpetraOpMultiVec that shares the selected contents of *this.
Epetra_DataAccess
OrdinalType numCols() const
MultiVec< double > * CloneCopy() const
Creates a new EpetraOpMultiVec and copies contents of *this into the new vector (deep copy)...
void MvAddMv(double alpha, const MultiVec< double > &A, double beta, const MultiVec< double > &B)
Replace *this with .
int GetNumberVecs() const
Obtain the vector length of *this.