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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6574 - (show annotations)
Wed May 17 00:53:31 2017 UTC (4 years, 2 months ago) by jfenwick
File size: 10736 byte(s)
Don't call openmp function if you aren't in openmp

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 //solverParams.set("DiagPivotThresh", sb.getDiagonalDominanceThreshold());
59 //solverParams.set("SymmetricMode", sb.isSymmetric());
60 extractParamIfSet<std::string>("Trans", pyParams, solverParams);
61 extractParamIfSet<bool>("Equil", pyParams, solverParams);
62 extractParamIfSet<std::string>("IterRefine", pyParams, solverParams);
63 extractParamIfSet<bool>("SymmetricMode", pyParams, solverParams);
64 extractParamIfSet<ST>("DiagPivotThresh", pyParams, solverParams);
65 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
66 amesosParams->set(solver->name(), solverParams);
67 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_MUMPS) &&
68 Amesos2::query("MUMPS")) {
69 solver = Amesos2::create<Matrix, Vector>("MUMPS", A, X, B);
70 Teuchos::ParameterList solverParams(solver->name());
71 if (sb.isVerbose()) {
72 solverParams.set("ICNTL(4)", 4);
73 }
74 extractParamIfSet<int>("ICNTL(1)", pyParams, solverParams);
75 extractParamIfSet<int>("ICNTL(2)", pyParams, solverParams);
76 extractParamIfSet<int>("ICNTL(3)", pyParams, solverParams);
77 extractParamIfSet<int>("ICNTL(4)", pyParams, solverParams);
78 extractParamIfSet<int>("ICNTL(6)", pyParams, solverParams);
79 extractParamIfSet<int>("ICNTL(9)", pyParams, solverParams);
80 extractParamIfSet<int>("ICNTL(11)", pyParams, solverParams);
81 amesosParams->set(solver->name(), solverParams);
82 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_TRILINOS) &&
83 Amesos2::query("Basker")) {
84 solver = Amesos2::create<Matrix, Vector>("Basker", A, X, B);
85 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_SUPERLU) &&
86 Amesos2::query("superludist")) {
87 solver = Amesos2::create<Matrix, Vector>("superludist", A, X, B);
88 Teuchos::ParameterList solverParams(solver->name());
89 extractParamIfSet<int>("npcol", pyParams, solverParams);
90 extractParamIfSet<int>("nprow", pyParams, solverParams);
91 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
92 extractParamIfSet<bool>("ReplaceTinyPivot", pyParams, solverParams);
93 amesosParams->set(solver->name(), solverParams);
94 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_SUPERLU) &&
95 Amesos2::query("superlu")) {
96 solver = Amesos2::create<Matrix, Vector>("superlu", A, X, B);
97 Teuchos::ParameterList solverParams(solver->name());
98 solverParams.set("DiagPivotThresh", sb.getDiagonalDominanceThreshold());
99 solverParams.set("ILU_DropTol", sb.getDropTolerance());
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("DiagPivotThresh", sb.getDiagonalDominanceThreshold());
124 solverParams.set("SymmetricMode", sb.isSymmetric());
125 extractParamIfSet<int>("nprocs", pyParams, solverParams);
126 extractParamIfSet<std::string>("trans", pyParams, solverParams);
127 extractParamIfSet<int>("panel_size", pyParams, solverParams);
128 extractParamIfSet<int>("relax", pyParams, solverParams);
129 extractParamIfSet<bool>("Equil", pyParams, solverParams);
130 extractParamIfSet<bool>("SymmetricMode", pyParams, solverParams);
131 extractParamIfSet<ST>("DiagPivotThresh", pyParams, solverParams);
132 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
133 amesosParams->set(solver->name(), solverParams);
134 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_PARDISO) &&
135 Amesos2::query("pardiso_mkl")) {
136 solver = Amesos2::create<Matrix, Vector>("pardiso_mkl", A, X, B);
137 Teuchos::ParameterList solverParams(solver->name());
138 extractParamIfSet<int>("IPARM(2)", pyParams, solverParams);
139 extractParamIfSet<int>("IPARM(4)", pyParams, solverParams);
140 extractParamIfSet<int>("IPARM(8)", pyParams, solverParams);
141 extractParamIfSet<int>("IPARM(10)", pyParams, solverParams);
142 extractParamIfSet<int>("IPARM(18)", pyParams, solverParams);
143 extractParamIfSet<int>("IPARM(24)", pyParams, solverParams);
144 extractParamIfSet<int>("IPARM(25)", pyParams, solverParams);
145 extractParamIfSet<int>("IPARM(60)", pyParams, solverParams);
146 amesosParams->set(solver->name(), solverParams);
147 } else if (Amesos2::query("amesos2_cholmod")) {
148 solver = Amesos2::create<Matrix, Vector>("amesos2_cholmod", A, X, B);
149 Teuchos::ParameterList solverParams(solver->name());
150 solverParams.set("DiagPivotThresh", sb.getDiagonalDominanceThreshold());
151 solverParams.set("SymmetricMode", sb.isSymmetric());
152 extractParamIfSet<std::string>("Trans", pyParams, solverParams);
153 extractParamIfSet<bool>("Equil", pyParams, solverParams);
154 extractParamIfSet<std::string>("IterRefine", pyParams, solverParams);
155 extractParamIfSet<bool>("SymmetricMode", pyParams, solverParams);
156 extractParamIfSet<ST>("DiagPivotThresh", pyParams, solverParams);
157 extractParamIfSet<std::string>("ColPerm", pyParams, solverParams);
158 amesosParams->set(solver->name(), solverParams);
159 } else if (Amesos2::query("lapack")) {
160 solver = Amesos2::create<Matrix, Vector>("lapack", A, X, B);
161 } else {
162 if (dontcare) {
163 throw TrilinosAdapterException("Could not find an Amesos2 direct solver!");
164 } else {
165 throw TrilinosAdapterException("The requested direct solver is not available!");
166 }
167 }
168 solver->setParameters(amesosParams);
169 return solver;
170 }
171
172 typedef Tpetra::CrsMatrix<real_t,LO,GO,NT> RealMatrix;
173 typedef Tpetra::CrsMatrix<cplx_t,LO,GO,NT> ComplexMatrix;
174
175 // instantiate
176 template
177 RCP<DirectSolverType<RealMatrix, RealVector> >
178 createDirectSolver<RealMatrix,RealVector>(const escript::SolverBuddy& sb,
179 RCP<const RealMatrix> A,
180 RCP<RealVector> X,
181 RCP<const RealVector> B);
182 template
183 RCP<DirectSolverType<ComplexMatrix, ComplexVector> >
184 createDirectSolver<ComplexMatrix, ComplexVector>(
185 const escript::SolverBuddy& sb,
186 RCP<const ComplexMatrix> A,
187 RCP<ComplexVector> X,
188 RCP<const ComplexVector> B);
189
190 /* Amesos2 does not currently support block matrices!
191 template
192 RCP<DirectSolverType<RealBlockMatrix, RealBlockVector> >
193 createDirectSolver<RealBlockMatrix,RealBlockVector>(
194 const escript::SolverBuddy& sb,
195 RCP<const RealBlockMatrix> A,
196 RCP<RealBlockVector> X,
197 RCP<const RealBlockVector> B);
198 template
199 RCP<DirectSolverType<ComplexBlockMatrix, ComplexBlockVector> >
200 createDirectSolver<ComplexBlockMatrix, ComplexBlockVector>(
201 const escript::SolverBuddy& sb,
202 RCP<const ComplexBlockMatrix> A,
203 RCP<ComplexBlockVector> X,
204 RCP<const ComplexBlockVector> B);
205 */
206
207 } // end of namespace
208

  ViewVC Help
Powered by ViewVC 1.1.26