/[escript]/trunk/trilinoswrap/src/Amesos2Wrapper.cpp
ViewVC logotype

Contents of /trunk/trilinoswrap/src/Amesos2Wrapper.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6799 - (show annotations)
Mon Mar 25 05:53:58 2019 UTC (3 weeks, 5 days ago) by aellery
File size: 10836 byte(s)
I have rewritten the solverbuddy. Briefly:

1. The remaining AMG code has been removed from PASO.
2. If Trilinos is available, escript will now use it by default.
3. eScript will use a direct solver by default, (if one is available,) when solving 2 dimensional problems and an iterative solver, by default, when solving 3 dimensional problems. This can be changed by a user by manually specifying which solver to use.
4. There is a new option available, setHermitian(), that allows a user to specify when a coefficient matrix is Hermitian.
5. Symmetry information is always passed onto the Trilinos solver when this information is relevant.
6. All tests have been updated, when relevant, to reflect these changes.
7. I fixed a couple of undocumented bugs.


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 Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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 #include <trilinoswrap/util.h>
20
21 #include <escript/SolverOptions.h>
22
23 #include <Amesos2.hpp>
24 #include <Tpetra_CrsMatrix.hpp>
25
26 #include <boost/python/dict.hpp>
27
28 using Teuchos::RCP;
29
30 namespace bp = boost::python;
31
32 namespace esys_trilinos {
33
34 template<class Matrix, class Vector>
35 RCP<DirectSolverType<Matrix,Vector> > createDirectSolver(
36 const escript::SolverBuddy& sb,
37 RCP<const Matrix> A,
38 RCP<Vector> X,
39 RCP<const Vector> B)
40 {
41 using util::extractParamIfSet;
42
43 typedef typename Matrix::scalar_type ST;
44
45 RCP<DirectSolverType<Matrix,Vector> > solver;
46 RCP<Teuchos::ParameterList> amesosParams = Teuchos::parameterList("Amesos2");
47 const bp::dict& pyParams = sb.getTrilinosParameters();
48
49 const escript::SolverOptions method = sb.getSolverMethod();
50 // did user request specific direct solver or not?
51 const bool dontcare = method == escript::SO_METHOD_DIRECT;
52
53 if ((dontcare || method == escript::SO_METHOD_DIRECT_TRILINOS) &&
54 Amesos2::query("klu2")) {
55 solver = Amesos2::create<Matrix, Vector>("klu2", A, X, B);
56 Teuchos::ParameterList solverParams(solver->name());
57 // the doco says these params exist but clearly they don't :-(
58 extractParamIfSet<std::string>("Trans", pyParams, solverParams);
59 extractParamIfSet<bool>("Equil", pyParams, solverParams);
60 extractParamIfSet<std::string>("IterRefine", pyParams, solverParams);
61 extractParamIfSet<bool>("SymmetricMode", pyParams, solverParams);
62 extractParamIfSet<ST>("DiagPivotThresh", pyParams, solverParams);
63 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
64 amesosParams->set(solver->name(), solverParams);
65 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_MUMPS) &&
66 Amesos2::query("MUMPS")) {
67 solver = Amesos2::create<Matrix, Vector>("MUMPS", A, X, B);
68 Teuchos::ParameterList solverParams(solver->name());
69 // solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
70 if (sb.isVerbose()) {
71 solverParams.set("ICNTL(4)", 4);
72 }
73 extractParamIfSet<int>("ICNTL(1)", pyParams, solverParams);
74 extractParamIfSet<int>("ICNTL(2)", pyParams, solverParams);
75 extractParamIfSet<int>("ICNTL(3)", pyParams, solverParams);
76 extractParamIfSet<int>("ICNTL(4)", pyParams, solverParams);
77 extractParamIfSet<int>("ICNTL(6)", pyParams, solverParams);
78 extractParamIfSet<int>("ICNTL(9)", pyParams, solverParams);
79 extractParamIfSet<int>("ICNTL(11)", pyParams, solverParams);
80 amesosParams->set(solver->name(), solverParams);
81 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_TRILINOS) &&
82 Amesos2::query("Basker")) {
83 solver = Amesos2::create<Matrix, Vector>("Basker", A, X, B);
84 Teuchos::ParameterList solverParams(solver->name());
85 solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
86 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_SUPERLU) &&
87 Amesos2::query("superludist")) {
88 solver = Amesos2::create<Matrix, Vector>("superludist", A, X, B);
89 Teuchos::ParameterList solverParams(solver->name());
90 solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
91 extractParamIfSet<int>("npcol", pyParams, solverParams);
92 extractParamIfSet<int>("nprow", pyParams, solverParams);
93 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
94 extractParamIfSet<bool>("ReplaceTinyPivot", pyParams, solverParams);
95 amesosParams->set(solver->name(), solverParams);
96 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_SUPERLU) &&
97 Amesos2::query("superlu")) {
98 solver = Amesos2::create<Matrix, Vector>("superlu", A, X, B);
99 Teuchos::ParameterList solverParams(solver->name());
100 solverParams.set("SymmetricMode", sb.isSymmetric());
101 extractParamIfSet<std::string>("Trans", pyParams, solverParams);
102 extractParamIfSet<bool>("Equil", pyParams, solverParams);
103 extractParamIfSet<std::string>("IterRefine", pyParams, solverParams);
104 extractParamIfSet<bool>("SymmetricMode", pyParams, solverParams);
105 extractParamIfSet<ST>("DiagPivotThresh", pyParams, solverParams);
106 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
107 extractParamIfSet<bool>("ILU_Flag", pyParams, solverParams);
108 extractParamIfSet<ST>("ILU_DropTol", pyParams, solverParams);
109 extractParamIfSet<ST>("ILU_FillFactor", pyParams, solverParams);
110 extractParamIfSet<std::string>("ILU_Norm", pyParams, solverParams);
111 extractParamIfSet<std::string>("ILU_MILU", pyParams, solverParams);
112 extractParamIfSet<ST>("ILU_FillTol", pyParams, solverParams);
113 amesosParams->set(solver->name(), solverParams);
114 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_SUPERLU) &&
115 Amesos2::query("superlumt")) {
116 solver = Amesos2::create<Matrix, Vector>("superlumt", A, X, B);
117 Teuchos::ParameterList solverParams(solver->name());
118 #ifdef _OPENMP
119 solverParams.set("nprocs", omp_get_max_threads());
120 #else
121 solverParams.set("nprocs", 1);
122 #endif
123 solverParams.set("SymmetricMode", sb.isSymmetric());
124 extractParamIfSet<int>("nprocs", pyParams, solverParams);
125 extractParamIfSet<std::string>("trans", pyParams, solverParams);
126 extractParamIfSet<int>("panel_size", pyParams, solverParams);
127 extractParamIfSet<int>("relax", pyParams, solverParams);
128 extractParamIfSet<bool>("Equil", pyParams, solverParams);
129 extractParamIfSet<bool>("SymmetricMode", pyParams, solverParams);
130 extractParamIfSet<ST>("DiagPivotThresh", pyParams, solverParams);
131 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
132 amesosParams->set(solver->name(), solverParams);
133 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_PARDISO) &&
134 Amesos2::query("pardiso_mkl")) {
135 solver = Amesos2::create<Matrix, Vector>("pardiso_mkl", A, X, B);
136 Teuchos::ParameterList solverParams(solver->name());
137 extractParamIfSet<int>("IPARM(2)", pyParams, solverParams);
138 extractParamIfSet<int>("IPARM(4)", pyParams, solverParams);
139 extractParamIfSet<int>("IPARM(8)", pyParams, solverParams);
140 extractParamIfSet<int>("IPARM(10)", pyParams, solverParams);
141 extractParamIfSet<int>("IPARM(18)", pyParams, solverParams);
142 extractParamIfSet<int>("IPARM(24)", pyParams, solverParams);
143 extractParamIfSet<int>("IPARM(25)", pyParams, solverParams);
144 extractParamIfSet<int>("IPARM(60)", pyParams, solverParams);
145 amesosParams->set(solver->name(), solverParams);
146 } else if (Amesos2::query("amesos2_cholmod")) {
147 solver = Amesos2::create<Matrix, Vector>("amesos2_cholmod", A, X, B);
148 Teuchos::ParameterList solverParams(solver->name());
149
150 solverParams.set("SymmetricMode", sb.isSymmetric());
151 extractParamIfSet<std::string>("Trans", pyParams, solverParams);
152 extractParamIfSet<bool>("Equil", pyParams, solverParams);
153 extractParamIfSet<std::string>("IterRefine", pyParams, solverParams);
154 extractParamIfSet<bool>("SymmetricMode", pyParams, solverParams);
155 extractParamIfSet<ST>("DiagPivotThresh", pyParams, solverParams);
156 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
157 amesosParams->set(solver->name(), solverParams);
158 } else if (Amesos2::query("lapack")) {
159 solver = Amesos2::create<Matrix, Vector>("lapack", A, X, B);
160 Teuchos::ParameterList solverParams(solver->name());
161 solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
162 } else {
163 if (dontcare) {
164 throw TrilinosAdapterException("Could not find an Amesos2 direct solver!");
165 } else {
166 throw TrilinosAdapterException("The requested direct solver is not available!");
167 }
168 }
169 solver->setParameters(amesosParams);
170 return solver;
171 }
172
173 typedef Tpetra::CrsMatrix<real_t,LO,GO,NT> RealMatrix;
174 typedef Tpetra::CrsMatrix<cplx_t,LO,GO,NT> ComplexMatrix;
175
176 // instantiate
177 template
178 RCP<DirectSolverType<RealMatrix, RealVector> >
179 createDirectSolver<RealMatrix,RealVector>(const escript::SolverBuddy& sb,
180 RCP<const RealMatrix> A,
181 RCP<RealVector> X,
182 RCP<const RealVector> B);
183 template
184 RCP<DirectSolverType<ComplexMatrix, ComplexVector> >
185 createDirectSolver<ComplexMatrix, ComplexVector>(
186 const escript::SolverBuddy& sb,
187 RCP<const ComplexMatrix> A,
188 RCP<ComplexVector> X,
189 RCP<const ComplexVector> B);
190
191 /* Amesos2 does not currently support block matrices!
192 template
193 RCP<DirectSolverType<RealBlockMatrix, RealBlockVector> >
194 createDirectSolver<RealBlockMatrix,RealBlockVector>(
195 const escript::SolverBuddy& sb,
196 RCP<const RealBlockMatrix> A,
197 RCP<RealBlockVector> X,
198 RCP<const RealBlockVector> B);
199 template
200 RCP<DirectSolverType<ComplexBlockMatrix, ComplexBlockVector> >
201 createDirectSolver<ComplexBlockMatrix, ComplexBlockVector>(
202 const escript::SolverBuddy& sb,
203 RCP<const ComplexBlockMatrix> A,
204 RCP<ComplexBlockVector> X,
205 RCP<const ComplexBlockVector> B);
206 */
207
208 } // end of namespace
209

  ViewVC Help
Powered by ViewVC 1.1.26