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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6799 - (hide annotations)
Mon Mar 25 05:53:58 2019 UTC (8 weeks, 1 day 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 caltinay 6104
2     /*****************************************************************************
3     *
4     * Copyright (c) 2016 by The University of Queensland
5     * http://www.uq.edu.au
6     *
7     * Primary Business: Queensland, Australia
8 caltinay 6116 * Licensed under the Apache License, version 2.0
9     * http://www.apache.org/licenses/LICENSE-2.0
10 caltinay 6104 *
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 caltinay 6375 #include <trilinoswrap/util.h>
20 caltinay 6104
21     #include <escript/SolverOptions.h>
22    
23     #include <Amesos2.hpp>
24 caltinay 6181 #include <Tpetra_CrsMatrix.hpp>
25 caltinay 6104
26 caltinay 6375 #include <boost/python/dict.hpp>
27    
28 caltinay 6104 using Teuchos::RCP;
29    
30 caltinay 6375 namespace bp = boost::python;
31    
32 caltinay 6104 namespace esys_trilinos {
33    
34 caltinay 6181 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 caltinay 6104 {
41 caltinay 6375 using util::extractParamIfSet;
42    
43     typedef typename Matrix::scalar_type ST;
44    
45 caltinay 6181 RCP<DirectSolverType<Matrix,Vector> > solver;
46 caltinay 6383 RCP<Teuchos::ParameterList> amesosParams = Teuchos::parameterList("Amesos2");
47 caltinay 6375 const bp::dict& pyParams = sb.getTrilinosParameters();
48 caltinay 6104
49 caltinay 6392 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 caltinay 6396 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 caltinay 6392 Amesos2::query("MUMPS")) {
67 caltinay 6181 solver = Amesos2::create<Matrix, Vector>("MUMPS", A, X, B);
68 caltinay 6383 Teuchos::ParameterList solverParams(solver->name());
69 aellery 6799 // solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
70     if (sb.isVerbose()) {
71 caltinay 6383 solverParams.set("ICNTL(4)", 4);
72 caltinay 6104 }
73 caltinay 6383 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 caltinay 6392 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_TRILINOS) &&
82     Amesos2::query("Basker")) {
83     solver = Amesos2::create<Matrix, Vector>("Basker", A, X, B);
84 aellery 6799 Teuchos::ParameterList solverParams(solver->name());
85     solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
86 caltinay 6392 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_SUPERLU) &&
87     Amesos2::query("superludist")) {
88 caltinay 6181 solver = Amesos2::create<Matrix, Vector>("superludist", A, X, B);
89 aellery 6799 Teuchos::ParameterList solverParams(solver->name());
90     solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
91 caltinay 6383 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 caltinay 6392 } else if ((dontcare || method == escript::SO_METHOD_DIRECT_SUPERLU) &&
97     Amesos2::query("superlu")) {
98 caltinay 6181 solver = Amesos2::create<Matrix, Vector>("superlu", A, X, B);
99 caltinay 6383 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 caltinay 6392 } 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 jfenwick 6574 #ifdef _OPENMP
119 caltinay 6392 solverParams.set("nprocs", omp_get_max_threads());
120 jfenwick 6574 #else
121     solverParams.set("nprocs", 1);
122     #endif
123 caltinay 6392 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 caltinay 6181 solver = Amesos2::create<Matrix, Vector>("pardiso_mkl", A, X, B);
136 caltinay 6383 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 caltinay 6104 } else if (Amesos2::query("amesos2_cholmod")) {
147 caltinay 6181 solver = Amesos2::create<Matrix, Vector>("amesos2_cholmod", A, X, B);
148 caltinay 6383 Teuchos::ParameterList solverParams(solver->name());
149 aellery 6799
150 caltinay 6383 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 caltinay 6392 } else if (Amesos2::query("lapack")) {
159     solver = Amesos2::create<Matrix, Vector>("lapack", A, X, B);
160 aellery 6799 Teuchos::ParameterList solverParams(solver->name());
161     solverParams.set("MatrixType", (sb.isSymmetric() || sb.isHermitian()) ? "symmetric" : "general");
162 caltinay 6104 } else {
163 caltinay 6392 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 caltinay 6104 }
169 caltinay 6383 solver->setParameters(amesosParams);
170 caltinay 6104 return solver;
171     }
172    
173 caltinay 6181 typedef Tpetra::CrsMatrix<real_t,LO,GO,NT> RealMatrix;
174     typedef Tpetra::CrsMatrix<cplx_t,LO,GO,NT> ComplexMatrix;
175    
176     // instantiate
177 caltinay 6104 template
178 caltinay 6181 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 caltinay 6104 template
184 caltinay 6181 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 caltinay 6104
191 caltinay 6181 /* 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 caltinay 6104 } // end of namespace
209    

  ViewVC Help
Powered by ViewVC 1.1.26