/[escript]/branches/trilinos_from_5897/trilinoswrap/src/Amesos2Wrapper.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/trilinoswrap/src/Amesos2Wrapper.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6104 - (show annotations)
Wed Mar 30 06:01:20 2016 UTC (3 years ago) by caltinay
File size: 3419 byte(s)
Factored out and templetized preconditioner,solver and direct solver creation.
The SystemMatrix constructor now takes an optional arg 'isComplex'.
Some complex operations are commented out as we need the complex getSampleData*
methods from trunk for them to work.

It looks like we have to modify the Abstract class in escript eventually as
there is a single method that takes a `double` argument (nullifyRowsAndCols).


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include <trilinoswrap/Amesos2Wrapper.h>
18 #include <trilinoswrap/TrilinosAdapterException.h>
19
20 #include <escript/SolverOptions.h>
21
22 #include <Amesos2.hpp>
23
24 using Teuchos::RCP;
25
26 namespace esys_trilinos {
27
28 template<typename ST>
29 RCP<DirectSolverType<ST> > createDirectSolver(const escript::SolverBuddy& sb,
30 RCP<const MatrixType<ST> > A,
31 RCP<VectorType<ST> > X,
32 RCP<const VectorType<ST> > B)
33 {
34 typedef MatrixType<ST> MT;
35 typedef VectorType<ST> VT;
36
37 RCP<DirectSolverType<ST> > solver;
38 RCP<Teuchos::ParameterList> params = Teuchos::parameterList();
39
40 // TODO: options!
41 if (Amesos2::query("MUMPS")) {
42 solver = Amesos2::create<MT, VT>("MUMPS", A, X, B);
43 if (sb.isVerbose()) {
44 params->set("ICNTL(4)", 4);
45 }
46 } else if (Amesos2::query("klu2")) {
47 solver = Amesos2::create<MT, VT>("klu2", A, X, B);
48 params->set("DiagPivotThresh", sb.getDiagonalDominanceThreshold());
49 params->set("SymmetricMode", sb.isSymmetric());
50 } else if (Amesos2::query("superludist")) {
51 solver = Amesos2::create<MT, VT>("superludist", A, X, B);
52 } else if (Amesos2::query("superlu")) {
53 solver = Amesos2::create<MT, VT>("superlu", A, X, B);
54 params->set("DiagPivotThresh", sb.getDiagonalDominanceThreshold());
55 params->set("ILU_DropTol", sb.getDropTolerance());
56 params->set("SymmetricMode", sb.isSymmetric());
57 } else if (Amesos2::query("pardiso_mkl")) {
58 solver = Amesos2::create<MT, VT>("pardiso_mkl", A, X, B);
59 } else if (Amesos2::query("lapack")) {
60 solver = Amesos2::create<MT, VT>("lapack", A, X, B);
61 } else if (Amesos2::query("amesos2_cholmod")) {
62 solver = Amesos2::create<MT, VT>("amesos2_cholmod", A, X, B);
63 } else if (Amesos2::query("superlumt")) {
64 solver = Amesos2::create<MT, VT>("superlumt", A, X, B);
65 params->set("nprocs", omp_get_max_threads());
66 params->set("DiagPivotThresh", sb.getDiagonalDominanceThreshold());
67 params->set("SymmetricMode", sb.isSymmetric());
68 } else {
69 throw TrilinosAdapterException("Could not find an Amesos2 direct solver!");
70 }
71 solver->setParameters(params);
72 return solver;
73 }
74
75 // instantiate our two supported versions
76 template
77 RCP<DirectSolverType<real_t> > createDirectSolver<real_t>(
78 const escript::SolverBuddy& sb, RCP<const RealMatrix> A,
79 RCP<RealVector> X, RCP<const RealVector> B);
80 template
81 RCP<DirectSolverType<cplx_t> > createDirectSolver<cplx_t>(
82 const escript::SolverBuddy& sb, RCP<const ComplexMatrix> A,
83 RCP<ComplexVector> X, RCP<const ComplexVector> B);
84
85 } // end of namespace
86

  ViewVC Help
Powered by ViewVC 1.1.26