/[escript]/release/3.4.2/pasowrap/src/SystemMatrixAdapter.cpp
ViewVC logotype

Contents of /release/3.4.2/pasowrap/src/SystemMatrixAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4919 - (show annotations)
Wed Apr 30 06:25:55 2014 UTC (6 years, 1 month ago) by jfenwick
File size: 15765 byte(s)
Because we've got to!

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 const int row_blocksize, const escript::FunctionSpace& row_fs,
32 const 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(const int option) {
69
70 switch (option) {
71 case escript::ESCRIPT_DEFAULT:
72 return PASO_DEFAULT;
73 case escript::ESCRIPT_DIRECT:
74 return PASO_DIRECT;
75 case escript::ESCRIPT_CHOLEVSKY:
76 return PASO_CHOLEVSKY;
77 case escript::ESCRIPT_PCG:
78 return PASO_PCG;
79 case escript::ESCRIPT_CR:
80 return PASO_CR;
81 case escript::ESCRIPT_CGS:
82 return PASO_CGS;
83 case escript::ESCRIPT_BICGSTAB:
84 return PASO_BICGSTAB;
85 case escript::ESCRIPT_ILU0:
86 return PASO_ILU0;
87 case escript::ESCRIPT_ILUT:
88 return PASO_ILUT;
89 case escript::ESCRIPT_JACOBI:
90 return PASO_JACOBI;
91 case escript::ESCRIPT_GMRES:
92 return PASO_GMRES;
93 case escript::ESCRIPT_PRES20:
94 return PASO_PRES20;
95 case escript::ESCRIPT_LUMPING:
96 return PASO_LUMPING;
97 case escript::ESCRIPT_NO_REORDERING:
98 return PASO_NO_REORDERING;
99 case escript::ESCRIPT_MINIMUM_FILL_IN:
100 return PASO_MINIMUM_FILL_IN;
101 case escript::ESCRIPT_NESTED_DISSECTION:
102 return PASO_NESTED_DISSECTION;
103 case escript::ESCRIPT_MKL:
104 return PASO_MKL;
105 case escript::ESCRIPT_UMFPACK:
106 return PASO_UMFPACK;
107 case escript::ESCRIPT_ITERATIVE:
108 return PASO_ITERATIVE;
109 case escript::ESCRIPT_PASO:
110 return PASO_PASO;
111 case escript::ESCRIPT_AMG:
112 return PASO_AMG;
113 case escript::ESCRIPT_AMLI:
114 return PASO_AMLI;
115 case escript::ESCRIPT_REC_ILU:
116 return PASO_REC_ILU;
117 case escript::ESCRIPT_TRILINOS:
118 return PASO_TRILINOS;
119 case escript::ESCRIPT_NONLINEAR_GMRES:
120 return PASO_NONLINEAR_GMRES;
121 case escript::ESCRIPT_TFQMR :
122 return PASO_TFQMR;
123 case escript::ESCRIPT_MINRES:
124 return PASO_MINRES;
125 case escript::ESCRIPT_GAUSS_SEIDEL:
126 return PASO_GAUSS_SEIDEL;
127 case escript::ESCRIPT_RILU:
128 return PASO_RILU;
129 case escript::ESCRIPT_DEFAULT_REORDERING:
130 return PASO_DEFAULT_REORDERING;
131 case escript::ESCRIPT_SUPER_LU:
132 return PASO_SUPER_LU;
133 case escript::ESCRIPT_PASTIX:
134 return PASO_PASTIX;
135 case escript::ESCRIPT_YAIR_SHAPIRA_COARSENING:
136 return PASO_YAIR_SHAPIRA_COARSENING;
137 case escript::ESCRIPT_RUGE_STUEBEN_COARSENING:
138 return PASO_RUGE_STUEBEN_COARSENING;
139 case escript::ESCRIPT_STANDARD_COARSENING:
140 return PASO_STANDARD_COARSENING;
141 case escript::ESCRIPT_AGGREGATION_COARSENING:
142 return PASO_AGGREGATION_COARSENING;
143 case escript::ESCRIPT_NO_PRECONDITIONER:
144 return PASO_NO_PRECONDITIONER;
145 case escript::ESCRIPT_CLASSIC_INTERPOLATION_WITH_FF_COUPLING:
146 return PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING;
147 case escript::ESCRIPT_CLASSIC_INTERPOLATION:
148 return PASO_CLASSIC_INTERPOLATION;
149 case escript::ESCRIPT_DIRECT_INTERPOLATION:
150 return PASO_DIRECT_INTERPOLATION;
151 case escript::ESCRIPT_BOOMERAMG:
152 return PASO_BOOMERAMG;
153 case escript::ESCRIPT_CIJP_FIXED_RANDOM_COARSENING:
154 return PASO_CIJP_FIXED_RANDOM_COARSENING;
155 case escript::ESCRIPT_CIJP_COARSENING:
156 return PASO_CIJP_COARSENING;
157 case escript::ESCRIPT_FALGOUT_COARSENING:
158 return PASO_FALGOUT_COARSENING;
159 case escript::ESCRIPT_PMIS_COARSENING:
160 return PASO_PMIS_COARSENING;
161 case escript::ESCRIPT_HMIS_COARSENING:
162 return PASO_HMIS_COARSENING;
163 case escript::ESCRIPT_LINEAR_CRANK_NICOLSON:
164 return PASO_LINEAR_CRANK_NICOLSON;
165 case escript::ESCRIPT_CRANK_NICOLSON:
166 return PASO_CRANK_NICOLSON;
167 case escript::ESCRIPT_BACKWARD_EULER:
168 return PASO_BACKWARD_EULER;
169 default:
170 stringstream temp;
171 temp << "Error - Cannot map option value "<< option << " onto Paso";
172 throw PasoException(temp.str());
173 }
174 }
175
176 int SystemMatrixAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner,
177 const int package, const bool symmetry, Esys_MPIInfo* mpiInfo)
178 {
179 int out=SystemMatrix::getSystemMatrixTypeId(mapOptionToPaso(solver),
180 mapOptionToPaso(preconditioner), mapOptionToPaso(package),
181 symmetry?1:0, mpiInfo);
182 checkPasoError();
183 return out;
184 }
185
186 void SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
187 {
188 int first_row_index = m_system_matrix->row_distribution->first_component[m_system_matrix->mpi_info->rank];
189 int last_row_index = m_system_matrix->row_distribution->first_component[m_system_matrix->mpi_info->rank+1]-1;
190 int first_col_index = m_system_matrix->col_distribution->first_component[m_system_matrix->mpi_info->rank];
191 int last_col_index = m_system_matrix->col_distribution->first_component[m_system_matrix->mpi_info->rank+1]-1;
192
193 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n",
194 m_system_matrix->mpi_info->rank, m_system_matrix->mpi_info->size);
195
196 switch (m_system_matrix->type) {
197 case MATRIX_FORMAT_DEFAULT:
198 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n");
199 break;
200 case MATRIX_FORMAT_CSC:
201 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n");
202 break;
203 case MATRIX_FORMAT_BLK1:
204 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n");
205 break;
206 case MATRIX_FORMAT_OFFSET1:
207 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n");
208 break;
209 case MATRIX_FORMAT_TRILINOS_CRS:
210 fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n");
211 break;
212 default:
213 fprintf(stdout, "\tMatrix type unknown\n");
214 break;
215 }
216
217 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
218 fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
219 fprintf(stdout, "\tmainBlock numRows %d\n", m_system_matrix->mainBlock->numRows);
220 fprintf(stdout, "\tmainBlock numCols %d\n", m_system_matrix->mainBlock->numCols);
221 fprintf(stdout, "\tmainBlock pattern numOutput %d\n", m_system_matrix->mainBlock->pattern->numOutput);
222 fprintf(stdout, "\tcol_coupleBlock numRows %d\n", m_system_matrix->col_coupleBlock->numRows);
223 fprintf(stdout, "\tcol_coupleBlock numCols %d\n", m_system_matrix->col_coupleBlock->numCols);
224 fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", m_system_matrix->col_coupleBlock->pattern->numOutput);
225 fprintf(stdout, "\trow_coupleBlock numRows %d\n", m_system_matrix->row_coupleBlock->numRows);
226 fprintf(stdout, "\trow_coupleBlock numCols %d\n", m_system_matrix->row_coupleBlock->numCols);
227 fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", m_system_matrix->row_coupleBlock->pattern->numOutput);
228 fprintf(stdout, "\trow_block_size %d\n", m_system_matrix->row_block_size);
229 fprintf(stdout, "\tcol_block_size %d\n", m_system_matrix->col_block_size);
230 fprintf(stdout, "\tblock_size %d\n", m_system_matrix->block_size);
231 fprintf(stdout, "\tlogical_row_block_size %d\n", m_system_matrix->logical_row_block_size);
232 fprintf(stdout, "\tlogical_col_block_size %d\n", m_system_matrix->logical_col_block_size);
233
234 }
235
236 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, boost::python::object& options) const
237 {
238 Options paso_options;
239 options.attr("resetDiagnostics")();
240 escriptToPasoOptions(&paso_options,options);
241 if ( out.getDataPointSize() != getColumnBlockSize()) {
242 throw PasoException("solve : column block size does not match the number of components of solution.");
243 } else if ( in.getDataPointSize() != getRowBlockSize()) {
244 throw PasoException("solve : row block size does not match the number of components of right hand side.");
245 } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
246 throw PasoException("solve : column function space and function space of solution don't match.");
247 } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
248 throw PasoException("solve : row function space and function space of right hand side don't match.");
249 }
250 out.expand();
251 in.expand();
252 double* out_dp=out.getSampleDataRW(0);
253 double* in_dp=in.getSampleDataRW(0);
254 paso::solve(m_system_matrix, out_dp, in_dp, &paso_options);
255 pasoToEscriptOptions(&paso_options,options);
256 checkPasoError();
257 }
258
259 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
260 {
261 if ( col_q.getDataPointSize() != getColumnBlockSize()) {
262 throw PasoException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
263 } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
264 throw PasoException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
265 } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
266 throw PasoException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
267 } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
268 throw PasoException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
269 }
270 row_q.expand();
271 col_q.expand();
272 row_q.requireWrite();
273 col_q.requireWrite();
274 double* row_q_dp=row_q.getSampleDataRW(0);
275 double* col_q_dp=col_q.getSampleDataRW(0);
276 m_system_matrix->nullifyRowsAndCols(row_q_dp, col_q_dp, mdv);
277 checkPasoError();
278 }
279
280 void SystemMatrixAdapter::saveMM(const std::string& filename) const
281 {
282 m_system_matrix->saveMM(filename.c_str());
283 checkPasoError();
284 }
285
286 void SystemMatrixAdapter::saveHB(const std::string& filename) const
287 {
288 m_system_matrix->saveHB(filename.c_str());
289 checkPasoError();
290 }
291
292 void SystemMatrixAdapter::resetValues() const
293 {
294 m_system_matrix->setValues(0.);
295 solve_free(m_system_matrix.get());
296 checkPasoError();
297 }
298
299 void SystemMatrixAdapter::pasoToEscriptOptions(const Options* paso_options,
300 boost::python::object& options)
301 {
302 #define SET(__key__,__val__,__type__) options.attr("_updateDiagnostics")(__key__,(__type__)paso_options->__val__)
303 SET("num_iter", num_iter, int);
304 SET("num_level", num_level, int);
305 SET("num_inner_iter", num_inner_iter, int);
306 SET("time", time, double);
307 SET("set_up_time", set_up_time, double);
308 SET("net_time", net_time, double);
309 SET("residual_norm", residual_norm, double);
310 SET("converged",converged, bool);
311 SET("time_step_backtracking_used", time_step_backtracking_used,bool);
312 SET("coarse_level_sparsity",coarse_level_sparsity,double);
313 SET("num_coarse_unknowns",num_coarse_unknowns,int);
314 #undef SET
315 }
316
317 void SystemMatrixAdapter::escriptToPasoOptions(Options* paso_options,
318 const boost::python::object& options)
319 {
320 escript::SolverBuddy sb = boost::python::extract<escript::SolverBuddy>(options);
321
322 paso_options->setDefaults();
323 paso_options->method = mapOptionToPaso(sb.getSolverMethod());
324 paso_options->package = mapOptionToPaso(sb.getPackage());
325 paso_options->verbose = sb.isVerbose();
326 paso_options->symmetric = sb.isSymmetric();
327 paso_options->tolerance = sb.getTolerance();
328 paso_options->absolute_tolerance = sb.getAbsoluteTolerance();
329 paso_options->inner_tolerance = sb.getInnerTolerance();
330 paso_options->adapt_inner_tolerance = sb.adaptInnerTolerance();
331 paso_options->reordering = mapOptionToPaso(sb.getReordering());
332 paso_options->preconditioner = mapOptionToPaso(sb.getPreconditioner());
333 paso_options->iter_max = sb.getIterMax();
334 paso_options->inner_iter_max = sb.getInnerIterMax();
335 paso_options->drop_tolerance = sb.getDropTolerance();
336 paso_options->drop_storage = sb.getDropStorage();
337 paso_options->truncation = sb.getTruncation();
338 paso_options->restart = sb._getRestartForC();
339 paso_options->sweeps = sb.getNumSweeps();
340 paso_options->pre_sweeps = sb.getNumPreSweeps();
341 paso_options->post_sweeps = sb.getNumPostSweeps();
342 paso_options->level_max = sb.getLevelMax();
343 paso_options->min_coarse_matrix_size = sb.getMinCoarseMatrixSize();
344 paso_options->coarsening_threshold = sb.getCoarseningThreshold();
345 paso_options->accept_failed_convergence = sb.acceptConvergenceFailure();
346 paso_options->coarsening_method = mapOptionToPaso(sb.getCoarsening());
347 paso_options->smoother = mapOptionToPaso(sb.getSmoother());
348 paso_options->relaxation_factor = sb.getRelaxationFactor();
349 paso_options->use_local_preconditioner = sb.useLocalPreconditioner();
350 paso_options->min_coarse_sparsity = sb.getMinCoarseMatrixSparsity();
351 paso_options->refinements = sb.getNumRefinements();
352 paso_options->coarse_matrix_refinements = sb.getNumCoarseMatrixRefinements();
353 paso_options->usePanel = sb.usePanel();
354 paso_options->interpolation_method = sb.getAMGInterpolation();
355 paso_options->diagonal_dominance_threshold = sb.getDiagonalDominanceThreshold();
356 paso_options->ode_solver = sb.getODESolver();
357 }
358
359
360 } // end of namespace
361

  ViewVC Help
Powered by ViewVC 1.1.26