Anasazi Version of the Day
Loading...
Searching...
No Matches
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
49namespace 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
66 EpetraOpMultiVec::EpetraOpMultiVec(const Teuchos::RCP<Epetra_Operator> &Op,
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
74 EpetraOpMultiVec::EpetraOpMultiVec(const Teuchos::RCP<Epetra_Operator> &Op,
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
109 MultiVec<double>* EpetraOpMultiVec::CloneCopy() const
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
171 TEUCHOS_TEST_FOR_EXCEPTION(
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
190 TEUCHOS_TEST_FOR_EXCEPTION(
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,
202 Teuchos::SerialDenseMatrix<int,double>& B
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 );
215 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, EpetraSpecializedMultiVecFailure,
216 "Anasazi::EpetraOpMultiVec::MvTransMv(): Error returned from Epetra_Operator::Apply()" );
217
218 TEUCHOS_TEST_FOR_EXCEPTION(
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 );
240 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, EpetraSpecializedMultiVecFailure,
241 "Anasazi::EpetraOpMultiVec::MvDot(): Error returned from Epetra_Operator::Apply()" );
242
243 if (( (int)b.size() >= A_vec->GetNumberVecs() ) ) {
244 TEUCHOS_TEST_FOR_EXCEPTION(
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 );
259 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, EpetraSpecializedMultiVecFailure,
260 "Anasazi::EpetraOpMultiVec::MvNorm(): Error returned from Epetra_Operator::Apply()" );
261
262 if (( (int)normvec.size() >= Epetra_MV->NumVectors() ) ) {
263 TEUCHOS_TEST_FOR_EXCEPTION(
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);
287 TEUCHOS_TEST_FOR_EXCEPTION(
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...
Specialized adapter class for Anasazi::MultiVec that uses Epetra_MultiVector and Epetra_Operator to d...
void MvScale(double alpha)
Scale each element of the vectors in *this with alpha.
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...
MultiVec< double > * Clone(const int numvecs) const
Creates a new empty EpetraOpMultiVec containing numvecs columns.
void MvAddMv(double alpha, const MultiVec< double > &A, double beta, const MultiVec< double > &B)
Replace *this with .
MultiVec< double > * CloneCopy() const
Creates a new EpetraOpMultiVec and copies contents of *this into the new vector (deep copy).
const MultiVec< double > * CloneView(const std::vector< int > &index) const
Creates a new EpetraOpMultiVec that shares the selected contents of *this.
EpetraOpMultiVec(const Teuchos::RCP< Epetra_Operator > &Op, const Epetra_BlockMap &Map_in, const int numvecs)
Basic EpetraOpMultiVec constructor.
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 MvTimesMatAddMv(double alpha, const MultiVec< double > &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta)
Update *this with .
void SetBlock(const MultiVec< double > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
int GetNumberVecs() const
Obtain the vector length of *this.
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 ...
MultiVec< double > * CloneViewNonConst(const std::vector< int > &index)
Creates a new EpetraOpMultiVec that shares the selected contents of *this.
EpetraSpecializedMultiVecFailure is thrown when a return value from an Epetra call on an Epetra_Multi...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ConjType
Enumerated types used to specify conjugation arguments.