Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziTpetraAdapter.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
42#ifndef ANASAZI_TPETRA_ADAPTER_HPP
43#define ANASAZI_TPETRA_ADAPTER_HPP
44
81
82#include <Tpetra_MultiVector.hpp>
83#include <Tpetra_Operator.hpp>
84
85#include <Teuchos_Array.hpp>
86#include <Teuchos_Assert.hpp>
87#include <Teuchos_DefaultSerialComm.hpp>
88#include <Teuchos_CommHelpers.hpp>
89#include <Teuchos_ScalarTraits.hpp>
90#include <Teuchos_FancyOStream.hpp>
91
92#include <AnasaziConfigDefs.hpp>
93#include <AnasaziTypes.hpp>
97
98#ifdef HAVE_ANASAZI_TSQR
99# include <Tpetra_TsqrAdaptor.hpp>
100#endif // HAVE_ANASAZI_TSQR
101
102
103namespace Anasazi {
104
116 template<class Scalar, class LO, class GO, class Node>
117 class MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node> > {
118 typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
119 public:
125 static Teuchos::RCP<MV> Clone (const MV& X, const int numVecs) {
126 Teuchos::RCP<MV> Y (new MV (X.getMap (), numVecs, false));
127 Y->setCopyOrView (Teuchos::View);
128 return Y;
129 }
130
132 static Teuchos::RCP<MV> CloneCopy (const MV& X)
133 {
134 // Make a deep copy of X. The one-argument copy constructor
135 // does a shallow copy by default; the second argument tells it
136 // to do a deep copy.
137 Teuchos::RCP<MV> X_copy (new MV (X, Teuchos::Copy));
138 // Make Tpetra::MultiVector use the new view semantics. This is
139 // a no-op for the Kokkos refactor version of Tpetra; it only
140 // does something for the "classic" version of Tpetra. This
141 // shouldn't matter because Belos only handles MV through RCP
142 // and through this interface anyway, but it doesn't hurt to set
143 // it and make sure that it works.
144 X_copy->setCopyOrView (Teuchos::View);
145 return X_copy;
146 }
147
159 static Teuchos::RCP<MV>
160 CloneCopy (const MV& mv, const std::vector<int>& index)
161 {
162#ifdef HAVE_TPETRA_DEBUG
163 const char fnName[] = "Anasazi::MultiVecTraits::CloneCopy(mv,index)";
164 const size_t inNumVecs = mv.getNumVectors ();
165 TEUCHOS_TEST_FOR_EXCEPTION(
166 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
167 std::runtime_error, fnName << ": All indices must be nonnegative.");
168 TEUCHOS_TEST_FOR_EXCEPTION(
169 index.size () > 0 &&
170 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
171 std::runtime_error,
172 fnName << ": All indices must be strictly less than the number of "
173 "columns " << inNumVecs << " of the input multivector mv.");
174#endif // HAVE_TPETRA_DEBUG
175
176 // Tpetra wants an array of size_t, not of int.
177 Teuchos::Array<size_t> columns (index.size ());
178 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
179 columns[j] = index[j];
180 }
181 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
182 // continuous column index range in MultiVector::subCopy, so we
183 // don't have to check here.
184 Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
185 X_copy->setCopyOrView (Teuchos::View);
186 return X_copy;
187 }
188
195 static Teuchos::RCP<MV>
196 CloneCopy (const MV& mv, const Teuchos::Range1D& index)
197 {
198 const bool validRange = index.size() > 0 &&
199 index.lbound() >= 0 &&
200 index.ubound() < GetNumberVecs(mv);
201 if (! validRange) { // invalid range; generate error message
202 std::ostringstream os;
203 os << "Anasazi::MultiVecTraits::CloneCopy(mv,index=["
204 << index.lbound() << "," << index.ubound() << "]): ";
205 TEUCHOS_TEST_FOR_EXCEPTION(
206 index.size() == 0, std::invalid_argument,
207 os.str() << "Empty index range is not allowed.");
208 TEUCHOS_TEST_FOR_EXCEPTION(
209 index.lbound() < 0, std::invalid_argument,
210 os.str() << "Index range includes negative index/ices, which is not "
211 "allowed.");
212 TEUCHOS_TEST_FOR_EXCEPTION(
213 index.ubound() >= GetNumberVecs(mv), std::invalid_argument,
214 os.str() << "Index range exceeds number of vectors "
215 << mv.getNumVectors() << " in the input multivector.");
216 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
217 os.str() << "Should never get here!");
218 }
219 Teuchos::RCP<MV> X_copy = mv.subCopy (index);
220 X_copy->setCopyOrView (Teuchos::View);
221 return X_copy;
222 }
223
224 static Teuchos::RCP<MV>
225 CloneViewNonConst (MV& mv, const std::vector<int>& index)
226 {
227#ifdef HAVE_TPETRA_DEBUG
228 const char fnName[] = "Anasazi::MultiVecTraits::CloneViewNonConst(mv,index)";
229 const size_t numVecs = mv.getNumVectors ();
230 TEUCHOS_TEST_FOR_EXCEPTION(
231 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
232 std::invalid_argument,
233 fnName << ": All indices must be nonnegative.");
234 TEUCHOS_TEST_FOR_EXCEPTION(
235 index.size () > 0 &&
236 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
237 std::invalid_argument,
238 fnName << ": All indices must be strictly less than the number of "
239 "columns " << numVecs << " in the input MultiVector mv.");
240#endif // HAVE_TPETRA_DEBUG
241
242 // Tpetra wants an array of size_t, not of int.
243 Teuchos::Array<size_t> columns (index.size ());
244 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
245 columns[j] = index[j];
246 }
247 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
248 // continuous column index range in
249 // MultiVector::subViewNonConst, so we don't have to check here.
250 Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
251 X_view->setCopyOrView (Teuchos::View);
252 return X_view;
253 }
254
255 static Teuchos::RCP<MV>
256 CloneViewNonConst (MV& mv, const Teuchos::Range1D& index)
257 {
258 // NOTE (mfh 11 Jan 2011) We really should check for possible
259 // overflow of int here. However, the number of columns in a
260 // multivector typically fits in an int.
261 const int numCols = static_cast<int> (mv.getNumVectors());
262 const bool validRange = index.size() > 0 &&
263 index.lbound() >= 0 && index.ubound() < numCols;
264 if (! validRange) {
265 std::ostringstream os;
266 os << "Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
267 << index.lbound() << ", " << index.ubound() << "]): ";
268 TEUCHOS_TEST_FOR_EXCEPTION(
269 index.size() == 0, std::invalid_argument,
270 os.str() << "Empty index range is not allowed.");
271 TEUCHOS_TEST_FOR_EXCEPTION(
272 index.lbound() < 0, std::invalid_argument,
273 os.str() << "Index range includes negative inde{x,ices}, which is "
274 "not allowed.");
275 TEUCHOS_TEST_FOR_EXCEPTION(
276 index.ubound() >= numCols, std::invalid_argument,
277 os.str() << "Index range exceeds number of vectors " << numCols
278 << " in the input multivector.");
279 TEUCHOS_TEST_FOR_EXCEPTION(
280 true, std::logic_error,
281 os.str() << "Should never get here!");
282 }
283 Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
284 X_view->setCopyOrView (Teuchos::View);
285 return X_view;
286 }
287
288 static Teuchos::RCP<const MV>
289 CloneView (const MV& mv, const std::vector<int>& index)
290 {
291#ifdef HAVE_TPETRA_DEBUG
292 const char fnName[] = "Belos::MultiVecTraits<Scalar, "
293 "Tpetra::MultiVector<...> >::CloneView(mv,index)";
294 const size_t numVecs = mv.getNumVectors ();
295 TEUCHOS_TEST_FOR_EXCEPTION(
296 *std::min_element (index.begin (), index.end ()) < 0,
297 std::invalid_argument,
298 fnName << ": All indices must be nonnegative.");
299 TEUCHOS_TEST_FOR_EXCEPTION(
300 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
301 std::invalid_argument,
302 fnName << ": All indices must be strictly less than the number of "
303 "columns " << numVecs << " in the input MultiVector mv.");
304#endif // HAVE_TPETRA_DEBUG
305
306 // Tpetra wants an array of size_t, not of int.
307 Teuchos::Array<size_t> columns (index.size ());
308 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
309 columns[j] = index[j];
310 }
311 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
312 // continuous column index range in MultiVector::subView, so we
313 // don't have to check here.
314 Teuchos::RCP<const MV> X_view = mv.subView (columns);
315 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
316 return X_view;
317 }
318
319 static Teuchos::RCP<const MV>
320 CloneView (const MV& mv, const Teuchos::Range1D& index)
321 {
322 // NOTE (mfh 11 Jan 2011) We really should check for possible
323 // overflow of int here. However, the number of columns in a
324 // multivector typically fits in an int.
325 const int numCols = static_cast<int> (mv.getNumVectors());
326 const bool validRange = index.size() > 0 &&
327 index.lbound() >= 0 && index.ubound() < numCols;
328 if (! validRange) {
329 std::ostringstream os;
330 os << "Anasazi::MultiVecTraits::CloneView(mv, index=["
331 << index.lbound () << ", " << index.ubound() << "]): ";
332 TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
333 os.str() << "Empty index range is not allowed.");
334 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
335 os.str() << "Index range includes negative index/ices, which is not "
336 "allowed.");
337 TEUCHOS_TEST_FOR_EXCEPTION(
338 index.ubound() >= numCols, std::invalid_argument,
339 os.str() << "Index range exceeds number of vectors " << numCols
340 << " in the input multivector.");
341 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
342 os.str() << "Should never get here!");
343 }
344 Teuchos::RCP<const MV> X_view = mv.subView (index);
345 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
346 return X_view;
347 }
348
349 static ptrdiff_t GetGlobalLength( const MV& mv ) {
350 return static_cast<ptrdiff_t> (mv.getGlobalLength ());
351 }
352
353 static int GetNumberVecs (const MV& mv) {
354 return static_cast<int> (mv.getNumVectors ());
355 }
356
357 static void
358 MvTimesMatAddMv (Scalar alpha,
359 const MV& A,
360 const Teuchos::SerialDenseMatrix<int, Scalar>& B,
361 Scalar beta,
362 MV& mv)
363 {
364 using Teuchos::ArrayView;
365 using Teuchos::Comm;
366 using Teuchos::rcpFromRef;
367 typedef Tpetra::Map<LO, GO, Node> map_type;
368
369#ifdef HAVE_ANASAZI_TPETRA_TIMERS
370 const std::string timerName ("Anasazi::MVT::MvTimesMatAddMv");
371 Teuchos::RCP<Teuchos::Time> timer =
372 Teuchos::TimeMonitor::lookupCounter (timerName);
373 if (timer.is_null ()) {
374 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
375 }
376 TEUCHOS_TEST_FOR_EXCEPTION(
377 timer.is_null (), std::logic_error,
378 "Anasazi::MultiVecTraits::MvTimesMatAddMv: "
379 "Failed to look up timer \"" << timerName << "\". "
380 "Please report this bug to the Belos developers.");
381
382 // This starts the timer. It will be stopped on scope exit.
383 Teuchos::TimeMonitor timeMon (*timer);
384#endif
385
386 // Check if B is 1-by-1, in which case we can just call update()
387 if (B.numRows () == 1 && B.numCols () == 1) {
388 mv.update (alpha*B(0,0), A, beta);
389 return;
390 }
391
392 // Create local map
393 Teuchos::SerialComm<int> serialComm;
394 map_type LocalMap (B.numRows (), A.getMap ()->getIndexBase (),
395 rcpFromRef<const Comm<int> > (serialComm),
396 Tpetra::LocallyReplicated);
397 // encapsulate Teuchos::SerialDenseMatrix data in ArrayView
398 ArrayView<const Scalar> Bvalues (B.values (), B.stride () * B.numCols ());
399 // create locally replicated MultiVector with a copy of this data
400 MV B_mv (rcpFromRef (LocalMap), Bvalues, B.stride (), B.numCols ());
401 mv.multiply (Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
402 }
403
411 static void
412 MvAddMv (Scalar alpha,
413 const MV& A,
414 Scalar beta,
415 const MV& B,
416 MV& mv)
417 {
418 mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
419 }
420
421 static void MvScale (MV& mv, Scalar alpha) {
422 mv.scale (alpha);
423 }
424
425 static void MvScale (MV& mv, const std::vector<Scalar>& alphas) {
426 mv.scale (alphas);
427 }
428
429 static void
430 MvTransMv (const Scalar alpha,
431 const MV& A,
432 const MV& B,
433 Teuchos::SerialDenseMatrix<int,Scalar>& C)
434 {
435 using Tpetra::LocallyReplicated;
436 using Teuchos::Comm;
437 using Teuchos::RCP;
438 using Teuchos::rcp;
439 using Teuchos::TimeMonitor;
440 typedef Tpetra::Map<LO,GO,Node> map_type;
441
442#ifdef HAVE_ANASAZI_TPETRA_TIMERS
443 const std::string timerName ("Anasazi::MVT::MvTransMv");
444 RCP<Teuchos::Time> timer = TimeMonitor::lookupCounter (timerName);
445 if (timer.is_null ()) {
446 timer = TimeMonitor::getNewCounter (timerName);
447 }
448 TEUCHOS_TEST_FOR_EXCEPTION(
449 timer.is_null (), std::logic_error, "Anasazi::MvTransMv: "
450 "Failed to look up timer \"" << timerName << "\". "
451 "Please report this bug to the Belos developers.");
452
453 // This starts the timer. It will be stopped on scope exit.
454 TimeMonitor timeMon (*timer);
455#endif // HAVE_ANASAZI_TPETRA_TIMERS
456
457 // Form alpha * A^H * B, then copy into the SerialDenseMatrix.
458 // We will create a multivector C_mv from a a local map. This
459 // map has a serial comm, the purpose being to short-circuit the
460 // MultiVector::reduce() call at the end of
461 // MultiVector::multiply(). Otherwise, the reduced multivector
462 // data would be copied back to the GPU, only to turn around and
463 // have to get it back here. This saves us a round trip for
464 // this data.
465 const int numRowsC = C.numRows ();
466 const int numColsC = C.numCols ();
467 const int strideC = C.stride ();
468
469 // Check if numRowsC == numColsC == 1, in which case we can call dot()
470 if (numRowsC == 1 && numColsC == 1) {
471 if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
472 // Short-circuit, as required by BLAS semantics.
473 C(0,0) = alpha;
474 return;
475 }
476 A.dot (B, Teuchos::ArrayView<Scalar> (C.values (), 1));
477 if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
478 C(0,0) *= alpha;
479 }
480 return;
481 }
482
483 // get comm
484 RCP<const Comm<int> > pcomm = A.getMap ()->getComm ();
485
486 // create local map with comm
487 RCP<const map_type> LocalMap =
488 rcp (new map_type (numRowsC, 0, pcomm, LocallyReplicated));
489 // create local multivector to hold the result
490 const bool INIT_TO_ZERO = true;
491 MV C_mv (LocalMap, numColsC, INIT_TO_ZERO);
492
493 // multiply result into local multivector
494 C_mv.multiply (Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, alpha, A, B,
495 Teuchos::ScalarTraits<Scalar>::zero ());
496
497 // create arrayview encapsulating the Teuchos::SerialDenseMatrix
498 Teuchos::ArrayView<Scalar> C_view (C.values (), strideC * numColsC);
499
500 // No accumulation to do; simply extract the multivector data
501 // into C. Extract a copy of the result into the array view
502 // (and therefore, the SerialDenseMatrix).
503 C_mv.get1dCopy (C_view, strideC);
504 }
505
507 static void
508 MvDot (const MV& A, const MV& B, std::vector<Scalar> &dots)
509 {
510 const size_t numVecs = A.getNumVectors ();
511 TEUCHOS_TEST_FOR_EXCEPTION(
512 numVecs != B.getNumVectors (), std::invalid_argument,
513 "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
514 "A and B must have the same number of columns. "
515 "A has " << numVecs << " column(s), "
516 "but B has " << B.getNumVectors () << " column(s).");
517#ifdef HAVE_TPETRA_DEBUG
518 TEUCHOS_TEST_FOR_EXCEPTION(
519 dots.size() < numVecs, std::invalid_argument,
520 "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
521 "The output array 'dots' must have room for all dot products. "
522 "A and B each have " << numVecs << " column(s), "
523 "but 'dots' only has " << dots.size() << " entry(/ies).");
524#endif // HAVE_TPETRA_DEBUG
525
526 Teuchos::ArrayView<Scalar> av (dots);
527 A.dot (B, av (0, numVecs));
528 }
529
531 static void
532 MvNorm (const MV& mv,
533 std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
534 {
535 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
536#ifdef HAVE_TPETRA_DEBUG
537 TEUCHOS_TEST_FOR_EXCEPTION(
538 normvec.size () < static_cast<std::vector<int>::size_type> (mv.getNumVectors ()),
539 std::invalid_argument,
540 "Anasazi::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
541 "argument must have at least as many entries as the number of vectors "
542 "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
543 << " < mv.getNumVectors() = " << mv.getNumVectors () << ".");
544#endif // HAVE_TPETRA_DEBUG
545 Teuchos::ArrayView<magnitude_type> av (normvec);
546 mv.norm2 (av (0, mv.getNumVectors ()));
547 }
548
549 static void
550 SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
551 {
552 using Teuchos::Range1D;
553 using Teuchos::RCP;
554 const size_t inNumVecs = A.getNumVectors ();
555#ifdef HAVE_TPETRA_DEBUG
556 TEUCHOS_TEST_FOR_EXCEPTION(
557 inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
558 "Anasazi::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
559 "have no more entries as the number of columns in the input MultiVector"
560 " A. A.getNumVectors() = " << inNumVecs << " < index.size () = "
561 << index.size () << ".");
562#endif // HAVE_TPETRA_DEBUG
563 RCP<MV> mvsub = CloneViewNonConst (mv, index);
564 if (inNumVecs > static_cast<size_t> (index.size ())) {
565 RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
566 Tpetra::deep_copy (*mvsub, *Asub);
567 } else {
568 Tpetra::deep_copy (*mvsub, A);
569 }
570 }
571
572 static void
573 SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
574 {
575 // Range1D bounds are signed; size_t is unsigned.
576 // Assignment of Tpetra::MultiVector is a deep copy.
577
578 // Tpetra::MultiVector::getNumVectors() returns size_t. It's
579 // fair to assume that the number of vectors won't overflow int,
580 // since the typical use case of multivectors involves few
581 // columns, but it's friendly to check just in case.
582 const size_t maxInt =
583 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
584 const bool overflow =
585 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
586 if (overflow) {
587 std::ostringstream os;
588 os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
589 << ", " << index.ubound () << "], mv): ";
590 TEUCHOS_TEST_FOR_EXCEPTION(
591 maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
592 "of columns (size_t) in the input MultiVector 'A' overflows int.");
593 TEUCHOS_TEST_FOR_EXCEPTION(
594 maxInt < mv.getNumVectors (), std::range_error, os.str () << "Number "
595 "of columns (size_t) in the output MultiVector 'mv' overflows int.");
596 }
597 // We've already validated the static casts above.
598 const int numColsA = static_cast<int> (A.getNumVectors ());
599 const int numColsMv = static_cast<int> (mv.getNumVectors ());
600 // 'index' indexes into mv; it's the index set of the target.
601 const bool validIndex =
602 index.lbound () >= 0 && index.ubound () < numColsMv;
603 // We can't take more columns out of A than A has.
604 const bool validSource = index.size () <= numColsA;
605
606 if (! validIndex || ! validSource) {
607 std::ostringstream os;
608 os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
609 << ", " << index.ubound () << "], mv): ";
610 TEUCHOS_TEST_FOR_EXCEPTION(
611 index.lbound() < 0, std::invalid_argument,
612 os.str() << "Range lower bound must be nonnegative.");
613 TEUCHOS_TEST_FOR_EXCEPTION(
614 index.ubound() >= numColsMv, std::invalid_argument,
615 os.str() << "Range upper bound must be less than the number of "
616 "columns " << numColsA << " in the 'mv' output argument.");
617 TEUCHOS_TEST_FOR_EXCEPTION(
618 index.size() > numColsA, std::invalid_argument,
619 os.str() << "Range must have no more elements than the number of "
620 "columns " << numColsA << " in the 'A' input argument.");
621 TEUCHOS_TEST_FOR_EXCEPTION(
622 true, std::logic_error, "Should never get here!");
623 }
624
625 // View of the relevant column(s) of the target multivector mv.
626 // We avoid view creation overhead by only creating a view if
627 // the index range is different than [0, (# columns in mv) - 1].
628 Teuchos::RCP<MV> mv_view;
629 if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
630 mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
631 } else {
632 mv_view = CloneViewNonConst (mv, index);
633 }
634
635 // View of the relevant column(s) of the source multivector A.
636 // If A has fewer columns than mv_view, then create a view of
637 // the first index.size() columns of A.
638 Teuchos::RCP<const MV> A_view;
639 if (index.size () == numColsA) {
640 A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
641 } else {
642 A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
643 }
644
645 Tpetra::deep_copy (*mv_view, *A_view);
646 }
647
648 static void Assign (const MV& A, MV& mv)
649 {
650 const char errPrefix[] = "Anasazi::MultiVecTraits::Assign(A, mv): ";
651
652 // Range1D bounds are signed; size_t is unsigned.
653 // Assignment of Tpetra::MultiVector is a deep copy.
654
655 // Tpetra::MultiVector::getNumVectors() returns size_t. It's
656 // fair to assume that the number of vectors won't overflow int,
657 // since the typical use case of multivectors involves few
658 // columns, but it's friendly to check just in case.
659 const size_t maxInt =
660 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
661 const bool overflow =
662 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
663 if (overflow) {
664 TEUCHOS_TEST_FOR_EXCEPTION(
665 maxInt < A.getNumVectors(), std::range_error,
666 errPrefix << "Number of columns in the input multivector 'A' "
667 "(a size_t) overflows int.");
668 TEUCHOS_TEST_FOR_EXCEPTION(
669 maxInt < mv.getNumVectors(), std::range_error,
670 errPrefix << "Number of columns in the output multivector 'mv' "
671 "(a size_t) overflows int.");
672 TEUCHOS_TEST_FOR_EXCEPTION(
673 true, std::logic_error, "Should never get here!");
674 }
675 // We've already validated the static casts above.
676 const int numColsA = static_cast<int> (A.getNumVectors ());
677 const int numColsMv = static_cast<int> (mv.getNumVectors ());
678 if (numColsA > numColsMv) {
679 TEUCHOS_TEST_FOR_EXCEPTION(
680 numColsA > numColsMv, std::invalid_argument,
681 errPrefix << "Input multivector 'A' has " << numColsA << " columns, "
682 "but output multivector 'mv' has only " << numColsMv << " columns.");
683 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
684 }
685 if (numColsA == numColsMv) {
686 Tpetra::deep_copy (mv, A);
687 } else {
688 Teuchos::RCP<MV> mv_view =
689 CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
690 Tpetra::deep_copy (*mv_view, A);
691 }
692 }
693
694 static void MvRandom (MV& mv) {
695 mv.randomize ();
696 }
697
698 static void
699 MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
700 {
701 mv.putScalar (alpha);
702 }
703
704 static void MvPrint (const MV& mv, std::ostream& os) {
705 mv.print (os);
706 }
707
708#ifdef HAVE_ANASAZI_TSQR
711 typedef Tpetra::TsqrAdaptor<Tpetra::MultiVector<Scalar, LO, GO, Node> > tsqr_adaptor_type;
712#endif // HAVE_ANASAZI_TSQR
713 };
714
715
717 template <class Scalar, class LO, class GO, class Node>
718 class OperatorTraits<Scalar,
719 Tpetra::MultiVector<Scalar,LO,GO,Node>,
720 Tpetra::Operator<Scalar,LO,GO,Node> >
721 {
722 public:
723 static void
724 Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
725 const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
726 Tpetra::MultiVector<Scalar,LO,GO,Node>& Y)
727 {
728 Op.apply (X, Y, Teuchos::NO_TRANS);
729 }
730 };
731
732
733template<class ST, class LO, class GO, class NT>
734struct OutputStreamTraits<Tpetra::Operator<ST, LO, GO, NT> > {
735 typedef Tpetra::Operator<ST, LO, GO, NT> operator_type;
736
737 static Teuchos::RCP<Teuchos::FancyOStream>
738 getOutputStream (const operator_type& op, int rootRank = 0)
739 {
740 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
741 Teuchos::RCP<const Teuchos::Comm<int> > comm = (op.getDomainMap())->getComm ();
742
743 // Select minimum MPI rank as the root rank for printing.
744 const int myRank = comm.is_null () ? 0 : comm->getRank ();
745 const int numProcs = comm.is_null () ? 1 : comm->getSize ();
746 if (rootRank < 0)
747 {
748 Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,1,&myRank,&rootRank);
749 }
750
751 // This is irreversible, but that's only a problem if the input std::ostream
752 // is actually a Teuchos::FancyOStream on which this method has been
753 // called before, with a different root rank.
754 fos->setProcRankAndSize (myRank, numProcs);
755 fos->setOutputToRootOnly (rootRank);
756 return fos;
757 }
758};
759
760
761} // end of Anasazi namespace
762
763#endif
764// end of file ANASAZI_TPETRA_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.
Abstract class definition for Anasazi output stream.
Types and exceptions used within Anasazi solvers and interfaces.
static void MvDot(const MV &A, const MV &B, std::vector< Scalar > &dots)
For all columns j of A, set dots[j] := A[j]^T * B[j].
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const Teuchos::Range1D &index)
Create and return a deep copy of the given columns of mv.
static Teuchos::RCP< MV > Clone(const MV &X, const int numVecs)
Create a new MultiVector with numVecs columns.
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const std::vector< int > &index)
Create and return a deep copy of the given columns of mv.
static void MvAddMv(Scalar alpha, const MV &A, Scalar beta, const MV &B, MV &mv)
mv := alpha*A + beta*B
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &normvec)
For all columns j of mv, set normvec[j] = norm(mv[j]).
static Teuchos::RCP< MV > CloneCopy(const MV &X)
Create and return a deep copy of X.
Traits class which defines basic operations on multivectors.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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 ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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.
Virtual base class which defines basic traits for the operator type.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Anasazi's templated virtual class for constructing an operator that can interface with the OperatorTr...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Output managers remove the need for the eigensolver to know any information about the required output...