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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26