/[escript]/branches/diaplayground/pasowrap/src/SystemMatrixAdapter.cpp
ViewVC logotype

Contents of /branches/diaplayground/pasowrap/src/SystemMatrixAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5136 - (show annotations)
Tue Sep 9 07:13:55 2014 UTC (4 years, 5 months ago) by caltinay
File size: 15953 byte(s)
ripley now supports paso solvers again and returns an appropriate matrix type
id. Changed the getSystemMatrixTypeId() method to take a full SolverBuddy
instance and made some other simplifications.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by 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 "SystemMatrixAdapter.h"
18 #include <escript/SolverOptions.h>
19 #include <paso/Solver.h>
20
21 using namespace std;
22
23 namespace paso {
24
25 SystemMatrixAdapter::SystemMatrixAdapter()
26 {
27 throw PasoException("Error - Illegal to generate default SystemMatrixAdapter.");
28 }
29
30 SystemMatrixAdapter::SystemMatrixAdapter(SystemMatrix_ptr system_matrix,
31 int row_blocksize, const escript::FunctionSpace& row_fs,
32 int column_blocksize, const escript::FunctionSpace& column_fs) :
33 AbstractSystemMatrix(row_blocksize,row_fs,column_blocksize,column_fs),
34 m_system_matrix(system_matrix)
35 {
36 }
37
38 SystemMatrixAdapter::~SystemMatrixAdapter()
39 {
40 }
41
42 SystemMatrix_ptr SystemMatrixAdapter::getPaso_SystemMatrix() const
43 {
44 return m_system_matrix;
45 }
46
47 void SystemMatrixAdapter::ypAx(escript::Data& y, escript::Data& x) const
48 {
49 if ( x.getDataPointSize() != getColumnBlockSize()) {
50 throw PasoException("matrix vector product : column block size does not match the number of components in input.");
51 } else if (y.getDataPointSize() != getRowBlockSize()) {
52 throw PasoException("matrix vector product : row block size does not match the number of components in output.");
53 } else if ( x.getFunctionSpace() != getColumnFunctionSpace()) {
54 throw PasoException("matrix vector product : column function space and function space of input don't match.");
55 } else if (y.getFunctionSpace() != getRowFunctionSpace()) {
56 throw PasoException("matrix vector product : row function space and function space of output don't match.");
57 }
58 x.expand();
59 y.expand();
60 x.requireWrite();
61 y.requireWrite();
62 double* x_dp=x.getSampleDataRW(0);
63 double* y_dp=y.getSampleDataRW(0);
64 SystemMatrix_MatrixVector(1., m_system_matrix, x_dp, 1., y_dp);
65 checkPasoError();
66 }
67
68 int SystemMatrixAdapter::mapOptionToPaso(int option)
69 {
70 switch (option) {
71 case escript::SO_DEFAULT:
72 return PASO_DEFAULT;
73
74 case escript::SO_PACKAGE_MKL:
75 return PASO_MKL;
76 case escript::SO_PACKAGE_PASO:
77 return PASO_PASO;
78 case escript::SO_PACKAGE_PASTIX:
79 return PASO_PASTIX;
80 case escript::SO_PACKAGE_SUPER_LU:
81 return PASO_SUPER_LU;
82 case escript::SO_PACKAGE_TRILINOS:
83 return PASO_TRILINOS;
84 case escript::SO_PACKAGE_UMFPACK:
85 return PASO_UMFPACK;
86
87 case escript::SO_METHOD_BICGSTAB:
88 return PASO_BICGSTAB;
89 case escript::SO_METHOD_CGS:
90 return PASO_CGS;
91 case escript::SO_METHOD_CHOLEVSKY:
92 return PASO_CHOLEVSKY;
93 case escript::SO_METHOD_CR:
94 return PASO_CR;
95 case escript::SO_METHOD_DIRECT:
96 return PASO_DIRECT;
97 case escript::SO_METHOD_GMRES:
98 return PASO_GMRES;
99 case escript::SO_METHOD_ITERATIVE:
100 return PASO_ITERATIVE;
101 case escript::SO_METHOD_MINRES:
102 return PASO_MINRES;
103 case escript::SO_METHOD_NONLINEAR_GMRES:
104 return PASO_NONLINEAR_GMRES;
105 case escript::SO_METHOD_PCG:
106 return PASO_PCG;
107 case escript::SO_METHOD_PRES20:
108 return PASO_PRES20;
109 case escript::SO_METHOD_TFQMR:
110 return PASO_TFQMR;
111
112 case escript::SO_PRECONDITIONER_AMG:
113 return PASO_AMG;
114 case escript::SO_PRECONDITIONER_AMLI:
115 return PASO_AMLI;
116 case escript::SO_PRECONDITIONER_BOOMERAMG:
117 return PASO_BOOMERAMG;
118 case escript::SO_PRECONDITIONER_GAUSS_SEIDEL:
119 return PASO_GAUSS_SEIDEL;
120 case escript::SO_PRECONDITIONER_ILU0:
121 return PASO_ILU0;
122 case escript::SO_PRECONDITIONER_ILUT:
123 return PASO_ILUT;
124 case escript::SO_PRECONDITIONER_JACOBI:
125 return PASO_JACOBI;
126 case escript::SO_PRECONDITIONER_NONE:
127 return PASO_NO_PRECONDITIONER;
128 case escript::SO_PRECONDITIONER_REC_ILU:
129 return PASO_REC_ILU;
130 case escript::SO_PRECONDITIONER_RILU:
131 return PASO_RILU;
132
133 case escript::SO_ODESOLVER_BACKWARD_EULER:
134 return PASO_BACKWARD_EULER;
135 case escript::SO_ODESOLVER_CRANK_NICOLSON:
136 return PASO_CRANK_NICOLSON;
137 case escript::SO_ODESOLVER_LINEAR_CRANK_NICOLSON:
138 return PASO_LINEAR_CRANK_NICOLSON;
139
140 case escript::SO_INTERPOLATION_CLASSIC:
141 return PASO_CLASSIC_INTERPOLATION;
142 case escript::SO_INTERPOLATION_CLASSIC_WITH_FF_COUPLING:
143 return PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING;
144 case escript::SO_INTERPOLATION_DIRECT:
145 return PASO_DIRECT_INTERPOLATION;
146
147 case escript::SO_COARSENING_AGGREGATION:
148 return PASO_AGGREGATION_COARSENING;
149 case escript::SO_COARSENING_CIJP:
150 return PASO_CIJP_COARSENING;
151 case escript::SO_COARSENING_CIJP_FIXED_RANDOM:
152 return PASO_CIJP_FIXED_RANDOM_COARSENING;
153 case escript::SO_COARSENING_FALGOUT:
154 return PASO_FALGOUT_COARSENING;
155 case escript::SO_COARSENING_HMIS:
156 return PASO_HMIS_COARSENING;
157 case escript::SO_COARSENING_PMIS:
158 return PASO_PMIS_COARSENING;
159 case escript::SO_COARSENING_RUGE_STUEBEN:
160 return PASO_RUGE_STUEBEN_COARSENING;
161 case escript::SO_COARSENING_STANDARD:
162 return PASO_STANDARD_COARSENING;
163 case escript::SO_COARSENING_YAIR_SHAPIRA:
164 return PASO_YAIR_SHAPIRA_COARSENING;
165
166 case escript::SO_REORDERING_DEFAULT:
167 return PASO_DEFAULT_REORDERING;
168 case escript::SO_REORDERING_MINIMUM_FILL_IN:
169 return PASO_MINIMUM_FILL_IN;
170 case escript::SO_REORDERING_NESTED_DISSECTION:
171 return PASO_NESTED_DISSECTION;
172 case escript::SO_REORDERING_NONE:
173 return PASO_NO_REORDERING;
174
175 default:
176 stringstream temp;
177 temp << "Error - Cannot map option value "<< option << " onto Paso";
178 throw PasoException(temp.str());
179 }
180 }
181
182 int SystemMatrixAdapter::getSystemMatrixTypeId(int solver, int preconditioner,
183 int package, bool symmetry, const esysUtils::JMPI& mpiInfo)
184 {
185 int out=SystemMatrix::getSystemMatrixTypeId(mapOptionToPaso(solver),
186 mapOptionToPaso(preconditioner), mapOptionToPaso(package),
187 symmetry?1:0, mpiInfo);
188 checkPasoError();
189 return out;
190 }
191
192 void SystemMatrixAdapter::Print_Matrix_Info(bool full=false) const
193 {
194 int first_row_index = m_system_matrix->row_distribution->first_component[m_system_matrix->mpi_info->rank];
195 int last_row_index = m_system_matrix->row_distribution->first_component[m_system_matrix->mpi_info->rank+1]-1;
196 int first_col_index = m_system_matrix->col_distribution->first_component[m_system_matrix->mpi_info->rank];
197 int last_col_index = m_system_matrix->col_distribution->first_component[m_system_matrix->mpi_info->rank+1]-1;
198
199 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n",
200 m_system_matrix->mpi_info->rank, m_system_matrix->mpi_info->size);
201
202 switch (m_system_matrix->type) {
203 case MATRIX_FORMAT_DEFAULT:
204 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n");
205 break;
206 case MATRIX_FORMAT_CSC:
207 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n");
208 break;
209 case MATRIX_FORMAT_BLK1:
210 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n");
211 break;
212 case MATRIX_FORMAT_OFFSET1:
213 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n");
214 break;
215 case MATRIX_FORMAT_TRILINOS_CRS:
216 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n");
217 break;
218 default:
219 fprintf(stdout, "\tMatrix type unknown\n");
220 break;
221 }
222
223 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
224 fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
225 fprintf(stdout, "\tmainBlock numRows %d\n", m_system_matrix->mainBlock->numRows);
226 fprintf(stdout, "\tmainBlock numCols %d\n", m_system_matrix->mainBlock->numCols);
227 fprintf(stdout, "\tmainBlock pattern numOutput %d\n", m_system_matrix->mainBlock->pattern->numOutput);
228 fprintf(stdout, "\tcol_coupleBlock numRows %d\n", m_system_matrix->col_coupleBlock->numRows);
229 fprintf(stdout, "\tcol_coupleBlock numCols %d\n", m_system_matrix->col_coupleBlock->numCols);
230 fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", m_system_matrix->col_coupleBlock->pattern->numOutput);
231 fprintf(stdout, "\trow_coupleBlock numRows %d\n", m_system_matrix->row_coupleBlock->numRows);
232 fprintf(stdout, "\trow_coupleBlock numCols %d\n", m_system_matrix->row_coupleBlock->numCols);
233 fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", m_system_matrix->row_coupleBlock->pattern->numOutput);
234 fprintf(stdout, "\trow_block_size %d\n", m_system_matrix->row_block_size);
235 fprintf(stdout, "\tcol_block_size %d\n", m_system_matrix->col_block_size);
236 fprintf(stdout, "\tblock_size %d\n", m_system_matrix->block_size);
237 fprintf(stdout, "\tlogical_row_block_size %d\n", m_system_matrix->logical_row_block_size);
238 fprintf(stdout, "\tlogical_col_block_size %d\n", m_system_matrix->logical_col_block_size);
239
240 }
241
242 void SystemMatrixAdapter::setToSolution(escript::Data& out, escript::Data& in,
243 boost::python::object& options) const
244 {
245 Options paso_options;
246 options.attr("resetDiagnostics")();
247 escriptToPasoOptions(&paso_options,options);
248 if ( out.getDataPointSize() != getColumnBlockSize()) {
249 throw PasoException("solve : column block size does not match the number of components of solution.");
250 } else if ( in.getDataPointSize() != getRowBlockSize()) {
251 throw PasoException("solve : row block size does not match the number of components of right hand side.");
252 } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
253 throw PasoException("solve : column function space and function space of solution don't match.");
254 } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
255 throw PasoException("solve : row function space and function space of right hand side don't match.");
256 }
257 out.expand();
258 in.expand();
259 double* out_dp=out.getSampleDataRW(0);
260 double* in_dp=in.getSampleDataRW(0);
261 paso::solve(m_system_matrix, out_dp, in_dp, &paso_options);
262 pasoToEscriptOptions(&paso_options,options);
263 checkPasoError();
264 }
265
266 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,
267 escript::Data& col_q, double mdv)
268 {
269 if ( col_q.getDataPointSize() != getColumnBlockSize()) {
270 throw PasoException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
271 } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
272 throw PasoException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
273 } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
274 throw PasoException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
275 } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
276 throw PasoException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
277 }
278 row_q.expand();
279 col_q.expand();
280 row_q.requireWrite();
281 col_q.requireWrite();
282 double* row_q_dp=row_q.getSampleDataRW(0);
283 double* col_q_dp=col_q.getSampleDataRW(0);
284 m_system_matrix->nullifyRowsAndCols(row_q_dp, col_q_dp, mdv);
285 checkPasoError();
286 }
287
288 void SystemMatrixAdapter::saveMM(const std::string& filename) const
289 {
290 m_system_matrix->saveMM(filename.c_str());
291 checkPasoError();
292 }
293
294 void SystemMatrixAdapter::saveHB(const std::string& filename) const
295 {
296 m_system_matrix->saveHB(filename.c_str());
297 checkPasoError();
298 }
299
300 void SystemMatrixAdapter::resetValues()
301 {
302 m_system_matrix->setValues(0.);
303 solve_free(m_system_matrix.get());
304 checkPasoError();
305 }
306
307 void SystemMatrixAdapter::pasoToEscriptOptions(const Options* paso_options,
308 boost::python::object& options)
309 {
310 #define SET(__key__,__val__,__type__) options.attr("_updateDiagnostics")(__key__,(__type__)paso_options->__val__)
311 SET("num_iter", num_iter, int);
312 SET("num_level", num_level, int);
313 SET("num_inner_iter", num_inner_iter, int);
314 SET("time", time, double);
315 SET("set_up_time", set_up_time, double);
316 SET("net_time", net_time, double);
317 SET("residual_norm", residual_norm, double);
318 SET("converged",converged, bool);
319 SET("time_step_backtracking_used", time_step_backtracking_used,bool);
320 SET("coarse_level_sparsity",coarse_level_sparsity,double);
321 SET("num_coarse_unknowns",num_coarse_unknowns,int);
322 #undef SET
323 }
324
325 void SystemMatrixAdapter::escriptToPasoOptions(Options* paso_options,
326 const boost::python::object& options)
327 {
328 escript::SolverBuddy sb = boost::python::extract<escript::SolverBuddy>(options);
329
330 paso_options->setDefaults();
331 paso_options->method = mapOptionToPaso(sb.getSolverMethod());
332 paso_options->package = mapOptionToPaso(sb.getPackage());
333 paso_options->verbose = sb.isVerbose();
334 paso_options->symmetric = sb.isSymmetric();
335 paso_options->tolerance = sb.getTolerance();
336 paso_options->absolute_tolerance = sb.getAbsoluteTolerance();
337 paso_options->inner_tolerance = sb.getInnerTolerance();
338 paso_options->adapt_inner_tolerance = sb.adaptInnerTolerance();
339 paso_options->reordering = mapOptionToPaso(sb.getReordering());
340 paso_options->preconditioner = mapOptionToPaso(sb.getPreconditioner());
341 paso_options->ode_solver = mapOptionToPaso(sb.getODESolver());
342 paso_options->iter_max = sb.getIterMax();
343 paso_options->inner_iter_max = sb.getInnerIterMax();
344 paso_options->drop_tolerance = sb.getDropTolerance();
345 paso_options->drop_storage = sb.getDropStorage();
346 paso_options->truncation = sb.getTruncation();
347 paso_options->restart = sb._getRestartForC();
348 paso_options->sweeps = sb.getNumSweeps();
349 paso_options->pre_sweeps = sb.getNumPreSweeps();
350 paso_options->post_sweeps = sb.getNumPostSweeps();
351 paso_options->level_max = sb.getLevelMax();
352 paso_options->min_coarse_matrix_size = sb.getMinCoarseMatrixSize();
353 paso_options->coarsening_threshold = sb.getCoarseningThreshold();
354 paso_options->accept_failed_convergence = sb.acceptConvergenceFailure();
355 paso_options->coarsening_method = mapOptionToPaso(sb.getCoarsening());
356 paso_options->smoother = mapOptionToPaso(sb.getSmoother());
357 paso_options->relaxation_factor = sb.getRelaxationFactor();
358 paso_options->use_local_preconditioner = sb.useLocalPreconditioner();
359 paso_options->min_coarse_sparsity = sb.getMinCoarseMatrixSparsity();
360 paso_options->refinements = sb.getNumRefinements();
361 paso_options->coarse_matrix_refinements = sb.getNumCoarseMatrixRefinements();
362 paso_options->usePanel = sb.usePanel();
363 paso_options->interpolation_method = sb.getAMGInterpolation();
364 paso_options->diagonal_dominance_threshold = sb.getDiagonalDominanceThreshold();
365 }
366
367
368 } // end of namespace
369

  ViewVC Help
Powered by ViewVC 1.1.26