Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziThyraAdapter.hpp
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
47#ifndef ANASAZI_THYRA_ADAPTER_HPP
48#define ANASAZI_THYRA_ADAPTER_HPP
49
52#include "AnasaziConfigDefs.hpp"
53
54#include <Thyra_DetachedMultiVectorView.hpp>
55#include <Thyra_MultiVectorBase.hpp>
56#include <Thyra_MultiVectorStdOps.hpp>
57#include <Thyra_VectorStdOps.hpp>
58
59#include <Teuchos_Ptr.hpp>
60#include <Teuchos_ArrayRCP.hpp>
61#include <Teuchos_ArrayView.hpp>
62
63namespace Anasazi {
64
66 //
67 // Implementation of the Anasazi::MultiVecTraits for Thyra::MultiVectorBase
68 //
70
79 template<class ScalarType>
80 class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
81 {
82 private:
83 typedef Thyra::MultiVectorBase<ScalarType> TMVB;
84 typedef Teuchos::ScalarTraits<ScalarType> ST;
85 typedef typename ST::magnitudeType magType;
86
87 public:
88
91
96 static Teuchos::RCP<TMVB> Clone( const TMVB & mv, const int numvecs )
97 {
98 Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs );
99 return c;
100 }
101
106 static Teuchos::RCP<TMVB>
107 CloneCopy (const TMVB& mv)
108 {
109 const int numvecs = mv.domain()->dim();
110 // create the new multivector
111 Teuchos::RCP< TMVB > cc = Thyra::createMembers (mv.range(), numvecs);
112 // copy the data from the source multivector to the new multivector
113 Thyra::assign (Teuchos::outArg (*cc), mv);
114 return cc;
115 }
116
122 static Teuchos::RCP< TMVB > CloneCopy( const TMVB & mv, const std::vector<int>& index )
123 {
124 const int numvecs = index.size();
125 // create the new multivector
126 Teuchos::RCP<TMVB > cc = Thyra::createMembers (mv.range(), numvecs);
127 // create a view to the relevant part of the source multivector
128 Teuchos::RCP<const TMVB > view = mv.subView ( Teuchos::arrayViewFromVector( index ) );
129 // copy the data from the relevant view to the new multivector
130 Thyra::assign (Teuchos::outArg (*cc), *view);
131 return cc;
132 }
133
134 static Teuchos::RCP<TMVB>
135 CloneCopy (const TMVB& mv, const Teuchos::Range1D& index)
136 {
137 const int numVecs = index.size();
138 // Create the new multivector
139 Teuchos::RCP<TMVB> cc = Thyra::createMembers (mv.range(), numVecs);
140 // Create a view to the relevant part of the source multivector
141 Teuchos::RCP<const TMVB> view = mv.subView (index);
142 // Copy the data from the view to the new multivector.
143 Thyra::assign (Teuchos::outArg (*cc), *view);
144 return cc;
145 }
146
152 static Teuchos::RCP< TMVB > CloneViewNonConst( TMVB & mv, const std::vector<int>& index )
153 {
154 int numvecs = index.size();
155
156 // We do not assume that the indices are sorted, nor do we check that
157 // index.size() > 0. This code is fail-safe, in the sense that a zero
158 // length index vector will pass the error on the Thyra.
159
160 // Thyra has two ways to create an indexed View:
161 // * contiguous (via a range of columns)
162 // * indexed (via a vector of column indices)
163 // The former is significantly more efficient than the latter, in terms of
164 // computations performed with/against the created view.
165 // We will therefore check to see if the given indices are contiguous, and
166 // if so, we will use the contiguous view creation method.
167
168 int lb = index[0];
169 bool contig = true;
170 for (int i=0; i<numvecs; i++) {
171 if (lb+i != index[i]) contig = false;
172 }
173
174 Teuchos::RCP< TMVB > cc;
175 if (contig) {
176 const Thyra::Range1D rng(lb,lb+numvecs-1);
177 // create a contiguous view to the relevant part of the source multivector
178 cc = mv.subView(rng);
179 }
180 else {
181 // create an indexed view to the relevant part of the source multivector
182 cc = mv.subView( Teuchos::arrayViewFromVector( index ) );
183 }
184 return cc;
185 }
186
187 static Teuchos::RCP<TMVB>
188 CloneViewNonConst (TMVB& mv, const Teuchos::Range1D& index)
189 {
190 // We let Thyra be responsible for checking that the index range
191 // is nonempty.
192 //
193 // Create and return a contiguous view to the relevant part of
194 // the source multivector.
195 return mv.subView (index);
196 }
197
203 static Teuchos::RCP<const TMVB > CloneView( const TMVB & mv, const std::vector<int>& index )
204 {
205 int numvecs = index.size();
206
207 // We do not assume that the indices are sorted, nor do we check that
208 // index.size() > 0. This code is fail-safe, in the sense that a zero
209 // length index vector will pass the error on the Thyra.
210
211 // Thyra has two ways to create an indexed View:
212 // * contiguous (via a range of columns)
213 // * indexed (via a vector of column indices)
214 // The former is significantly more efficient than the latter, in terms of
215 // computations performed with/against the created view.
216 // We will therefore check to see if the given indices are contiguous, and
217 // if so, we will use the contiguous view creation method.
218
219 int lb = index[0];
220 bool contig = true;
221 for (int i=0; i<numvecs; i++) {
222 if (lb+i != index[i]) contig = false;
223 }
224
225 Teuchos::RCP< const TMVB > cc;
226 if (contig) {
227 const Thyra::Range1D rng(lb,lb+numvecs-1);
228 // create a contiguous view to the relevant part of the source multivector
229 cc = mv.subView(rng);
230 }
231 else {
232 // create an indexed view to the relevant part of the source multivector
233 cc = mv.subView(Teuchos::arrayViewFromVector( index ) );
234 }
235 return cc;
236 }
237
238 static Teuchos::RCP<const TMVB>
239 CloneView (const TMVB& mv, const Teuchos::Range1D& index)
240 {
241 // We let Thyra be responsible for checking that the index range
242 // is nonempty.
243 //
244 // Create and return a contiguous view to the relevant part of
245 // the source multivector.
246 return mv.subView (index);
247 }
248
250
253
255 static ptrdiff_t GetGlobalLength( const TMVB & mv )
256 { return mv.range()->dim(); }
257
259 static int GetNumberVecs( const TMVB & mv )
260 { return mv.domain()->dim(); }
261
263
266
269 static void MvTimesMatAddMv( const ScalarType alpha, const TMVB & A,
270 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
271 const ScalarType beta, TMVB & mv )
272 {
273 int m = B.numRows();
274 int n = B.numCols();
275 // Create a view of the B object!
276 Teuchos::RCP< const TMVB >
277 B_thyra = Thyra::createMembersView(
278 A.domain(),
279 RTOpPack::ConstSubMultiVectorView<ScalarType>(
280 0, m, 0, n,
281 Teuchos::arcpFromArrayView(Teuchos::arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
282 )
283 );
284 // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
285 A.apply(Thyra::NOTRANS,*B_thyra,Teuchos::outArg (mv),alpha,beta);
286 }
287
290 static void MvAddMv( const ScalarType alpha, const TMVB & A,
291 const ScalarType beta, const TMVB & B, TMVB & mv )
292 {
293 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::inoutArg;
294
295 Thyra::linear_combination<ScalarType>(
296 tuple(alpha, beta)(), tuple(ptrInArg(A), ptrInArg(B))(), ST::zero(), inoutArg(mv));
297 }
298
301 static void MvTransMv( const ScalarType alpha, const TMVB & A, const TMVB & mv,
302 Teuchos::SerialDenseMatrix<int,ScalarType>& B )
303 {
304 // Create a multivector to hold the result (m by n)
305 int m = A.domain()->dim();
306 int n = mv.domain()->dim();
307 // Create a view of the B object!
308 Teuchos::RCP< TMVB >
309 B_thyra = Thyra::createMembersView(
310 A.domain(),
311 RTOpPack::SubMultiVectorView<ScalarType>(
312 0, m, 0, n,
313 Teuchos::arcpFromArrayView(Teuchos::arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
314 )
315 );
316 A.apply(Thyra::CONJTRANS,mv,B_thyra.ptr(),alpha,Teuchos::ScalarTraits<ScalarType>::zero());
317 }
318
322 static void MvDot( const TMVB & mv, const TMVB & A, std::vector<ScalarType> &b )
323 { Thyra::dots(mv,A,Teuchos::arrayViewFromVector( b )); }
324
327 static void
328 MvScale (TMVB& mv,
329 const ScalarType alpha)
330 {
331 Thyra::scale (alpha, Teuchos::inOutArg (mv));
332 }
333
336 static void
337 MvScale (TMVB& mv,
338 const std::vector<ScalarType>& alpha)
339 {
340 for (unsigned int i=0; i<alpha.size(); i++) {
341 Thyra::scale (alpha[i], mv.col(i).ptr());
342 }
343 }
344
346
349
353 static void MvNorm( const TMVB & mv, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec )
354 { Thyra::norms_2(mv,Teuchos::arrayViewFromVector( normvec )); }
355
357
360
363 static void SetBlock( const TMVB & A, const std::vector<int>& index, TMVB & mv )
364 {
365 // Extract the "numvecs" columns of mv indicated by the index vector.
366 int numvecs = index.size();
367 std::vector<int> indexA(numvecs);
368 int numAcols = A.domain()->dim();
369 for (int i=0; i<numvecs; i++) {
370 indexA[i] = i;
371 }
372 // Thyra::assign requires that both arguments have the same number of
373 // vectors. Enforce this, by shrinking one to match the other.
374 if ( numAcols < numvecs ) {
375 // A does not have enough columns to satisfy index_plus. Shrink
376 // index_plus.
377 numvecs = numAcols;
378 }
379 else if ( numAcols > numvecs ) {
380 numAcols = numvecs;
381 indexA.resize( numAcols );
382 }
383 // create a view to the relevant part of the source multivector
384 Teuchos::RCP< const TMVB > relsource = A.subView( Teuchos::arrayViewFromVector( indexA ) );
385 // create a view to the relevant part of the destination multivector
386 Teuchos::RCP< TMVB > reldest = mv.subView( Teuchos::arrayViewFromVector( index ) );
387 // copy the data to the destination multivector subview
388 Thyra::assign (Teuchos::outArg (*reldest), *relsource);
389 }
390
391 static void
392 SetBlock (const TMVB& A, const Teuchos::Range1D& index, TMVB& mv)
393 {
394 const int numColsA = A.domain()->dim();
395 const int numColsMv = mv.domain()->dim();
396 // 'index' indexes into mv; it's the index set of the target.
397 const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
398 // We can't take more columns out of A than A has.
399 const bool validSource = index.size() <= numColsA;
400
401 if (! validIndex || ! validSource)
402 {
403 std::ostringstream os;
404 os << "Anasazi::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> "
405 ">::SetBlock(A, [" << index.lbound() << ", " << index.ubound()
406 << "], mv): ";
407 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
408 os.str() << "Range lower bound must be nonnegative.");
409 TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
410 os.str() << "Range upper bound must be less than "
411 "the number of columns " << numColsA << " in the "
412 "'mv' output argument.");
413 TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
414 os.str() << "Range must have no more elements than"
415 " the number of columns " << numColsA << " in the "
416 "'A' input argument.");
417 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
418 }
419
420 // View of the relevant column(s) of the target multivector mv.
421 // We avoid view creation overhead by only creating a view if
422 // the index range is different than [0, (# columns in mv) - 1].
423 Teuchos::RCP<TMVB> mv_view;
424 if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
425 mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
426 else
427 mv_view = mv.subView (index);
428
429 // View of the relevant column(s) of the source multivector A.
430 // If A has fewer columns than mv_view, then create a view of
431 // the first index.size() columns of A.
432 Teuchos::RCP<const TMVB> A_view;
433 if (index.size() == numColsA)
434 A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
435 else
436 A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
437
438 // Copy the data to the destination multivector.
439 Thyra::assign (Teuchos::outArg (*mv_view), *A_view);
440 }
441
442 static void
443 Assign (const TMVB& A, TMVB& mv)
444 {
445 using Teuchos::ptr;
446 using Teuchos::Range1D;
447 using Teuchos::RCP;
448
449 const int numColsA = A.domain()->dim();
450 const int numColsMv = mv.domain()->dim();
451 if (numColsA > numColsMv) {
452 const std::string prefix ("Anasazi::MultiVecTraits<Scalar, "
453 "Thyra::MultiVectorBase<Scalar>"
454 " >::Assign(A, mv): ");
455 TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
456 prefix << "Input multivector 'A' has "
457 << numColsA << " columns, but output multivector "
458 "'mv' has only " << numColsMv << " columns.");
459 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
460 }
461 // Copy the data to the destination multivector.
462 if (numColsA == numColsMv) {
463 Thyra::assign (outArg (mv), A);
464 }
465 else {
466 RCP<TMVB> mv_view = CloneViewNonConst (mv, Range1D(0, numColsA-1));
467 Thyra::assign (outArg (*mv_view), A);
468 }
469 }
470
473 static void MvRandom( TMVB & mv )
474 {
475 // Thyra::randomize generates via a uniform distribution on [l,u]
476 // We will use this to generate on [-1,1]
477 Thyra::randomize(-Teuchos::ScalarTraits<ScalarType>::one(),
478 Teuchos::ScalarTraits<ScalarType>::one(),
479 Teuchos::outArg (mv));
480 }
481
484 static void
485 MvInit (TMVB& mv,
486 ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
487 {
488 Thyra::assign (Teuchos::outArg (mv), alpha);
489 }
490
492
495
498 static void MvPrint( const TMVB & mv, std::ostream& os )
499 {
500 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&os,false));
501 out->setf(std::ios_base::scientific);
502 out->precision(16);
503 mv.describe(*out,Teuchos::VERB_EXTREME);
504 }
505
507
508 };
509
510
512 //
513 // Implementation of the Anasazi::OperatorTraits for Thyra::LinearOpBase
514 //
516
526 template <class ScalarType>
527 class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
528 {
529 public:
530
534 static void Apply ( const Thyra::LinearOpBase< ScalarType >& Op, const Thyra::MultiVectorBase< ScalarType > & x, Thyra::MultiVectorBase< ScalarType > & y )
535 {
536 Op.apply(Thyra::NOTRANS,x,Teuchos::outArg (y), Teuchos::ScalarTraits<ScalarType>::one(),Teuchos::ScalarTraits<ScalarType>::zero());
537 }
538
539 };
540
541} // end of Anasazi namespace
542
543#endif
544// end of file ANASAZI_THYRA_ADAPTER_HPP
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
static void MvTimesMatAddMv(const ScalarType alpha, const TMVB &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, TMVB &mv)
Update mv with .
static void MvRandom(TMVB &mv)
Replace the vectors in mv with random vectors.
static Teuchos::RCP< const TMVB > CloneView(const TMVB &mv, const std::vector< int > &index)
Creates a new const MultiVectorBase that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< TMVB > CloneViewNonConst(TMVB &mv, const std::vector< int > &index)
Creates a new MultiVectorBase that shares the selected contents of mv (shallow copy).
static void SetBlock(const TMVB &A, const std::vector< int > &index, TMVB &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index.
static void MvTransMv(const ScalarType alpha, const TMVB &A, const TMVB &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv, const std::vector< int > &index)
Creates a new MultiVectorBase and copies the selected contents of mv into the new vector (deep copy).
static Teuchos::RCP< TMVB > Clone(const TMVB &mv, const int numvecs)
Creates a new empty MultiVectorBase containing numvecs columns.
static void MvInit(TMVB &mv, ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static ptrdiff_t GetGlobalLength(const TMVB &mv)
Obtain the vector length of mv.
static void MvAddMv(const ScalarType alpha, const TMVB &A, const ScalarType beta, const TMVB &B, TMVB &mv)
Replace mv with .
static void MvScale(TMVB &mv, const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv)
Creates a new MultiVectorBase and copies contents of mv into the new vector (deep copy).
static void MvDot(const TMVB &mv, const TMVB &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void MvScale(TMVB &mv, const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
static int GetNumberVecs(const TMVB &mv)
Obtain the number of vectors in mv.
static void MvPrint(const TMVB &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
static void MvNorm(const TMVB &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ,...
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void Assign(const MV &A, MV &mv)
mv := A
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index.
static void Apply(const Thyra::LinearOpBase< ScalarType > &Op, const Thyra::MultiVectorBase< ScalarType > &x, Thyra::MultiVectorBase< ScalarType > &y)
This method takes the MultiVectorBase x and applies the LinearOpBase Op to it resulting in the MultiV...
Virtual base class which defines basic traits for the operator type.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.