Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_SolverCore_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Amesos2: Templated Direct Sparse Solver Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//
42// @HEADER
43
53#ifndef AMESOS2_SOLVERCORE_DEF_HPP
54#define AMESOS2_SOLVERCORE_DEF_HPP
55
56#include "Kokkos_ArithTraits.hpp"
57
58#include "Amesos2_MatrixAdapter_def.hpp"
59#include "Amesos2_MultiVecAdapter_def.hpp"
60
61#include "Amesos2_Util.hpp"
62
63#include "KokkosSparse_spmv.hpp"
64#include "KokkosBlas.hpp"
65
66namespace Amesos2 {
67
68
69template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
71 Teuchos::RCP<const Matrix> A,
72 Teuchos::RCP<Vector> X,
73 Teuchos::RCP<const Vector> B )
74 : matrixA_(createConstMatrixAdapter<Matrix>(A))
75 , multiVecX_(X) // may be null
76 , multiVecB_(B) // may be null
77 , globalNumRows_(matrixA_->getGlobalNumRows())
78 , globalNumCols_(matrixA_->getGlobalNumCols())
79 , globalNumNonZeros_(matrixA_->getGlobalNNZ())
80 , rowIndexBase_(matrixA_->getRowIndexBase())
81 , columnIndexBase_(matrixA_->getColumnIndexBase())
82 , rank_(Teuchos::rank(*this->getComm()))
83 , root_(rank_ == 0)
84 , nprocs_(Teuchos::size(*this->getComm()))
85{
86 TEUCHOS_TEST_FOR_EXCEPTION(
88 std::invalid_argument,
89 "Matrix shape inappropriate for this solver");
90}
91
92
94template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
99
100
101template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
102Solver<Matrix,Vector>&
104{
105#ifdef HAVE_AMESOS2_TIMERS
106 Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
107#endif
108
109 loadA(PREORDERING);
110
111 static_cast<solver_type*>(this)->preOrdering_impl();
112 ++status_.numPreOrder_;
113 status_.last_phase_ = PREORDERING;
114
115 return *this;
116}
117
118
119template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
120Solver<Matrix,Vector>&
122{
123#ifdef HAVE_AMESOS2_TIMERS
124 Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
125#endif
126
127 if( !status_.preOrderingDone() ){
128 preOrdering();
129 if( !matrix_loaded_ ) loadA(SYMBFACT);
130 } else {
131 loadA(SYMBFACT);
132 }
133
134 static_cast<solver_type*>(this)->symbolicFactorization_impl();
135 ++status_.numSymbolicFact_;
136 status_.last_phase_ = SYMBFACT;
137
138 return *this;
139}
141
142template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
143Solver<Matrix,Vector>&
145{
146#ifdef HAVE_AMESOS2_TIMERS
147 Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
148#endif
149
150 if( !status_.symbolicFactorizationDone() ){
151 symbolicFactorization();
152 if( !matrix_loaded_ ) loadA(NUMFACT);
153 } else {
154 loadA(NUMFACT);
155 }
156
157 static_cast<solver_type*>(this)->numericFactorization_impl();
158 ++status_.numNumericFact_;
159 status_.last_phase_ = NUMFACT;
160
161 return *this;
162}
163
164
165template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
166void
168{
169 solve(multiVecX_.ptr(), multiVecB_.ptr());
170}
171
172template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
173void
175 const Teuchos::Ptr<const Vector> B) const
176{
177#ifdef HAVE_AMESOS2_TIMERS
178 Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
179#endif
180
181 X.assert_not_null();
182 B.assert_not_null();
183
184 if (control_.useIterRefine_) {
185 solve_ir(X, B, control_.maxNumIterRefines_, control_.verboseIterRefine_);
186 return;
187 }
188
189 const Teuchos::RCP<MultiVecAdapter<Vector> > x =
190 createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(X));
191 const Teuchos::RCP<const MultiVecAdapter<Vector> > b =
192 createConstMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(B));
193
194#ifdef HAVE_AMESOS2_DEBUG
195 // Check some required properties of X and B
196 TEUCHOS_TEST_FOR_EXCEPTION
197 (x->getGlobalLength() != matrixA_->getGlobalNumCols(),
198 std::invalid_argument,
199 "MultiVector X must have length equal to the number of "
200 "global columns in A. X->getGlobalLength() = "
201 << x->getGlobalLength() << " != A->getGlobalNumCols() = "
202 << matrixA_->getGlobalNumCols() << ".");
203
204 TEUCHOS_TEST_FOR_EXCEPTION(b->getGlobalLength() != matrixA_->getGlobalNumRows(),
205 std::invalid_argument,
206 "MultiVector B must have length equal to the number of "
207 "global rows in A");
208
209 TEUCHOS_TEST_FOR_EXCEPTION(x->getGlobalNumVectors() != b->getGlobalNumVectors(),
210 std::invalid_argument,
211 "X and B MultiVectors must have the same number of vectors");
212#endif // HAVE_AMESOS2_DEBUG
213
214 if( !status_.numericFactorizationDone() ){
215 // This casting-away of constness is probably OK because this
216 // function is meant to be "logically const"
217 const_cast<type&>(*this).numericFactorization();
218 }
219
220 static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*x), Teuchos::ptrInArg(*b));
221 ++status_.numSolve_;
222 status_.last_phase_ = SOLVE;
223}
224
225template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
226void
227SolverCore<ConcreteSolver,Matrix,Vector>::solve(Vector* X, const Vector* B) const
228{
229 solve(Teuchos::ptr(X), Teuchos::ptr(B));
230}
231
232
233template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
234int
235SolverCore<ConcreteSolver,Matrix,Vector>::solve_ir(const int maxNumIters, const bool verbose)
237 return solve_ir(multiVecX_.ptr(), multiVecB_.ptr(), maxNumIters, verbose);
238}
239
240template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
241int
242SolverCore<ConcreteSolver,Matrix,Vector>::solve_ir(Vector* X, const Vector* B, const int maxNumIters, const bool verbose) const
243{
244 return solve_ir(Teuchos::ptr(X), Teuchos::ptr(B), maxNumIters, verbose);
245}
246
247template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
248int
249SolverCore<ConcreteSolver,Matrix,Vector>::solve_ir(const Teuchos::Ptr< Vector> x,
250 const Teuchos::Ptr<const Vector> b,
251 const int maxNumIters,
252 const bool verbose) const
253{
254 using KAT = Kokkos::ArithTraits<scalar_type>;
255 using impl_scalar_type = typename KAT::val_type;
256 using magni_type = typename KAT::mag_type;
257 using host_execution_space = Kokkos::DefaultHostExecutionSpace;
258 using host_crsmat_t = KokkosSparse::CrsMatrix<impl_scalar_type, int, host_execution_space, void, int>;
259 using host_graph_t = typename host_crsmat_t::StaticCrsGraphType;
260 using host_values_t = typename host_crsmat_t::values_type::non_const_type;
261 using host_row_map_t = typename host_graph_t::row_map_type::non_const_type;
262 using host_colinds_t = typename host_graph_t::entries_type::non_const_type;
263 using host_mvector_t = Kokkos::View<impl_scalar_type **, Kokkos::LayoutLeft, host_execution_space>;
264 using host_vector_t = Kokkos::View<impl_scalar_type *, Kokkos::LayoutLeft, host_execution_space>;
265 using host_magni_view = Kokkos::View<magni_type *, Kokkos::LayoutLeft, host_execution_space>;
266
267 const impl_scalar_type one(1.0);
268 const impl_scalar_type mone = impl_scalar_type(-one);
269 const magni_type eps = KAT::eps ();
270
271 // get data needed for IR
272 using MVAdapter = MultiVecAdapter<Vector>;
273 Teuchos::RCP< MVAdapter> X = createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(x));
274 Teuchos::RCP<const MVAdapter> B = createConstMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(b));
275
276 auto r_ = B->clone();
277 auto e_ = X->clone();
278 auto r = r_.ptr();
279 auto e = e_.ptr();
280 Teuchos::RCP< MVAdapter> R = createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(r));
281 Teuchos::RCP< MVAdapter> E = createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(e));
282
283 const size_t nrhs = X->getGlobalNumVectors();
284 const int nnz = this->matrixA_->getGlobalNNZ();
285 const int nrows = this->matrixA_->getGlobalNumRows();
286
287 // get local matriix
288 host_crsmat_t crsmat;
289 host_graph_t static_graph;
290 host_row_map_t rowmap_view;
291 host_colinds_t colind_view;
292 host_values_t values_view;
293 if (this->root_) {
294 Kokkos::resize(rowmap_view, 1+nrows);
295 Kokkos::resize(colind_view, nnz);
296 Kokkos::resize(values_view, nnz);
297 } else {
298 Kokkos::resize(rowmap_view, 1);
299 Kokkos::resize(colind_view, 0);
300 Kokkos::resize(values_view, 0);
301 }
302
303 int nnz_ret = 0;
305 MatrixAdapter<Matrix>, host_values_t, host_colinds_t, host_row_map_t>::do_get(
306 this->matrixA_.ptr(),
307 values_view, colind_view, rowmap_view,
308 nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_);
309
310 if (this->root_) {
311 static_graph = host_graph_t(colind_view, rowmap_view);
312 crsmat = host_crsmat_t("CrsMatrix", nrows, values_view, static_graph);
313 }
314
315 //
316 // ** First Solve **
317 static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*X), Teuchos::ptrInArg(*B));
318
319
320 // auxiliary scalar Kokkos views
321 const int ldx = (this->root_ ? X->getGlobalLength() : 0);
322 const int ldb = (this->root_ ? B->getGlobalLength() : 0);
323 const int ldr = (this->root_ ? R->getGlobalLength() : 0);
324 const int lde = (this->root_ ? E->getGlobalLength() : 0);
325 const bool initialize_data = true;
326 const bool not_initialize_data = true;
327 host_mvector_t X_view;
328 host_mvector_t B_view;
329 host_mvector_t R_view;
330 host_mvector_t E_view;
331
332 global_size_type rowIndexBase = this->rowIndexBase_;
333 auto Xptr = Teuchos::Ptr< MVAdapter>(X.ptr());
334 auto Bptr = Teuchos::Ptr<const MVAdapter>(B.ptr());
335 auto Rptr = Teuchos::Ptr< MVAdapter>(R.ptr());
336 auto Eptr = Teuchos::Ptr< MVAdapter>(E.ptr());
337 Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
338 do_get( initialize_data, Xptr, X_view, ldx, CONTIGUOUS_AND_ROOTED, rowIndexBase);
339 Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
340 do_get( initialize_data, Bptr, B_view, ldb, CONTIGUOUS_AND_ROOTED, rowIndexBase);
341 Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
342 do_get(not_initialize_data, Rptr, R_view, ldr, CONTIGUOUS_AND_ROOTED, rowIndexBase);
343 Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
344 do_get(not_initialize_data, Eptr, E_view, lde, CONTIGUOUS_AND_ROOTED, rowIndexBase);
345
346
347 host_magni_view x0norms("x0norms", nrhs);
348 host_magni_view bnorms("bnorms", nrhs);
349 host_magni_view enorms("enorms", nrhs);
350 if (this->root_) {
351 // compute initial solution norms (used for stopping criteria)
352 for (size_t j = 0; j < nrhs; j++) {
353 auto x_subview = Kokkos::subview(X_view, Kokkos::ALL(), j);
354 host_vector_t x_1d (const_cast<impl_scalar_type*>(x_subview.data()), x_subview.extent(0));
355 x0norms(j) = KokkosBlas::nrm2(x_1d);
356 }
357 if (verbose) {
358 std::cout << std::endl
359 << " SolverCore :: solve_ir (maxNumIters = " << maxNumIters
360 << ", tol = " << x0norms(0) << " * " << eps << " = " << x0norms(0)*eps
361 << ")" << std::endl;
362 }
363
364 // compute residual norm
365 if (verbose) {
366 std::cout << " bnorm = ";
367 for (size_t j = 0; j < nrhs; j++) {
368 auto b_subview = Kokkos::subview(B_view, Kokkos::ALL(), j);
369 host_vector_t b_1d (const_cast<impl_scalar_type*>(b_subview.data()), b_subview.extent(0));
370 bnorms(j) = KokkosBlas::nrm2(b_1d);
371 std::cout << bnorms(j) << ", ";
372 }
373 std::cout << std::endl;
374 }
375 }
376
378 //
379 // ** Iterative Refinement **
380 int numIters = 0;
381 int converged = 0; // 0 = has not converged, 1 = converged
382 for (numIters = 0; numIters < maxNumIters && converged == 0; ++numIters) {
383 // r = b - Ax (on rank-0)
384 if (this->root_) {
385 Kokkos::deep_copy(R_view, B_view);
386 KokkosSparse::spmv("N", mone, crsmat, X_view, one, R_view);
387 Kokkos::fence();
388
389 if (verbose) {
390 // compute residual norm
391 std::cout << " > " << numIters << " : norm(r,x,e) = ";
392 for (size_t j = 0; j < nrhs; j++) {
393 auto r_subview = Kokkos::subview(R_view, Kokkos::ALL(), j);
394 auto x_subview = Kokkos::subview(X_view, Kokkos::ALL(), j);
395 host_vector_t r_1d (const_cast<impl_scalar_type*>(r_subview.data()), r_subview.extent(0));
396 host_vector_t x_1d (const_cast<impl_scalar_type*>(x_subview.data()), x_subview.extent(0));
397 impl_scalar_type rnorm = KokkosBlas::nrm2(r_1d);
398 impl_scalar_type xnorm = KokkosBlas::nrm2(x_1d);
399 std::cout << rnorm << " -> " << rnorm/bnorms(j) << " " << xnorm << " " << enorms(j) << ", ";
400 }
401 std::cout << std::endl;
403 }
404
405 // e = A^{-1} r
406 Util::put_1d_data_helper_kokkos_view<MVAdapter, host_mvector_t>::
407 do_put(Rptr, R_view, ldr, CONTIGUOUS_AND_ROOTED, rowIndexBase);
408 static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*E), Teuchos::ptrInArg(*R));
409 Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
410 do_get(initialize_data, Eptr, E_view, lde, CONTIGUOUS_AND_ROOTED, rowIndexBase);
411
412 // x = x + e (on rank-0)
413 if (this->root_) {
414 KokkosBlas::axpy(one, E_view, X_view);
416 if (numIters < maxNumIters-1) {
417 // compute norm of corrections for "convergence" check
418 converged = 1;
419 for (size_t j = 0; j < nrhs; j++) {
420 auto e_subview = Kokkos::subview(E_view, Kokkos::ALL(), j);
421 host_vector_t e_1d (const_cast<impl_scalar_type*>(e_subview.data()), e_subview.extent(0));
422 enorms(j) = KokkosBlas::nrm2(e_1d);
423 if (enorms(j) > eps * x0norms(j)) {
424 converged = 0;
426 }
427 if (verbose && converged) {
428 std::cout << " converged " << std::endl;
429 }
430 }
431 }
432
433 // broadcast "converged"
434 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &converged);
435 } // end of for-loop for IR iteration
436
437 if (verbose && this->root_) {
438 // r = b - Ax
439 Kokkos::deep_copy(R_view, B_view);
440 KokkosSparse::spmv("N", mone, crsmat, X_view, one, R_view);
441 Kokkos::fence();
442 std::cout << " > final residual norm = ";
443 for (size_t j = 0; j < nrhs; j++) {
444 auto r_subview = Kokkos::subview(R_view, Kokkos::ALL(), j);
445 host_vector_t r_1d (const_cast<impl_scalar_type*>(r_subview.data()), r_subview.extent(0));
446 scalar_type rnorm = KokkosBlas::nrm2(r_1d);
447 std::cout << rnorm << " -> " << rnorm/bnorms(j) << ", ";
448 }
449 std::cout << std::endl << std::endl;
450 }
451
452 // put X for output
453 Util::put_1d_data_helper_kokkos_view<MVAdapter, host_mvector_t>::
454 do_put(Xptr, X_view, ldx, CONTIGUOUS_AND_ROOTED, rowIndexBase);
455
456 return numIters;
457}
458
459template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
460bool
462{
463#ifdef HAVE_AMESOS2_TIMERS
464 Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
465#endif
466
467 return( static_cast<solver_type*>(this)->matrixShapeOK_impl() );
468}
469
470 // The RCP should probably be to a const Matrix, but that throws a
471 // wrench in some of the traits mechanisms that aren't expecting
472 // const types
473template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
474void
475SolverCore<ConcreteSolver,Matrix,Vector>::setA( const Teuchos::RCP<const Matrix> a,
476 EPhase keep_phase )
477{
478 matrixA_ = createConstMatrixAdapter(a);
479
480#ifdef HAVE_AMESOS2_DEBUG
481 TEUCHOS_TEST_FOR_EXCEPTION( (keep_phase != CLEAN) &&
482 (globalNumRows_ != matrixA_->getGlobalNumRows() ||
483 globalNumCols_ != matrixA_->getGlobalNumCols()),
484 std::invalid_argument,
485 "Dimensions of new matrix be the same as the old matrix if "
486 "keeping any solver phase" );
487#endif
488
489 status_.last_phase_ = keep_phase;
490
491 // Reset phase counters
492 switch( status_.last_phase_ ){
493 case CLEAN:
494 status_.numPreOrder_ = 0;
495 // Intentional fallthrough.
496 case PREORDERING:
497 status_.numSymbolicFact_ = 0;
498 // Intentional fallthrough.
499 case SYMBFACT:
500 status_.numNumericFact_ = 0;
501 // Intentional fallthrough.
502 case NUMFACT: // probably won't ever happen by itself
503 status_.numSolve_ = 0;
504 // Intentional fallthrough.
505 case SOLVE: // probably won't ever happen
506 break;
507 }
508
509 // Re-get the matrix dimensions in case they have changed
510 globalNumNonZeros_ = matrixA_->getGlobalNNZ();
511 globalNumCols_ = matrixA_->getGlobalNumCols();
512 globalNumRows_ = matrixA_->getGlobalNumRows();
513}
514
515
516template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
517Solver<Matrix,Vector>&
519 const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
520{
521#ifdef HAVE_AMESOS2_TIMERS
522 Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
523#endif
524
525 if( parameterList->name() == "Amesos2" ){
526 // First, validate the parameterList
527 Teuchos::RCP<const Teuchos::ParameterList> valid_params = getValidParameters();
528 parameterList->validateParameters(*valid_params);
529
530 // Do everything here that is for generic status and control parameters
531 control_.setControlParameters(parameterList);
532
533 // Finally, hook to the implementation's parameter list parser
534 // First check if there is a dedicated sublist for this solver and use that if there is
535 if( parameterList->isSublist(name()) ){
536 // Have control look through the solver's parameter list to see if
537 // there is anything it recognizes (mostly the "Transpose" parameter)
538 control_.setControlParameters(Teuchos::sublist(parameterList, name()));
539
540 static_cast<solver_type*>(this)->setParameters_impl(Teuchos::sublist(parameterList, name()));
541 }
542 }
543
544 return *this;
545}
546
547
548template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
549Teuchos::RCP<const Teuchos::ParameterList>
551{
552#ifdef HAVE_AMESOS2_TIMERS
553 Teuchos::TimeMonitor LocalTimer1( timers_.totalTime_ );
554#endif
555
556 using Teuchos::ParameterList;
557 using Teuchos::RCP;
558 using Teuchos::rcp;
559
560 //RCP<ParameterList> control_params = rcp(new ParameterList("Amesos2 Control"));
561 RCP<ParameterList> control_params = rcp(new ParameterList("Amesos2"));
562 control_params->set("Transpose", false, "Whether to solve with the matrix transpose");
563 control_params->set("Iterative refinement", false, "Whether to solve with iterative refinement");
564 control_params->set("Number of iterative refinements", 2, "Number of iterative refinements");
565 control_params->set("Verboes for iterative refinement", false, "Verbosity for iterative refinements");
566 // control_params->set("AddToDiag", "");
567 // control_params->set("AddZeroToDiag", false);
568 // control_params->set("MatrixProperty", "general");
569 // control_params->set("Reindex", false);
570
571 RCP<const ParameterList>
572 solver_params = static_cast<const solver_type*>(this)->getValidParameters_impl();
573 // inject the "Transpose" parameter into the solver's valid parameters
574 Teuchos::rcp_const_cast<ParameterList>(solver_params)->set("Transpose", false,
575 "Whether to solve with the "
576 "matrix transpose");
577
578 RCP<ParameterList> amesos2_params = rcp(new ParameterList("Amesos2"));
579 amesos2_params->setParameters(*control_params);
580 amesos2_params->set(name(), *solver_params);
581
582 return amesos2_params;
583}
584
585
586template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
587std::string
589{
590 std::ostringstream oss;
591 oss << name() << " solver interface";
592 return oss.str();
593}
594
595
596template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
597void
599 Teuchos::FancyOStream &out,
600 const Teuchos::EVerbosityLevel verbLevel) const
601{
602 if( matrixA_.is_null() || (rank_ != 0) ){ return; }
603 using std::endl;
604 using std::setw;
605 using Teuchos::VERB_DEFAULT;
606 using Teuchos::VERB_NONE;
607 using Teuchos::VERB_LOW;
608 using Teuchos::VERB_MEDIUM;
609 using Teuchos::VERB_HIGH;
610 using Teuchos::VERB_EXTREME;
611 Teuchos::EVerbosityLevel vl = verbLevel;
612 if (vl == VERB_DEFAULT) vl = VERB_LOW;
613 Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getComm();
614 size_t width = 1;
615 for( size_t dec = 10; dec < globalNumRows_; dec *= 10 ) {
616 ++width;
617 }
618 width = std::max<size_t>(width,size_t(11)) + 2;
619 Teuchos::OSTab tab(out);
620 // none: print nothing
621 // low: print O(1) info from node 0
622 // medium: print O(P) info, num entries per node
623 // high: print O(N) info, num entries per row
624 // extreme: print O(NNZ) info: print indices and values
625 //
626 // for medium and higher, print constituent objects at specified verbLevel
627 if( vl != VERB_NONE ) {
628 std::string p = name();
629 Util::printLine(out);
630 out << this->description() << std::endl << std::endl;
631
632 out << p << "Matrix has " << globalNumRows_ << " rows"
633 << " and " << globalNumNonZeros_ << " nonzeros"
634 << std::endl;
635 if( vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME ){
636 out << p << "Nonzero elements per row = "
637 << globalNumNonZeros_ / globalNumRows_
638 << std::endl;
639 out << p << "Percentage of nonzero elements = "
640 << 100.0 * globalNumNonZeros_ / (globalNumRows_ * globalNumCols_)
641 << std::endl;
642 }
643 if( vl == VERB_HIGH || vl == VERB_EXTREME ){
644 out << p << "Use transpose = " << control_.useTranspose_
645 << std::endl;
646 out << p << "Use iterative refinement = " << control_.useIterRefine_
647 << std::endl;
648 }
649 if ( vl == VERB_EXTREME ){
650 printTiming(out,vl);
651 }
652 Util::printLine(out);
653 }
654}
655
656
657template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
658void
660 Teuchos::FancyOStream &out,
661 const Teuchos::EVerbosityLevel verbLevel) const
662{
663 if( matrixA_.is_null() || (rank_ != 0) ){ return; }
664
665 double preTime = timers_.preOrderTime_.totalElapsedTime();
666 double symTime = timers_.symFactTime_.totalElapsedTime();
667 double numTime = timers_.numFactTime_.totalElapsedTime();
668 double solTime = timers_.solveTime_.totalElapsedTime();
669 double totTime = timers_.totalTime_.totalElapsedTime();
670 double overhead = totTime - (preTime + symTime + numTime + solTime);
671
672 std::string p = name() + " : ";
673 Util::printLine(out);
674
675 if(verbLevel != Teuchos::VERB_NONE)
676 {
677 out << p << "Time to convert matrix to implementation format = "
678 << timers_.mtxConvTime_.totalElapsedTime() << " (s)"
679 << std::endl;
680 out << p << "Time to redistribute matrix = "
681 << timers_.mtxRedistTime_.totalElapsedTime() << " (s)"
682 << std::endl;
683
684 out << p << "Time to convert vectors to implementation format = "
685 << timers_.vecConvTime_.totalElapsedTime() << " (s)"
686 << std::endl;
687 out << p << "Time to redistribute vectors = "
688 << timers_.vecRedistTime_.totalElapsedTime() << " (s)"
689 << std::endl;
690
691 out << p << "Number of pre-orderings = "
692 << status_.getNumPreOrder()
693 << std::endl;
694 out << p << "Time for pre-ordering = "
695 << preTime << " (s), avg = "
696 << preTime / status_.getNumPreOrder() << " (s)"
697 << std::endl;
698
699 out << p << "Number of symbolic factorizations = "
700 << status_.getNumSymbolicFact()
701 << std::endl;
702 out << p << "Time for sym fact = "
703 << symTime << " (s), avg = "
704 << symTime / status_.getNumSymbolicFact() << " (s)"
705 << std::endl;
706
707 out << p << "Number of numeric factorizations = "
708 << status_.getNumNumericFact()
709 << std::endl;
710 out << p << "Time for num fact = "
711 << numTime << " (s), avg = "
712 << numTime / status_.getNumNumericFact() << " (s)"
713 << std::endl;
714
715 out << p << "Number of solve phases = "
716 << status_.getNumSolve()
717 << std::endl;
718 out << p << "Time for solve = "
719 << solTime << " (s), avg = "
720 << solTime / status_.getNumSolve() << " (s)"
721 << std::endl;
722
723 out << p << "Total time spent in Amesos2 = "
724 << totTime << " (s)"
725 << std::endl;
726 out << p << "Total time spent in the Amesos2 interface = "
727 << overhead << " (s)"
728 << std::endl;
729 out << p << " (the above time does not include solver time)"
730 << std::endl;
731 out << p << "Amesos2 interface time / total time = "
732 << overhead / totTime
733 << std::endl;
734 Util::printLine(out);
735 }
736}
737
738
739template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
740void
742 Teuchos::ParameterList& timingParameterList) const
743{
744 Teuchos::ParameterList temp;
745 timingParameterList = temp.setName("NULL");
746}
747
748
749template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
750std::string
752{
753 std::string solverName = solver_type::name;
754 return solverName;
755}
756
757template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
758void
760{
761 matrix_loaded_ = static_cast<solver_type*>(this)->loadA_impl(current_phase);
762}
763
764
765} // end namespace Amesos2
766
767#endif // AMESOS2_SOLVERCORE_DEF_HPP
Utility functions for Amesos2.
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition Amesos2_SolverCore_decl.hpp:106
super_type & numericFactorization()
Performs numeric factorization on the matrix A.
Definition Amesos2_SolverCore_def.hpp:144
SolverCore(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize a Solver instance.
Definition Amesos2_SolverCore_def.hpp:70
~SolverCore()
Destructor.
Definition Amesos2_SolverCore_def.hpp:95
bool matrixShapeOK()
Returns true if the solver can handle this matrix shape.
Definition Amesos2_SolverCore_def.hpp:461
EPhase
Used to indicate a phase in the direct solution.
Definition Amesos2_TypeDecl.hpp:65
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition Amesos2_Util.hpp:663