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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
43 #include "Teuchos_ScalarTraits.hpp"
44 
49 namespace Anasazi {
50 
52  //
53  //--------Anasazi::EpetraOpMultiVec Implementation-------------------------------
54  //
56 
57  // Construction/Destruction
58 
59  EpetraOpMultiVec::EpetraOpMultiVec(const Teuchos::RCP<Epetra_Operator> &Op, const Epetra_BlockMap& Map_in, const int numvecs)
60  : Epetra_OP( Op )
61  {
62  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector(Map_in, numvecs) );
63  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector(Map_in, numvecs) );
64  }
65 
67  const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride)
68  : Epetra_OP( Op )
69  {
70  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::Copy, Map_in, array, stride, numvecs) );
71  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector(Map_in, numvecs) );
72  }
73 
75  Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index)
76  : Epetra_OP( Op )
77  {
78  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector(CV, P_vec, &(const_cast<std::vector<int> &>(index))[0], index.size()) );
79  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector( P_vec.Map(), index.size()) );
80  }
81 
83  : Epetra_OP( P_vec.Epetra_OP )
84  {
85  Epetra_MV = Teuchos::rcp( new Epetra_MultiVector( *(P_vec.Epetra_MV) ) );
86  Epetra_MV_Temp = Teuchos::rcp( new Epetra_MultiVector( *(P_vec.Epetra_MV_Temp) ) );
87  }
88 
89  //
90  // member functions inherited from Anasazi::MultiVec
91  //
92  //
93  // Simulating a virtual copy constructor. If we could rely on the co-variance
94  // of virtual functions, we could return a pointer to EpetraOpMultiVec
95  // (the derived type) instead of a pointer to the pure virtual base class.
96  //
97 
98  MultiVec<double>* EpetraOpMultiVec::Clone ( const int numvecs ) const
99  {
100  EpetraOpMultiVec *ptr_apv = new EpetraOpMultiVec( Epetra_OP, Epetra_MV->Map(), numvecs );
101  return ptr_apv; // safe upcast.
102  }
103  //
104  // the following is a virtual copy constructor returning
105  // a pointer to the pure virtual class. vector values are
106  // copied.
107  //
108 
110  {
111  EpetraOpMultiVec *ptr_apv = new EpetraOpMultiVec(*this);
112  return ptr_apv; // safe upcast
113  }
114 
115 
116  MultiVec<double>* EpetraOpMultiVec::CloneCopy ( const std::vector<int>& index ) const
117  {
118  EpetraOpMultiVec * ptr_apv = new EpetraOpMultiVec( Epetra_OP, Copy, *Epetra_MV, index);
119  return ptr_apv; // safe upcast.
120  }
121 
122 
123  MultiVec<double>* EpetraOpMultiVec::CloneViewNonConst ( const std::vector<int>& index )
124  {
125  EpetraOpMultiVec * ptr_apv = new EpetraOpMultiVec( Epetra_OP, Epetra_DataAccess::View, *Epetra_MV, index);
126  return ptr_apv; // safe upcast.
127  }
128 
129  const MultiVec<double>* EpetraOpMultiVec::CloneView ( const std::vector<int>& index ) const
130  {
131  EpetraOpMultiVec * ptr_apv = new EpetraOpMultiVec( Epetra_OP, Epetra_DataAccess::View, *Epetra_MV, index);
132  return ptr_apv; // safe upcast.
133  }
134 
135 
136  void EpetraOpMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index )
137  {
138  // this should be revisited to e
139  EpetraOpMultiVec temp_vec(Epetra_OP, Epetra_DataAccess::View, *Epetra_MV, index);
140 
141  int numvecs = index.size();
142  if ( A.GetNumberVecs() != numvecs ) {
143  std::vector<int> index2( numvecs );
144  for(int i=0; i<numvecs; i++)
145  index2[i] = i;
146  EpetraOpMultiVec *tmp_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
147  TEUCHOS_TEST_FOR_EXCEPTION( tmp_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::SetBlocks() cast of MultiVec<double> to EpetraOpMultiVec failed.");
148  EpetraOpMultiVec A_vec(Epetra_OP, Epetra_DataAccess::View, *(tmp_vec->GetEpetraMultiVector()), index2);
149  temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
150  }
151  else {
152  temp_vec.MvAddMv( 1.0, A, 0.0, A );
153  }
154  }
155 
156  //-------------------------------------------------------------
157  //
158  // *this <- alpha * A * B + beta * (*this)
159  //
160  //-------------------------------------------------------------
161 
162  void EpetraOpMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A,
163  const Teuchos::SerialDenseMatrix<int,double>& B, double beta )
164  {
165  Epetra_LocalMap LocalMap(B.numRows(), 0, Epetra_MV->Map().Comm());
166  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
167 
168  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
169  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::SetBlocks() cast of MultiVec<double> to EpetraOpMultiVec failed.");
170 
172  Epetra_MV->Multiply( 'N', 'N', alpha, *(A_vec->GetEpetraMultiVector()), B_Pvec, beta ) != 0,
173  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
174  }
175 
176  //-------------------------------------------------------------
177  //
178  // *this <- alpha * A + beta * B
179  //
180  //-------------------------------------------------------------
181 
182  void EpetraOpMultiVec::MvAddMv ( double alpha , const MultiVec<double>& A,
183  double beta, const MultiVec<double>& B)
184  {
185  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
186  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::MvAddMv() cast of MultiVec<double> to EpetraOpMultiVec failed.");
187  EpetraOpMultiVec *B_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(B));
188  TEUCHOS_TEST_FOR_EXCEPTION( B_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::MvAddMv() cast of MultiVec<double> to EpetraOpMultiVec failed.");
189 
191  Epetra_MV->Update( alpha, *(A_vec->GetEpetraMultiVector()), beta, *(B_vec->GetEpetraMultiVector()), 0.0 ) != 0,
192  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvAddMv() call to Epetra_MultiVec::Update() returned a nonzero value.");
193  }
194 
195  //-------------------------------------------------------------
196  //
197  // dense B <- alpha * A^T * OP * (*this)
198  //
199  //-------------------------------------------------------------
200 
201  void EpetraOpMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
203 #ifdef HAVE_ANASAZI_EXPERIMENTAL
204  , ConjType conj
205 #endif
206  ) const
207  {
208  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
209 
210  if (A_vec) {
211  Epetra_LocalMap LocalMap(B.numRows(), 0, Epetra_MV->Map().Comm());
212  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
213 
214  int info = Epetra_OP->Apply( *Epetra_MV, *Epetra_MV_Temp );
216  "Anasazi::EpetraOpMultiVec::MvTransMv(): Error returned from Epetra_Operator::Apply()" );
217 
219  B_Pvec.Multiply( 'T', 'N', alpha, *(A_vec->GetEpetraMultiVector()), *Epetra_MV_Temp, 0.0 ) != 0,
220  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvTransMv() call to Epetra_MultiVector::Multiply() returned a nonzero value.");
221  }
222  }
223 
224  //-------------------------------------------------------------
225  //
226  // b[i] = A[i]^T * OP * this[i]
227  //
228  //-------------------------------------------------------------
229 
230  void EpetraOpMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double> & b
231 #ifdef HAVE_ANASAZI_EXPERIMENTAL
232  , ConjType conj
233 #endif
234  ) const
235  {
236  EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
237  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraOpMultiVec::MvDot() cast of MultiVec<double> to EpetraOpMultiVec failed.");
238 
239  int info = Epetra_OP->Apply( *Epetra_MV, *Epetra_MV_Temp );
241  "Anasazi::EpetraOpMultiVec::MvDot(): Error returned from Epetra_Operator::Apply()" );
242 
243  if (( (int)b.size() >= A_vec->GetNumberVecs() ) ) {
245  Epetra_MV_Temp->Dot( *(A_vec->GetEpetraMultiVector()), &b[0] ) != 0,
246  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvDot() call to Epetra_MultiVector::Dot() returned an error.");
247  }
248  }
249 
250  //-------------------------------------------------------------
251  //
252  // normvec[i] = || this[i] ||_OP
253  //
254  //-------------------------------------------------------------
255 
256  void EpetraOpMultiVec::MvNorm ( std::vector<double> & normvec ) const
257  {
258  int info = Epetra_OP->Apply( *Epetra_MV, *Epetra_MV_Temp );
260  "Anasazi::EpetraOpMultiVec::MvNorm(): Error returned from Epetra_Operator::Apply()" );
261 
262  if (( (int)normvec.size() >= Epetra_MV->NumVectors() ) ) {
264  Epetra_MV_Temp->Dot( *Epetra_MV, &normvec[0] ) != 0,
265  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvNorm() call to Epetra_MultiVector::Dot() returned an error.");
266  }
267 
268  for (int i=0; i<Epetra_MV->NumVectors(); ++i)
269  normvec[i] = Teuchos::ScalarTraits<double>::squareroot( normvec[i] );
270  }
271 
272  //-------------------------------------------------------------
273  //
274  // this[i] = alpha[i] * this[i]
275  //
276  //-------------------------------------------------------------
277  void EpetraOpMultiVec::MvScale ( const std::vector<double>& alpha )
278  {
279  // Check to make sure the vector is as long as the multivector has columns.
280  int numvecs = this->GetNumberVecs();
281  TEUCHOS_TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument,
282  "Anasazi::EpetraOpMultiVec::MvScale() alpha argument size was inconsistent with number of vectors in mv.");
283 
284  std::vector<int> tmp_index( 1, 0 );
285  for (int i=0; i<numvecs; i++) {
286  Epetra_MultiVector temp_vec(Epetra_DataAccess::View, *Epetra_MV, &tmp_index[0], 1);
288  temp_vec.Scale( alpha[i] ) != 0,
289  EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvScale() call to Epetra_MultiVector::Scale() returned a nonzero value.");
290  tmp_index[0]++;
291  }
292  }
293 
294 } // namespace Anasazi
Declarations of specialized Anasazi multi-vector and operator classes using Epetra_MultiVector and Ep...
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...
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...
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
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...
MultiVec< double > * CloneViewNonConst(const std::vector< int > &index)
Creates a new EpetraOpMultiVec that shares the selected contents of *this.
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.
Copy
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.
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.