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