Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Details_TriDiSolver_def.hpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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
43#ifndef IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP
44#define IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP
45
46#include "Ifpack2_LocalFilter.hpp"
47#include "Teuchos_LAPACK.hpp"
48#include "Tpetra_MultiVector.hpp"
49#include "Tpetra_Map.hpp"
50#include "Tpetra_Import.hpp"
51#include "Tpetra_Export.hpp"
52
53#ifdef HAVE_MPI
54# include <mpi.h>
55# include "Teuchos_DefaultMpiComm.hpp"
56#else
57# include "Teuchos_DefaultSerialComm.hpp"
58#endif // HAVE_MPI
59
60
61namespace Ifpack2 {
62namespace Details {
63
65// Non-stub (full) implementation
67
68template<class MatrixType>
70TriDiSolver (const Teuchos::RCP<const row_matrix_type>& A) :
71 A_ (A),
72 initializeTime_ (0.0),
73 computeTime_ (0.0),
74 applyTime_ (0.0),
75 numInitialize_ (0),
76 numCompute_ (0),
77 numApply_ (0),
78 isInitialized_ (false),
79 isComputed_ (false)
80{}
81
82
83template<class MatrixType>
84Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::map_type>
86 TEUCHOS_TEST_FOR_EXCEPTION(
87 A_.is_null (), std::runtime_error, "Ifpack2::Details::TriDiSolver::"
88 "getDomainMap: The input matrix A is null. Please call setMatrix() with a "
89 "nonnull input matrix before calling this method.");
90 // For an input matrix A, TriDiSolver solves Ax=b for x.
91 // Thus, its Maps are reversed from those of the input matrix.
92 return A_->getRangeMap ();
93}
94
95
96template<class MatrixType>
97Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::map_type>
99 TEUCHOS_TEST_FOR_EXCEPTION(
100 A_.is_null (), std::runtime_error, "Ifpack2::Details::TriDiSolver::"
101 "getRangeMap: The input matrix A is null. Please call setMatrix() with a "
102 "nonnull input matrix before calling this method.");
103 // For an input matrix A, TriDiSolver solves Ax=b for x.
104 // Thus, its Maps are reversed from those of the input matrix.
105 return A_->getDomainMap ();
106}
107
108
109template<class MatrixType>
110void
112setParameters (const Teuchos::ParameterList& params) {
113 (void) params; // this preconditioner doesn't currently take any parameters
114}
115
116
117template<class MatrixType>
118bool
120 return isInitialized_;
121}
122
123
124template<class MatrixType>
125bool
127 return isComputed_;
128}
129
130
131template<class MatrixType>
132int
134 return numInitialize_;
135}
136
137
138template<class MatrixType>
139int
141 return numCompute_;
142}
143
144
145template<class MatrixType>
146int
148 return numApply_;
149}
150
151
152template<class MatrixType>
153double
155 return initializeTime_;
156}
157
158
159template<class MatrixType>
160double
162 return computeTime_;
163}
164
165
166template<class MatrixType>
167double
169 return applyTime_;
170}
171
172
173template<class MatrixType>
174Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::row_matrix_type>
176 return A_;
177}
178
179
180template<class MatrixType>
182reset ()
183{
184 isInitialized_ = false;
185 isComputed_ = false;
186 A_local_ = Teuchos::null;
187 A_local_tridi_.reshape (0);
188 ipiv_.resize (0);
189}
190
191
192template<class MatrixType>
194setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
195{
196 // It's legitimate to call setMatrix() with a null input. This has
197 // the effect of resetting the preconditioner's internal state.
198 if (! A_.is_null ()) {
199 const global_size_t numRows = A->getRangeMap ()->getGlobalNumElements ();
200 const global_size_t numCols = A->getDomainMap ()->getGlobalNumElements ();
201 TEUCHOS_TEST_FOR_EXCEPTION(
202 numRows != numCols, std::invalid_argument, "Ifpack2::Details::TriDiSolver::"
203 "setMatrix: Input matrix must be (globally) square. "
204 "The matrix you provided is " << numRows << " by " << numCols << ".");
205 }
206 // Clear any previously computed objects.
207 reset ();
208
209 // Now that we've cleared the state, we can keep the matrix.
210 A_ = A;
211}
212
213
214template<class MatrixType>
216{
217 using Teuchos::Comm;
218 using Teuchos::null;
219 using Teuchos::RCP;
220 using Teuchos::rcp;
221 using Teuchos::Time;
222 using Teuchos::TimeMonitor;
223 const std::string timerName ("Ifpack2::Details::TriDiSolver::initialize");
224
225 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
226 if (timer.is_null ()) {
227 timer = TimeMonitor::getNewCounter (timerName);
228 }
229
230 double startTime = timer->wallTime();
231
232 { // Begin timing here.
233 Teuchos::TimeMonitor timeMon (*timer);
234
235 TEUCHOS_TEST_FOR_EXCEPTION(
236 A_.is_null (), std::runtime_error, "Ifpack2::Details::TriDiSolver::"
237 "initialize: The input matrix A is null. Please call setMatrix() "
238 "with a nonnull input before calling this method.");
239
240 TEUCHOS_TEST_FOR_EXCEPTION(
241 ! A_->hasColMap (), std::invalid_argument, "Ifpack2::Details::TriDiSolver: "
242 "The constructor's input matrix must have a column Map, "
243 "so that it has local indices.");
244
245 // Clear any previously computed objects.
246 reset ();
247
248 // Make the local filter of the input matrix A.
249 if (A_->getComm ()->getSize () > 1) {
250 A_local_ = rcp (new LocalFilter<row_matrix_type> (A_));
251 } else {
252 A_local_ = A_;
253 }
254
255 TEUCHOS_TEST_FOR_EXCEPTION(
256 A_local_.is_null (), std::logic_error, "Ifpack2::Details::TriDiSolver::"
257 "initialize: A_local_ is null after it was supposed to have been "
258 "initialized. Please report this bug to the Ifpack2 developers.");
259
260 // Allocate the TriDi local matrix and the pivot array.
261 const size_t numRows = A_local_->getLocalNumRows ();
262 const size_t numCols = A_local_->getLocalNumCols ();
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 numRows != numCols, std::logic_error, "Ifpack2::Details::TriDiSolver::"
265 "initialize: Local filter matrix is not square. This should never happen. "
266 "Please report this bug to the Ifpack2 developers.");
267 A_local_tridi_.reshape (numRows);
268 ipiv_.resize (numRows);
269 std::fill (ipiv_.begin (), ipiv_.end (), 0);
270
271 isInitialized_ = true;
272 ++numInitialize_;
273 }
274
275 initializeTime_ += (timer->wallTime() - startTime);
276}
277
278
279template<class MatrixType>
282
283
284template<class MatrixType>
286{
287 using Teuchos::RCP;
288 const std::string timerName ("Ifpack2::Details::TriDiSolver::compute");
289
290 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
291 if (timer.is_null ()) {
292 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
293 }
294
295 double startTime = timer->wallTime();
296
297 // Begin timing here.
298 {
299 Teuchos::TimeMonitor timeMon (*timer);
300 TEUCHOS_TEST_FOR_EXCEPTION(
301 A_.is_null (), std::runtime_error, "Ifpack2::Details::TriDiSolver::"
302 "compute: The input matrix A is null. Please call setMatrix() with a "
303 "nonnull input, then call initialize(), before calling this method.");
304
305 TEUCHOS_TEST_FOR_EXCEPTION(
306 A_local_.is_null (), std::logic_error, "Ifpack2::Details::TriDiSolver::"
307 "compute: A_local_ is null. Please report this bug to the Ifpack2 "
308 "developers.");
309
310 isComputed_ = false;
311 if (! this->isInitialized ()) {
312 this->initialize ();
313 }
314 extract (A_local_tridi_, *A_local_); // extract the tridi local matrix
315
316 factor (A_local_tridi_, ipiv_ ()); // factor the tridi local matrix
317
318 isComputed_ = true;
319 ++numCompute_;
320 }
321 computeTime_ += (timer->wallTime() - startTime);
322}
323
324template<class MatrixType>
325void TriDiSolver<MatrixType, false>::factor (Teuchos::SerialTriDiMatrix<int, scalar_type>& A,
326 const Teuchos::ArrayView<int>& ipiv)
327{
328 // Fill the LU permutation array with zeros.
329 std::fill (ipiv.begin (), ipiv.end (), 0);
330
331 Teuchos::LAPACK<int, scalar_type> lapack;
332 int INFO = 0;
333 lapack.GTTRF (A.numRowsCols (),
334 A.DL(),
335 A.D(),
336 A.DU(),
337 A.DU2(),
338 ipiv.getRawPtr (), &INFO);
339 // INFO < 0 is a bug.
340 TEUCHOS_TEST_FOR_EXCEPTION(
341 INFO < 0, std::logic_error, "Ifpack2::Details::TriDiSolver::factor: "
342 "LAPACK's _GTTRF (tridiagonal LU factorization with partial pivoting) "
343 "was called incorrectly. INFO = " << INFO << " < 0. "
344 "Please report this bug to the Ifpack2 developers.");
345 // INFO > 0 means the matrix is singular. This is probably an issue
346 // either with the choice of rows the rows we extracted, or with the
347 // input matrix itself.
348 TEUCHOS_TEST_FOR_EXCEPTION(
349 INFO > 0, std::runtime_error, "Ifpack2::Details::TriDiSolver::factor: "
350 "LAPACK's _GTTRF (tridiagonal LU factorization with partial pivoting) "
351 "reports that the computed U factor is exactly singular. U(" << INFO <<
352 "," << INFO << ") (one-based index i) is exactly zero. This probably "
353 "means that the input matrix has a singular diagonal block.");
354}
355
356
357template<class MatrixType>
358void TriDiSolver<MatrixType, false>::
359applyImpl (const MV& X,
360 MV& Y,
361 const Teuchos::ETransp mode,
362 const scalar_type alpha,
363 const scalar_type beta) const
364{
365 using Teuchos::ArrayRCP;
366 using Teuchos::RCP;
367 using Teuchos::rcp;
368 using Teuchos::rcpFromRef;
369 using Teuchos::CONJ_TRANS;
370 using Teuchos::TRANS;
371
372 const int numVecs = static_cast<int> (X.getNumVectors ());
373 if (alpha == STS::zero ()) { // don't need to solve the linear system
374 if (beta == STS::zero ()) {
375 // Use BLAS AXPY semantics for beta == 0: overwrite, clobbering
376 // any Inf or NaN values in Y (rather than multiplying them by
377 // zero, resulting in NaN values).
378 Y.putScalar (STS::zero ());
379 }
380 else { // beta != 0
381 Y.scale (STS::zero ());
382 }
383 }
384 else { // alpha != 0; must solve the linear system
385 Teuchos::LAPACK<int, scalar_type> lapack;
386 // If beta is nonzero, Y is not constant stride, or alpha != 1, we
387 // have to use a temporary output multivector Y_tmp. It gets a
388 // copy of alpha*X, since GETRS overwrites its (multi)vector input
389 // with its output.
390 RCP<MV> Y_tmp;
391 if (beta == STS::zero () && Y.isConstantStride () && alpha == STS::one ()) {
392 deep_copy(Y, X);
393 Y_tmp = rcpFromRef (Y);
394 }
395 else {
396 Y_tmp = rcp (new MV (createCopy(X))); // constructor copies X
397 if (alpha != STS::one ()) {
398 Y_tmp->scale (alpha);
399 }
400 }
401 const int Y_stride = static_cast<int> (Y_tmp->getStride ());
402 ArrayRCP<scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
403 scalar_type* const Y_ptr = Y_view.getRawPtr ();
404 int INFO = 0;
405 const char trans =
406 (mode == CONJ_TRANS ? 'C' : (mode == TRANS ? 'T' : 'N'));
407 lapack.GTTRS (trans, A_local_tridi_.numRowsCols(), numVecs,
408 A_local_tridi_.DL(),
409 A_local_tridi_.D(),
410 A_local_tridi_.DU(),
411 A_local_tridi_.DU2(),
412 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
413 TEUCHOS_TEST_FOR_EXCEPTION(
414 INFO != 0, std::runtime_error, "Ifpack2::Details::TriDiSolver::"
415 "applyImpl: LAPACK's _GTTRS (tridiagonal solve using LU factorization "
416 "with partial pivoting) failed with INFO = " << INFO << " != 0.");
417
418 if (beta != STS::zero ()) {
419 Y.update (alpha, *Y_tmp, beta);
420 }
421 else if (! Y.isConstantStride ()) {
422 deep_copy(Y, *Y_tmp);
423 }
424 }
425}
426
427
428template<class MatrixType>
430apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
431 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
432 Teuchos::ETransp mode,
433 scalar_type alpha,
434 scalar_type beta) const
435{
436 using Teuchos::ArrayView;
437 using Teuchos::as;
438 using Teuchos::RCP;
439 using Teuchos::rcp;
440 using Teuchos::rcpFromRef;
441
442 const std::string timerName ("Ifpack2::Details::TriDiSolver::apply");
443 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
444 if (timer.is_null ()) {
445 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
446 }
447
448 double startTime = timer->wallTime();
449
450 // Begin timing here.
451 {
452 Teuchos::TimeMonitor timeMon (*timer);
453
454 TEUCHOS_TEST_FOR_EXCEPTION(
455 ! isComputed_, std::runtime_error, "Ifpack2::Details::TriDiSolver::apply: "
456 "You must have called the compute() method before you may call apply(). "
457 "You may call the apply() method as many times as you want after calling "
458 "compute() once, but you must have called compute() at least once.");
459
460 const size_t numVecs = X.getNumVectors ();
461
462 TEUCHOS_TEST_FOR_EXCEPTION(
463 numVecs != Y.getNumVectors (), std::runtime_error,
464 "Ifpack2::Details::TriDiSolver::apply: X and Y have different numbers "
465 "of vectors. X has " << X.getNumVectors () << ", but Y has "
466 << X.getNumVectors () << ".");
467
468 if (numVecs == 0) {
469 return; // done! nothing to do
470 }
471
472 // Set up "local" views of X and Y.
473 RCP<const MV> X_local;
474 RCP<MV> Y_local;
475 const bool multipleProcs = (A_->getRowMap ()->getComm ()->getSize () >= 1);
476 if (multipleProcs) {
477 // Interpret X and Y as "local" multivectors, that is, in the
478 // local filter's domain resp. range Maps. "Interpret" means that
479 // we create views with different Maps; we don't have to copy.
480 X_local = X.offsetView (A_local_->getDomainMap (), 0);
481 Y_local = Y.offsetViewNonConst (A_local_->getRangeMap (), 0);
482 }
483 else { // only one process in A_'s communicator
484 // X and Y are already "local"; no need to set up local views.
485 X_local = rcpFromRef (X);
486 Y_local = rcpFromRef (Y);
487 }
488
489 // Apply the local operator:
490 // Y_local := beta*Y_local + alpha*M^{-1}*X_local
491 this->applyImpl (*X_local, *Y_local, mode, alpha, beta);
492
493 ++numApply_; // We've successfully finished the work of apply().
494 }
495
496 applyTime_ += (timer->wallTime() - startTime);
497}
498
499
500template<class MatrixType>
502{
503 std::ostringstream out;
504
505 // Output is a valid YAML dictionary in flow style. If you don't
506 // like everything on a single line, you should call describe()
507 // instead.
508 out << "\"Ifpack2::Details::TriDiSolver\": ";
509 out << "{";
510 if (this->getObjectLabel () != "") {
511 out << "Label: \"" << this->getObjectLabel () << "\", ";
512 }
513 out << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
514 << "Computed: " << (isComputed () ? "true" : "false") << ", ";
515
516 if (A_.is_null ()) {
517 out << "Matrix: null";
518 }
519 else {
520 out << "Matrix: not null"
521 << ", Global matrix dimensions: ["
522 << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]";
523 }
524
525 out << "}";
526 return out.str ();
527}
528
529
530template<class MatrixType>
531void TriDiSolver<MatrixType, false>::describeLocal (Teuchos::FancyOStream& out,
532 const Teuchos::EVerbosityLevel verbLevel) const {
533 using Teuchos::FancyOStream;
534 using Teuchos::OSTab;
535 using Teuchos::RCP;
536 using Teuchos::rcpFromRef;
537 using std::endl;
538
539 if (verbLevel == Teuchos::VERB_NONE) {
540 return;
541 }
542 else {
543 RCP<FancyOStream> ptrOut = rcpFromRef (out);
544 OSTab tab1 (ptrOut);
545 if (this->getObjectLabel () != "") {
546 out << "label: " << this->getObjectLabel () << endl;
547 }
548 out << "initialized: " << (isInitialized_ ? "true" : "false") << endl
549 << "computed: " << (isComputed_ ? "true" : "false") << endl
550 << "number of initialize calls: " << numInitialize_ << endl
551 << "number of compute calls: " << numCompute_ << endl
552 << "number of apply calls: " << numApply_ << endl
553 << "total time in seconds in initialize: " << initializeTime_ << endl
554 << "total time in seconds in compute: " << computeTime_ << endl
555 << "total time in seconds in apply: " << applyTime_ << endl;
556 if (verbLevel >= Teuchos::VERB_EXTREME) {
557 out << "A_local_tridi_:" << endl;
558 A_local_tridi_.print(out);
559 }
560 out << "ipiv_: " << Teuchos::toString (ipiv_) << endl;
561 }
562}
563
564template<class MatrixType>
565void TriDiSolver<MatrixType, false>::describe (Teuchos::FancyOStream& out,
566 const Teuchos::EVerbosityLevel verbLevel) const
567{
568 using Teuchos::FancyOStream;
569 using Teuchos::OSTab;
570 using Teuchos::RCP;
571 using Teuchos::rcpFromRef;
572 using std::endl;
573
574 RCP<FancyOStream> ptrOut = rcpFromRef (out);
575 OSTab tab0 (ptrOut);
576 if (A_.is_null ()) {
577 // If A_ is null, we don't have a communicator, so we can't
578 // safely print local data on all processes. Just print the
579 // local data without arbitration between processes, and hope
580 // for the best.
581 if (verbLevel > Teuchos::VERB_NONE) {
582 out << "Ifpack2::Details::TriDiSolver:" << endl;
583 }
584 describeLocal (out, verbLevel);
585 }
586 else {
587 // If A_ is not null, we have a communicator, so we can
588 // arbitrate among all processes to print local data.
589 const Teuchos::Comm<int>& comm = * (A_->getRowMap ()->getComm ());
590 const int myRank = comm.getRank ();
591 const int numProcs = comm.getSize ();
592 if (verbLevel > Teuchos::VERB_NONE && myRank == 0) {
593 out << "Ifpack2::Details::TriDiSolver:" << endl;
594 }
595 OSTab tab1 (ptrOut);
596 for (int p = 0; p < numProcs; ++p) {
597 if (myRank == p) {
598 out << "Process " << myRank << ":" << endl;
599 describeLocal (out, verbLevel);
600 }
601 comm.barrier ();
602 comm.barrier ();
603 comm.barrier ();
604 } // for p = 0 .. numProcs-1
605 }
606}
607
608template<class MatrixType>
609void TriDiSolver<MatrixType, false>::extract (Teuchos::SerialTriDiMatrix<int, scalar_type>& A_local_tridi,
610 const row_matrix_type& A_local)
611{
612 using Teuchos::Array;
613 using Teuchos::ArrayView;
614 typedef local_ordinal_type LO;
615 typedef typename Teuchos::ArrayView<LO>::size_type size_type;
616
617 // Fill the local tridi matrix with zeros.
618 A_local_tridi.putScalar (STS::zero ());
619
620 //
621 // Map both row and column indices to local indices. We can use the
622 // row Map's local indices for row indices, and the column Map's
623 // local indices for column indices. It doesn't really matter;
624 // either way is just a permutation of rows and columns.
625 //
626 const map_type& rowMap = * (A_local.getRowMap ());
627
628 // Temporary arrays to hold the indices and values of the entries in
629 // each row of A_local.
630 const size_type maxNumRowEntries =
631 static_cast<size_type> (A_local.getLocalMaxNumRowEntries ());
632 nonconst_local_inds_host_view_type localIndices("localIndices",maxNumRowEntries);
633 nonconst_values_host_view_type values ("values",maxNumRowEntries);
634
635 const LO numLocalRows = static_cast<LO> (rowMap.getLocalNumElements ());
636 const LO minLocalRow = rowMap.getMinLocalIndex ();
637 // This slight complication of computing the upper loop bound avoids
638 // issues if the row Map has zero entries on the calling process.
639 const LO maxLocalRow = minLocalRow + numLocalRows; // exclusive bound
640 for (LO localRow = minLocalRow; localRow < maxLocalRow; ++localRow) {
641 // The LocalFilter automatically excludes "off-process" entries.
642 // That means all the column indices in this row belong to the
643 // domain Map. We can, therefore, just use the local row and
644 // column indices to put each entry directly in the tridi matrix.
645 // It's OK if the column Map puts the local indices in a different
646 // order; the Import will bring them into the correct order.
647 const size_type numEntriesInRow =
648 static_cast<size_type> (A_local.getNumEntriesInLocalRow (localRow));
649 size_t numEntriesOut = 0; // ignored
650 A_local.getLocalRowCopy (localRow,
651 localIndices,
652 values,
653 numEntriesOut);
654 for (LO k = 0; k < numEntriesInRow; ++k) {
655 const LO localCol = localIndices[k];
656 const scalar_type val = values[k];
657 // We use += instead of =, in case there are duplicate entries
658 // in the row. There should not be, but why not be general?
659 // NOTE: we only extract the TriDi part of the row matrix. Do not extract DU2
660 if( localCol >= localRow-1 && localCol <= localRow+1 )
661 A_local_tridi(localRow, localCol) += val;
662 }
663 }
664}
665
667// Stub implementation
669
670template<class MatrixType>
671TriDiSolver<MatrixType, true>::TriDiSolver (const Teuchos::RCP<const row_matrix_type>& A) {
672 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
673}
674
675
676template<class MatrixType>
677Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::map_type>
679 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
680}
681
682
683template<class MatrixType>
684Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::map_type>
686 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
687}
688
689
690template<class MatrixType>
691void
692TriDiSolver<MatrixType, true>::setParameters (const Teuchos::ParameterList& params) {
693 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
694}
695
696
697template<class MatrixType>
698bool
700 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
701}
702
703
704template<class MatrixType>
705bool
707 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
708}
709
710
711template<class MatrixType>
712int
714 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
715}
716
717
718template<class MatrixType>
719int
721 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
722}
723
724
725template<class MatrixType>
726int
728 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
729}
730
731
732template<class MatrixType>
733double
735 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
736}
737
738
739template<class MatrixType>
740double
742 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
743}
744
745
746template<class MatrixType>
747double
749 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
750}
751
752
753template<class MatrixType>
754Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::row_matrix_type>
756 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
757}
758
759
760template<class MatrixType>
761void TriDiSolver<MatrixType, true>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
762{
763 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
764}
765
766
767template<class MatrixType>
769{
770 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
771}
772
773
774template<class MatrixType>
775TriDiSolver<MatrixType, true>::~TriDiSolver ()
776{
777 // Destructors should never throw exceptions.
778}
779
780
781template<class MatrixType>
783{
784 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
785}
786
787
788template<class MatrixType>
790 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
791 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
792 Teuchos::ETransp mode,
793 scalar_type alpha,
794 scalar_type beta) const
795{
796 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
797}
798
799
800template<class MatrixType>
801std::string
802TriDiSolver<MatrixType, true>::description () const
803{
804 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
805}
806
807
808template<class MatrixType>
809void TriDiSolver<MatrixType, true>::describe(Teuchos::FancyOStream& out,
810 const Teuchos::EVerbosityLevel verbLevel) const
811{
812 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented");
813}
814
815}// namespace Details
816} // namespace Ifpack2
817
818#define IFPACK2_DETAILS_TRIDISOLVER_INSTANT(S,LO,GO,N) \
819 template class Ifpack2::Details::TriDiSolver< Tpetra::RowMatrix<S, LO, GO, N> >;
820
821#endif // IFPACK2_DETAILS_TRIDISOLVER_HPP
virtual void setMatrix(const Teuchos::RCP< const Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > &A)=0
Set the new matrix.
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition Ifpack2_Details_TriDiSolver_decl.hpp:109
"Preconditioner" that uses LAPACK's tridi LU.
Definition Ifpack2_Details_TriDiSolver_decl.hpp:85
Access only local rows and columns of a sparse matrix.
Definition Ifpack2_LocalFilter_decl.hpp:163
virtual bool isInitialized() const=0
True if the preconditioner has been successfully initialized, else false.
virtual Teuchos::RCP< const Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > getMatrix() const=0
The input matrix given to the constructor.
virtual void apply(const Tpetra::MultiVector< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > &X, Tpetra::MultiVector< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, MatrixType::scalar_type alpha=Teuchos::ScalarTraits< MatrixType::scalar_type >::one(), MatrixType::scalar_type beta=Teuchos::ScalarTraits< MatrixType::scalar_type >::zero()) const=0
Apply the preconditioner to X, putting the result in Y.
virtual Teuchos::RCP< const Tpetra::Map< MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > getRangeMap() const=0
The range Map of this operator.
virtual Teuchos::RCP< const Tpetra::Map< MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > getDomainMap() const=0
The domain Map of this operator.
virtual bool isComputed() const=0
True if the preconditioner has been successfully computed, else false.
TRANS
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74