/[escript]/trunk/pasowrap/src/SystemMatrixAdapter.cpp
ViewVC logotype

Contents of /trunk/pasowrap/src/SystemMatrixAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4821 - (show annotations)
Tue Apr 1 04:58:33 2014 UTC (5 years, 4 months ago) by sshaw
File size: 16786 byte(s)
moved SolverOptions to c++, split into SolverOptions for the options and SolverBuddy as the state as a precursor to per-pde solving... does break some use cases (e.g. pde.getSolverOptions().DIRECT will now fail, new value access is with SolverOptions.DIRECT), examples and documentation updated to match
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
20 using namespace std;
21
22 namespace paso {
23
24 struct null_deleter
25 {
26 void operator()(void const *ptr) const
27 {
28 }
29 };
30
31 PASOWRAP_DLL_API
32 SystemMatrixAdapter::SystemMatrixAdapter()
33 {
34 throw PasoException("Error - Illegal to generate default SystemMatrixAdapter.");
35 }
36
37 PASOWRAP_DLL_API
38 SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
39 const int row_blocksize,
40 const escript::FunctionSpace& row_functionspace,
41 const int column_blocksize,
42 const escript::FunctionSpace& column_functionspace):
43 AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
44 {
45 m_system_matrix.reset(system_matrix,null_deleter());
46 }
47
48 PASOWRAP_DLL_API
49 SystemMatrixAdapter::~SystemMatrixAdapter()
50 {
51 if (m_system_matrix.unique()) {
52 Paso_SystemMatrix* mat=m_system_matrix.get();
53 Paso_SystemMatrix_free(mat);
54 }
55 }
56
57 PASOWRAP_DLL_API
58 Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const
59 {
60 return m_system_matrix.get();
61 }
62
63 PASOWRAP_DLL_API
64 void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const
65 {
66 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
67
68 if ( x.getDataPointSize() != getColumnBlockSize()) {
69 throw PasoException("matrix vector product : column block size does not match the number of components in input.");
70 } else if (y.getDataPointSize() != getRowBlockSize()) {
71 throw PasoException("matrix vector product : row block size does not match the number of components in output.");
72 } else if ( x.getFunctionSpace() != getColumnFunctionSpace()) {
73 throw PasoException("matrix vector product : column function space and function space of input don't match.");
74 } else if (y.getFunctionSpace() != getRowFunctionSpace()) {
75 throw PasoException("matrix vector product : row function space and function space of output don't match.");
76 }
77 x.expand();
78 y.expand();
79 x.requireWrite();
80 y.requireWrite();
81 double* x_dp=x.getSampleDataRW(0);
82 double* y_dp=y.getSampleDataRW(0);
83 Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
84 checkPasoError();
85 }
86
87 PASOWRAP_DLL_API
88 int SystemMatrixAdapter::mapOptionToPaso(const int option) {
89
90 switch (option) {
91 case escript::ESCRIPT_DEFAULT:
92 return PASO_DEFAULT;
93 case escript::ESCRIPT_DIRECT:
94 return PASO_DIRECT;
95 case escript::ESCRIPT_CHOLEVSKY:
96 return PASO_CHOLEVSKY;
97 case escript::ESCRIPT_PCG:
98 return PASO_PCG;
99 case escript::ESCRIPT_CR:
100 return PASO_CR;
101 case escript::ESCRIPT_CGS:
102 return PASO_CGS;
103 case escript::ESCRIPT_BICGSTAB:
104 return PASO_BICGSTAB;
105 case escript::ESCRIPT_ILU0:
106 return PASO_ILU0;
107 case escript::ESCRIPT_ILUT:
108 return PASO_ILUT;
109 case escript::ESCRIPT_JACOBI:
110 return PASO_JACOBI;
111 case escript::ESCRIPT_GMRES:
112 return PASO_GMRES;
113 case escript::ESCRIPT_PRES20:
114 return PASO_PRES20;
115 case escript::ESCRIPT_LUMPING:
116 return PASO_LUMPING;
117 case escript::ESCRIPT_NO_REORDERING:
118 return PASO_NO_REORDERING;
119 case escript::ESCRIPT_MINIMUM_FILL_IN:
120 return PASO_MINIMUM_FILL_IN;
121 case escript::ESCRIPT_NESTED_DISSECTION:
122 return PASO_NESTED_DISSECTION;
123 case escript::ESCRIPT_MKL:
124 return PASO_MKL;
125 case escript::ESCRIPT_UMFPACK:
126 return PASO_UMFPACK;
127 case escript::ESCRIPT_ITERATIVE:
128 return PASO_ITERATIVE;
129 case escript::ESCRIPT_PASO:
130 return PASO_PASO;
131 case escript::ESCRIPT_AMG:
132 return PASO_AMG;
133 case escript::ESCRIPT_AMLI:
134 return PASO_AMLI;
135 case escript::ESCRIPT_REC_ILU:
136 return PASO_REC_ILU;
137 case escript::ESCRIPT_TRILINOS:
138 return PASO_TRILINOS;
139 case escript::ESCRIPT_NONLINEAR_GMRES:
140 return PASO_NONLINEAR_GMRES;
141 case escript::ESCRIPT_TFQMR :
142 return PASO_TFQMR;
143 case escript::ESCRIPT_MINRES:
144 return PASO_MINRES;
145 case escript::ESCRIPT_GAUSS_SEIDEL:
146 return PASO_GAUSS_SEIDEL;
147 case escript::ESCRIPT_RILU:
148 return PASO_RILU;
149 case escript::ESCRIPT_DEFAULT_REORDERING:
150 return PASO_DEFAULT_REORDERING;
151 case escript::ESCRIPT_SUPER_LU:
152 return PASO_SUPER_LU;
153 case escript::ESCRIPT_PASTIX:
154 return PASO_PASTIX;
155 case escript::ESCRIPT_YAIR_SHAPIRA_COARSENING:
156 return PASO_YAIR_SHAPIRA_COARSENING;
157 case escript::ESCRIPT_RUGE_STUEBEN_COARSENING:
158 return PASO_RUGE_STUEBEN_COARSENING;
159 case escript::ESCRIPT_STANDARD_COARSENING:
160 return PASO_STANDARD_COARSENING;
161 case escript::ESCRIPT_AGGREGATION_COARSENING:
162 return PASO_AGGREGATION_COARSENING;
163 case escript::ESCRIPT_NO_PRECONDITIONER:
164 return PASO_NO_PRECONDITIONER;
165 case escript::ESCRIPT_CLASSIC_INTERPOLATION_WITH_FF_COUPLING:
166 return PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING;
167 case escript::ESCRIPT_CLASSIC_INTERPOLATION:
168 return PASO_CLASSIC_INTERPOLATION;
169 case escript::ESCRIPT_DIRECT_INTERPOLATION:
170 return PASO_DIRECT_INTERPOLATION;
171 case escript::ESCRIPT_BOOMERAMG:
172 return PASO_BOOMERAMG;
173 case escript::ESCRIPT_CIJP_FIXED_RANDOM_COARSENING:
174 return PASO_CIJP_FIXED_RANDOM_COARSENING;
175 case escript::ESCRIPT_CIJP_COARSENING:
176 return PASO_CIJP_COARSENING;
177 case escript::ESCRIPT_FALGOUT_COARSENING:
178 return PASO_FALGOUT_COARSENING;
179 case escript::ESCRIPT_PMIS_COARSENING:
180 return PASO_PMIS_COARSENING;
181 case escript::ESCRIPT_HMIS_COARSENING:
182 return PASO_HMIS_COARSENING;
183 case escript::ESCRIPT_LINEAR_CRANK_NICOLSON:
184 return PASO_LINEAR_CRANK_NICOLSON;
185 case escript::ESCRIPT_CRANK_NICOLSON:
186 return PASO_CRANK_NICOLSON;
187 case escript::ESCRIPT_BACKWARD_EULER:
188 return PASO_BACKWARD_EULER;
189 default:
190 stringstream temp;
191 temp << "Error - Cannot map option value "<< option << " onto Paso";
192 throw PasoException(temp.str());
193 }
194 }
195
196 PASOWRAP_DLL_API
197 int SystemMatrixAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner,
198 const int package, const bool symmetry, Esys_MPIInfo* mpiInfo)
199 {
200 int out=Paso_SystemMatrix_getSystemMatrixTypeId(mapOptionToPaso(solver),
201 mapOptionToPaso(preconditioner), mapOptionToPaso(package),
202 symmetry?1:0, mpiInfo);
203 checkPasoError();
204 return out;
205 }
206
207 PASOWRAP_DLL_API
208 void SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
209 {
210 Paso_SystemMatrix* mat=m_system_matrix.get();
211 int first_row_index = mat->row_distribution->first_component[mat->mpi_info->rank];
212 int last_row_index = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
213 int first_col_index = mat->col_distribution->first_component[mat->mpi_info->rank];
214 int last_col_index = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;
215
216 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
217
218 switch (mat->type) {
219 case MATRIX_FORMAT_DEFAULT:
220 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n");
221 break;
222 case MATRIX_FORMAT_CSC:
223 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n");
224 break;
225 case MATRIX_FORMAT_BLK1:
226 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n");
227 break;
228 case MATRIX_FORMAT_OFFSET1:
229 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n");
230 break;
231 case MATRIX_FORMAT_TRILINOS_CRS:
232 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n");
233 break;
234 default:
235 fprintf(stdout, "\tMatrix type unknown\n");
236 break;
237 }
238
239 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
240 fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
241 fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
242 fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
243 fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
244 fprintf(stdout, "\tcol_coupleBlock numRows %d\n", mat->col_coupleBlock->numRows);
245 fprintf(stdout, "\tcol_coupleBlock numCols %d\n", mat->col_coupleBlock->numCols);
246 fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", mat->col_coupleBlock->pattern->numOutput);
247 fprintf(stdout, "\trow_coupleBlock numRows %d\n", mat->row_coupleBlock->numRows);
248 fprintf(stdout, "\trow_coupleBlock numCols %d\n", mat->row_coupleBlock->numCols);
249 fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", mat->row_coupleBlock->pattern->numOutput);
250 fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
251 fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
252 fprintf(stdout, "\tblock_size %d\n", mat->block_size);
253 fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
254 fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
255
256 }
257
258 PASOWRAP_DLL_API
259 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, boost::python::object& options) const
260 {
261 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
262 Paso_Options paso_options;
263 options.attr("resetDiagnostics")();
264 escriptToPasoOptions(&paso_options,options);
265 if ( out.getDataPointSize() != getColumnBlockSize()) {
266 throw PasoException("solve : column block size does not match the number of components of solution.");
267 } else if ( in.getDataPointSize() != getRowBlockSize()) {
268 throw PasoException("solve : row block size does not match the number of components of right hand side.");
269 } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
270 throw PasoException("solve : column function space and function space of solution don't match.");
271 } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
272 throw PasoException("solve : row function space and function space of right hand side don't match.");
273 }
274 out.expand();
275 in.expand();
276 double* out_dp=out.getSampleDataRW(0);
277 double* in_dp=in.getSampleDataRW(0);
278 Paso_solve(mat,out_dp,in_dp,&paso_options);
279 pasoToEscriptOptions(&paso_options,options);
280 checkPasoError();
281
282 }
283
284 PASOWRAP_DLL_API
285 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
286 {
287 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
288 if ( col_q.getDataPointSize() != getColumnBlockSize()) {
289 throw PasoException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
290 } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
291 throw PasoException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
292 } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
293 throw PasoException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
294 } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
295 throw PasoException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
296 }
297 row_q.expand();
298 col_q.expand();
299 row_q.requireWrite();
300 col_q.requireWrite();
301 double* row_q_dp=row_q.getSampleDataRW(0);
302 double* col_q_dp=col_q.getSampleDataRW(0);
303 Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
304 checkPasoError();
305 }
306
307 PASOWRAP_DLL_API
308 void SystemMatrixAdapter::saveMM(const std::string& fileName) const
309 {
310 if( fileName.size() == 0 )
311 {
312 throw PasoException("Null file name!");
313 }
314
315 char *fName = TMPMEMALLOC(fileName.size()+1,char);
316
317 strcpy(fName,fileName.c_str());
318 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
319 Paso_SystemMatrix_saveMM(mat,fName);
320 checkPasoError();
321 TMPMEMFREE(fName);
322
323 }
324
325 PASOWRAP_DLL_API
326 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
327 {
328 if( fileName.size() == 0 )
329 {
330 throw PasoException("Null file name!");
331 }
332
333 char *fName = TMPMEMALLOC(fileName.size()+1,char);
334
335 strcpy(fName,fileName.c_str());
336 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
337 Paso_SystemMatrix_saveHB(mat,fName);
338 checkPasoError();
339 TMPMEMFREE(fName);
340
341 }
342
343 PASOWRAP_DLL_API
344 void SystemMatrixAdapter::resetValues() const
345 {
346 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
347 Paso_SystemMatrix_setValues(mat,0.);
348 Paso_solve_free(mat);
349 checkPasoError();
350 }
351
352 PASOWRAP_DLL_API
353 void SystemMatrixAdapter::pasoToEscriptOptions(const Paso_Options* paso_options,boost::python::object& options)
354 {
355 #define SET(__key__,__val__,__type__) options.attr("_updateDiagnostics")(__key__,(__type__)paso_options->__val__)
356 SET("num_iter", num_iter, int);
357 SET("num_level", num_level, int);
358 SET("num_inner_iter", num_inner_iter, int);
359 SET("time", time, double);
360 SET("set_up_time", set_up_time, double);
361 SET("net_time", net_time, double);
362 SET("residual_norm", residual_norm, double);
363 SET("converged",converged, bool);
364 SET("time_step_backtracking_used", time_step_backtracking_used,bool);
365 SET("coarse_level_sparsity",coarse_level_sparsity,double);
366 SET("num_coarse_unknowns",num_coarse_unknowns,int);
367 #undef SET
368 }
369
370 PASOWRAP_DLL_API
371 void SystemMatrixAdapter::escriptToPasoOptions(Paso_Options* paso_options, const boost::python::object& options)
372 {
373
374 escript::SolverBuddy sb = boost::python::extract<escript::SolverBuddy>(options);
375
376 Paso_Options_setDefaults(paso_options);
377 paso_options->method = mapOptionToPaso(sb.getSolverMethod());
378 paso_options->package = mapOptionToPaso(sb.getPackage());
379 paso_options->verbose = sb.isVerbose();
380 paso_options->symmetric = sb.isSymmetric();
381 paso_options->tolerance = sb.getTolerance();
382 paso_options->absolute_tolerance = sb.getAbsoluteTolerance();
383 paso_options->inner_tolerance = sb.getInnerTolerance();
384 paso_options->adapt_inner_tolerance = sb.adaptInnerTolerance();
385 paso_options->reordering = mapOptionToPaso(sb.getReordering());
386 paso_options->preconditioner = mapOptionToPaso(sb.getPreconditioner());
387 paso_options->iter_max = sb.getIterMax();
388 paso_options->inner_iter_max = sb.getInnerIterMax();
389 paso_options->drop_tolerance = sb.getDropTolerance();
390 paso_options->drop_storage = sb.getDropStorage();
391 paso_options->truncation = sb.getTruncation();
392 paso_options->restart = sb._getRestartForC();
393 paso_options->sweeps = sb.getNumSweeps();
394 paso_options->pre_sweeps = sb.getNumPreSweeps();
395 paso_options->post_sweeps = sb.getNumPostSweeps();
396 paso_options->level_max = sb.getLevelMax();
397 paso_options->min_coarse_matrix_size = sb.getMinCoarseMatrixSize();
398 paso_options->coarsening_threshold = sb.getCoarseningThreshold();
399 paso_options->accept_failed_convergence = sb.acceptConvergenceFailure();
400 paso_options->coarsening_method = mapOptionToPaso(sb.getCoarsening());
401 paso_options->smoother = mapOptionToPaso(sb.getSmoother());
402 paso_options->relaxation_factor = sb.getRelaxationFactor();
403 paso_options->use_local_preconditioner = sb.useLocalPreconditioner();
404 paso_options->min_coarse_sparsity = sb.getMinCoarseMatrixSparsity();
405 paso_options->refinements = sb.getNumRefinements();
406 paso_options->coarse_matrix_refinements = sb.getNumCoarseMatrixRefinements();
407 paso_options->usePanel = sb.usePanel();
408 paso_options->interpolation_method = sb.getAMGInterpolation();
409 paso_options->diagonal_dominance_threshold = sb.getDiagonalDominanceThreshold();
410 paso_options->ode_solver = sb.getODESolver();
411 }
412
413
414 } // end of namespace

  ViewVC Help
Powered by ViewVC 1.1.26