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
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
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...
Interface for multivectors used by Anasazi's linear solvers.
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ConjType
Enumerated types used to specify conjugation arguments.