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 
46  EpetraOpMultiVec::EpetraOpMultiVec(const Teuchos::RCP<Epetra_Operator> &Op, const Epetra_BlockMap& Map_in, const int numvecs)
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 
53  EpetraOpMultiVec::EpetraOpMultiVec(const Teuchos::RCP<Epetra_Operator> &Op,
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 
61  EpetraOpMultiVec::EpetraOpMultiVec(const Teuchos::RCP<Epetra_Operator> &Op,
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 
158  TEUCHOS_TEST_FOR_EXCEPTION(
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 
177  TEUCHOS_TEST_FOR_EXCEPTION(
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,
189  Teuchos::SerialDenseMatrix<int,double>& B
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 );
202  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, EpetraSpecializedMultiVecFailure,
203  "Anasazi::EpetraOpMultiVec::MvTransMv(): Error returned from Epetra_Operator::Apply()" );
204 
205  TEUCHOS_TEST_FOR_EXCEPTION(
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 );
227  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, EpetraSpecializedMultiVecFailure,
228  "Anasazi::EpetraOpMultiVec::MvDot(): Error returned from Epetra_Operator::Apply()" );
229 
230  if (( (int)b.size() >= A_vec->GetNumberVecs() ) ) {
231  TEUCHOS_TEST_FOR_EXCEPTION(
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 );
246  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, EpetraSpecializedMultiVecFailure,
247  "Anasazi::EpetraOpMultiVec::MvNorm(): Error returned from Epetra_Operator::Apply()" );
248 
249  if (( (int)normvec.size() >= Epetra_MV->NumVectors() ) ) {
250  TEUCHOS_TEST_FOR_EXCEPTION(
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);
274  TEUCHOS_TEST_FOR_EXCEPTION(
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
Declarations of specialized Anasazi multi-vector and operator classes using Epetra_MultiVector and Ep...
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...
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.
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.
void MvTimesMatAddMv(double alpha, const MultiVec< double > &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta)
Update *this with .
const MultiVec< double > * CloneView(const std::vector< int > &index) const
Creates a new EpetraOpMultiVec that shares the selected contents of *this.
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.