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