/[escript]/trunk/dudley/src/CPPAdapter/SystemMatrixAdapter.cpp
ViewVC logotype

Contents of /trunk/dudley/src/CPPAdapter/SystemMatrixAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3317 - (show annotations)
Thu Oct 28 00:50:41 2010 UTC (9 years ago) by caltinay
File size: 13933 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26