Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziEpetraAdapter.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
48namespace Anasazi {
49
51 //
52 //--------Anasazi::EpetraMultiVec Implementation-------------------------------
53 //
55
56 // Construction/Destruction
57
58 EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array,
59 const int numvecs, const int stride)
60 : Epetra_MultiVector(Epetra_DataAccess::Copy, Map_in, array, stride, numvecs)
61 {
62 }
63
64
65 EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs)
66 : Epetra_MultiVector(Map_in, numvecs)
67 {
68 }
69
70
71 EpetraMultiVec::EpetraMultiVec(Epetra_DataAccess CV,
72 const Epetra_MultiVector& P_vec,
73 const std::vector<int>& index )
74 : Epetra_MultiVector(CV, P_vec, &(const_cast<std::vector<int> &>(index))[0], index.size())
75 {
76 }
77
78
79 EpetraMultiVec::EpetraMultiVec(const Epetra_MultiVector& P_vec)
80 : Epetra_MultiVector(P_vec)
81 {
82 }
83
84
85 //
86 // member functions inherited from Anasazi::MultiVec
87 //
88 //
89 // Simulating a virtual copy constructor. If we could rely on the co-variance
90 // of virtual functions, we could return a pointer to EpetraMultiVec
91 // (the derived type) instead of a pointer to the pure virtual base class.
92 //
93
94 MultiVec<double>* EpetraMultiVec::Clone ( const int numvecs ) const
95 {
96 EpetraMultiVec * ptr_apv = new EpetraMultiVec(Map(), numvecs);
97 return ptr_apv; // safe upcast.
98 }
99 //
100 // the following is a virtual copy constructor returning
101 // a pointer to the pure virtual class. vector values are
102 // copied.
103 //
104
106 {
107 EpetraMultiVec *ptr_apv = new EpetraMultiVec(*this);
108 return ptr_apv; // safe upcast
109 }
110
111
112 MultiVec<double>* EpetraMultiVec::CloneCopy ( const std::vector<int>& index ) const
113 {
114 EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::Copy, *this, index);
115 return ptr_apv; // safe upcast.
116 }
117
118
119 MultiVec<double>* EpetraMultiVec::CloneViewNonConst ( const std::vector<int>& index )
120 {
121 EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::View, *this, index);
122 return ptr_apv; // safe upcast.
123 }
124
125 const MultiVec<double>* EpetraMultiVec::CloneView ( const std::vector<int>& index ) const
126 {
127 EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::View, *this, index);
128 return ptr_apv; // safe upcast.
129 }
130
131
132 void EpetraMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index )
133 {
134 // this should be revisited to e
135 EpetraMultiVec temp_vec(Epetra_DataAccess::View, *this, index);
136
137 int numvecs = index.size();
138 if ( A.GetNumberVecs() != numvecs ) {
139 std::vector<int> index2( numvecs );
140 for(int i=0; i<numvecs; i++)
141 index2[i] = i;
142 EpetraMultiVec *tmp_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
143 TEUCHOS_TEST_FOR_EXCEPTION( tmp_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
144 EpetraMultiVec A_vec(Epetra_DataAccess::View, *tmp_vec, index2);
145 temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
146 }
147 else {
148 temp_vec.MvAddMv( 1.0, A, 0.0, A );
149 }
150 }
151
152 //-------------------------------------------------------------
153 //
154 // *this <- alpha * A * B + beta * (*this)
155 //
156 //-------------------------------------------------------------
157
158 void EpetraMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A,
159 const Teuchos::SerialDenseMatrix<int,double>& B, double beta )
160 {
161 Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
162 Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
163
164 EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
165 TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
166
167 TEUCHOS_TEST_FOR_EXCEPTION(
168 Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta ) != 0,
169 EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
170 }
171
172 //-------------------------------------------------------------
173 //
174 // *this <- alpha * A + beta * B
175 //
176 //-------------------------------------------------------------
177
178 void EpetraMultiVec::MvAddMv ( double alpha , const MultiVec<double>& A,
179 double beta, const MultiVec<double>& B)
180 {
181 EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
182 TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
183 EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B));
184 TEUCHOS_TEST_FOR_EXCEPTION( B_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
185
186 TEUCHOS_TEST_FOR_EXCEPTION(
187 Update( alpha, *A_vec, beta, *B_vec, 0.0 ) != 0,
188 EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvAddMv() call to Epetra_MultiVec::Update() returned a nonzero value.");
189 }
190
191 //-------------------------------------------------------------
192 //
193 // dense B <- alpha * A^T * (*this)
194 //
195 //-------------------------------------------------------------
196
197 void EpetraMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
198 Teuchos::SerialDenseMatrix<int,double>& B
199#ifdef HAVE_ANASAZI_EXPERIMENTAL
200 , ConjType conj
201#endif
202 ) const
203 {
204 EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
205
206 if (A_vec) {
207 Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
208 Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
209
210 TEUCHOS_TEST_FOR_EXCEPTION(
211 B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 ) != 0,
212 EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTransMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
213 }
214 }
215
216 //-------------------------------------------------------------
217 //
218 // b[i] = A[i]^T * this[i]
219 //
220 //-------------------------------------------------------------
221
222 void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double> & b
223#ifdef HAVE_ANASAZI_EXPERIMENTAL
224 , ConjType conj
225#endif
226 ) const
227 {
228 EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
229 TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvDot() cast of MultiVec<double> to EpetraMultiVec failed.");
230
231 if (( (int)b.size() >= A_vec->NumVectors() ) ) {
232 TEUCHOS_TEST_FOR_EXCEPTION(
233 this->Dot( *A_vec, &b[0] ) != 0,
234 EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvDot() call to Epetra_MultiVec::Dot() returned a nonzero value.");
235 }
236 }
237
238 //-------------------------------------------------------------
239 //
240 // this[i] = alpha[i] * this[i]
241 //
242 //-------------------------------------------------------------
243 void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
244 {
245 // Check to make sure the vector is as long as the multivector has columns.
246 int numvecs = this->NumVectors();
247 TEUCHOS_TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument,
248 "Anasazi::EpetraMultiVec::MvScale() alpha argument size was inconsistent with number of vectors in mv.");
249
250 std::vector<int> tmp_index( 1, 0 );
251 for (int i=0; i<numvecs; i++) {
252 Epetra_MultiVector temp_vec(Epetra_DataAccess::View, *this, &tmp_index[0], 1);
253 TEUCHOS_TEST_FOR_EXCEPTION(
254 temp_vec.Scale( alpha[i] ) != 0,
255 EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvScale() call to Epetra_MultiVec::Scale() returned a nonzero value.");
256 tmp_index[0]++;
257 }
258 }
259
261 //
262 //--------Anasazi::EpetraOp Implementation-------------------------------------
263 //
265
266 //
267 // AnasaziOperator constructors
268 //
269 EpetraOp::EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op)
270 : Epetra_Op(Op)
271 {
272 }
273
275 {
276 }
277 //
278 // AnasaziOperator applications
279 //
281 MultiVec<double>& Y ) const
282 {
283 //
284 // This standard operator computes Y = A*X
285 //
286 MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
287 Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
288 Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
289
290 TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
291 TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
292
293 int info = Epetra_Op->Apply( *vec_X, *vec_Y );
294 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
295 "Anasazi::EpetraOp::Apply(): Error returned from Epetra_Operator::Apply()" );
296 }
297
299 //
300 //--------Anasazi::EpetraGenOp Implementation----------------------------------
301 //
303
304 //
305 // AnasaziOperator constructors
306 //
307
308 EpetraGenOp::EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp,
309 const Teuchos::RCP<Epetra_Operator> &MOp,
310 bool isAInverse_)
311 : isAInverse( isAInverse_ ), Epetra_AOp(AOp), Epetra_MOp(MOp)
312 {
313 }
314
316 {
317 }
318 //
319 // AnasaziOperator applications
320 //
322 {
323 //
324 // This generalized operator computes Y = A^{-1}*M*X
325 //
326 int info=0;
327 MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
328 Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
329 Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
330 Epetra_MultiVector temp_Y(*vec_Y);
331
332 TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
333 TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
334 //
335 // Need to cast away constness because the member function Apply is not declared const.
336 // Change the transpose setting for the operator if necessary and change it back when done.
337 //
338 // Apply M
339 info = Epetra_MOp->Apply( *vec_X, temp_Y );
340 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
341 "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
342 // Apply A or A^{-1}
343 if (isAInverse) {
344 info = Epetra_AOp->ApplyInverse( temp_Y, *vec_Y );
345 }
346 else {
347 info = Epetra_AOp->Apply( temp_Y, *vec_Y );
348 }
349 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
350 "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
351 }
352
353 int EpetraGenOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
354 {
355 //
356 // This generalized operator computes Y = A^{-1}*M*X
357 //
358 int info=0;
359 Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors());
360
361 // Apply M
362 info = Epetra_MOp->Apply( X, temp_Y );
363 if (info!=0) return info;
364
365 // Apply A or A^{-1}
366 if (isAInverse)
367 info = Epetra_AOp->ApplyInverse( temp_Y, Y );
368 else
369 info = Epetra_AOp->Apply( temp_Y, Y );
370
371 return info;
372 }
373
374 int EpetraGenOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
375 {
376 //
377 // This generalized operator computes Y = M^{-1}*A*X
378 //
379 int info=0;
380 Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors());
381
382 // Apply A or A^{-1}
383 if (isAInverse)
384 info = Epetra_AOp->Apply( X, temp_Y );
385 else
386 info = Epetra_AOp->ApplyInverse( X, temp_Y );
387 if (info!=0) return info;
388
389 // Apply M^{-1}
390 info = Epetra_MOp->ApplyInverse( temp_Y, Y );
391
392 return info;
393 }
394
396 //
397 //--------Anasazi::EpetraSymOp Implementation----------------------------------
398 //
400
401 //
402 // AnasaziOperator constructors
403 //
404 EpetraSymOp::EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op,
405 bool isTrans)
406 : Epetra_Op(Op), isTrans_(isTrans)
407 {
408 }
409
411 {
412 }
413 //
414 // AnasaziOperator applications
415 //
417 MultiVec<double>& Y ) const
418 {
419 int info=0;
420 MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
421 Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
422 Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
423 Epetra_MultiVector* temp_vec = new Epetra_MultiVector(
424 (isTrans_) ? Epetra_Op->OperatorDomainMap()
425 : Epetra_Op->OperatorRangeMap(),
426 vec_X->NumVectors() );
427
428 TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
429 TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
430 TEUCHOS_TEST_FOR_EXCEPTION( temp_vec==NULL, std::invalid_argument, "Anasazi::EpetraSymOp::Apply() allocation Epetra_MultiVector failed.");
431 //
432 // Need to cast away constness because the member function Apply
433 // is not declared const.
434 //
435 // Transpose the operator (if isTrans_ = true)
436 if (isTrans_) {
437 info = Epetra_Op->SetUseTranspose( isTrans_ );
438 if (info != 0) {
439 delete temp_vec;
440 TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
441 "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
442 }
443 }
444 //
445 // Compute A*X or A'*X
446 //
447 info=Epetra_Op->Apply( *vec_X, *temp_vec );
448 if (info!=0) {
449 delete temp_vec;
450 TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
451 "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
452 }
453 //
454 // Transpose/Un-transpose the operator based on value of isTrans_
455 info=Epetra_Op->SetUseTranspose( !isTrans_ );
456 if (info!=0) {
457 delete temp_vec;
458 TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
459 "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
460 }
461
462 // Compute A^T*(A*X) or A*A^T
463 info=Epetra_Op->Apply( *temp_vec, *vec_Y );
464 if (info!=0) {
465 delete temp_vec;
466 TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
467 "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
468 }
469
470 // Un-transpose the operator
471 info=Epetra_Op->SetUseTranspose( false );
472 delete temp_vec;
473 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
474 "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
475 }
476
477 int EpetraSymOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
478 {
479 int info=0;
480 Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors());
481 //
482 // This generalized operator computes Y = A^T*A*X or Y = A*A^T*X
483 //
484 // Transpose the operator (if isTrans_ = true)
485 if (isTrans_) {
486 info=Epetra_Op->SetUseTranspose( isTrans_ );
487 if (info!=0) { return info; }
488 }
489 //
490 // Compute A*X or A^T*X
491 //
492 info=Epetra_Op->Apply( X, temp_vec );
493 if (info!=0) { return info; }
494 //
495 // Transpose/Un-transpose the operator based on value of isTrans_
496 info=Epetra_Op->SetUseTranspose( !isTrans_ );
497 if (info!=0) { return info; }
498
499 // Compute A^T*(A*X) or A*A^T
500 info=Epetra_Op->Apply( temp_vec, Y );
501 if (info!=0) { return info; }
502
503 // Un-transpose the operator
504 info=Epetra_Op->SetUseTranspose( false );
505 return info;
506 }
507
508 int EpetraSymOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
509 {
510 int info=0;
511 Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors());
512 //
513 // This generalized operator computes Y = (A^T*A)^{-1}*X or Y = (A*A^T)^{-1}*X
514 //
515 // Transpose the operator (if isTrans_ = true)
516 if (!isTrans_) {
517 info=Epetra_Op->SetUseTranspose( !isTrans_ );
518 if (info!=0) { return info; }
519 }
520 //
521 // Compute A^{-1}*X or A^{-T}*X
522 //
523 info=Epetra_Op->ApplyInverse( X, temp_vec );
524 if (info!=0) { return info; }
525 //
526 // Transpose/Un-transpose the operator based on value of isTrans_
527 info=Epetra_Op->SetUseTranspose( isTrans_ );
528 if (info!=0) { return info; }
529
530 // Compute A^{-T}*(A^{-1}*X) or A^{-1}*A^{-T}
531 info=Epetra_Op->Apply( temp_vec, Y );
532 if (info!=0) { return info; }
533
534 // Un-transpose the operator
535 info=Epetra_Op->SetUseTranspose( false );
536 return info;
537 }
538
540 //
541 //--------Anasazi::EpetraSymMVOp Implementation--------------------------------
542 //
544
545 //
546 // Anasazi::Operator constructors
547 //
548 EpetraSymMVOp::EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
549 bool isTrans)
550 : Epetra_MV(MV), isTrans_(isTrans)
551 {
552 if (isTrans)
553 MV_localmap = Teuchos::rcp( new Epetra_LocalMap( Epetra_MV->NumVectors(), 0, Epetra_MV->Map().Comm() ) );
554 else
555 MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
556 }
557
558 //
559 // AnasaziOperator applications
560 //
562 {
563 int info=0;
564 MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
565 Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
566 Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
567
568 if (isTrans_) {
569
570 Epetra_MultiVector temp_vec( *MV_localmap, temp_X.GetNumberVecs() );
571
572 /* A'*X */
573 info = temp_vec.Multiply( 'T', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
574 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
575 "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
576
577 /* A*(A'*X) */
578 info = vec_Y->Multiply( 'N', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
579 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
580 "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
581 }
582 else {
583
584 Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
585
586 /* A*X */
587 info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
588 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
589 "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
590
591 /* A'*(A*X) */
592 info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
593 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
594 "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
595 }
596 }
597
599 //
600 //--------Anasazi::EpetraWSymMVOp Implementation--------------------------------
601 //
603
604 //
605 // Anasazi::Operator constructors
606 //
607 EpetraWSymMVOp::EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
608 const Teuchos::RCP<Epetra_Operator> &OP )
609 : Epetra_MV(MV), Epetra_OP(OP)
610 {
611 MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
612 Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) );
613 Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
614 }
615
616 //
617 // AnasaziOperator applications
618 //
620 {
621 int info=0;
622 MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
623 Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
624 Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
625
626 Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
627
628 /* WA*X */
629 info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
630 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
631 "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
632
633 /* A'*(WA*X) */
634 info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
635 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
636 "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
637 }
638
640 //
641 //--------Anasazi::EpetraW2SymMVOp Implementation--------------------------------
642 //
644
645 //
646 // Anasazi::Operator constructors
647 //
648 EpetraW2SymMVOp::EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
649 const Teuchos::RCP<Epetra_Operator> &OP )
650 : Epetra_MV(MV), Epetra_OP(OP)
651 {
652 MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
653 Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) );
654 Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
655 }
656
657 //
658 // AnasaziOperator applications
659 //
661 {
662 int info=0;
663 MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
664 Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
665 Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
666
667 Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
668
669 /* WA*X */
670 info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
671 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
672 "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
673
674 /* (WA)'*(WA*X) */
675 info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_WMV, temp_vec, 0.0 );
676 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
677 "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
678
679 }
680} // end namespace Anasazi
Declarations of Anasazi multi-vector and operator classes using Epetra_MultiVector and Epetra_Operato...
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method [inherited from Anasazi::Operator class].
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Apply inverse method [inherited from Epetra_Operator class].
EpetraGenOp(const Teuchos::RCP< Epetra_Operator > &AOp, const Teuchos::RCP< Epetra_Operator > &MOp, bool isAInverse=true)
Basic constructor for applying operator [default] or .
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
EpetraMultiVecAccessor is an interfaceto allow any Anasazi::MultiVec implementation that is based on ...
EpetraMultiVecFailure is thrown when a return value from an Epetra call on an Epetra_MultiVector is n...
Basic adapter class for Anasazi::MultiVec that uses Epetra_MultiVector.
MultiVec< double > * Clone(const int numvecs) const
Creates a new empty EpetraMultiVec containing numvecs columns.
void MvScale(double alpha)
Scale each element of the vectors in *this with alpha.
void MvAddMv(double alpha, const MultiVec< double > &A, double beta, const MultiVec< double > &B)
Replace *this with .
const MultiVec< double > * CloneView(const std::vector< int > &index) const
Creates a new EpetraMultiVec that shares the selected contents of *this.
void SetBlock(const MultiVec< double > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
EpetraMultiVec(const Epetra_BlockMap &Map_in, const int numvecs)
Basic EpetraMultiVec 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 .
MultiVec< double > * CloneViewNonConst(const std::vector< int > &index)
Creates a new EpetraMultiVec that shares the selected contents of *this.
MultiVec< double > * CloneCopy() const
Creates a new EpetraMultiVec and copies contents of *this into the new vector (deep copy).
void MvTimesMatAddMv(double alpha, const MultiVec< double > &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta)
Update *this with .
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 Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
This method takes the Anasazi::MultiVec X and applies the operator to it resulting in the Anasazi::Mu...
EpetraOp(const Teuchos::RCP< Epetra_Operator > &Op)
Basic constructor. Accepts reference-counted pointer to an Epetra_Operator.
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
EpetraSymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, bool isTrans=false)
Basic constructor for applying operator [default] or .
EpetraSymOp(const Teuchos::RCP< Epetra_Operator > &Op, bool isTrans=false)
Basic constructor for applying operator [default] or .
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Apply inverse method [inherited from Epetra_Operator class].
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method [inherited from Anasazi::Operator class].
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
EpetraW2SymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, const Teuchos::RCP< Epetra_Operator > &OP)
Basic constructor for applying operator .
EpetraWSymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, const Teuchos::RCP< Epetra_Operator > &OP)
Basic constructor for applying operator .
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
Interface for multivectors used by Anasazi's linear solvers.
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
Exceptions thrown to signal error in operator application.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ConjType
Enumerated types used to specify conjugation arguments.