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

Contents of /branches/symbolic_from_3470/dudley/src/CPPAdapter/SystemMatrixAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3517 - (show annotations)
Fri May 20 01:16:41 2011 UTC (7 years, 11 months ago) by caltinay
File size: 15162 byte(s)
Merged up to and including revision 3514 from trunk and implemented
symbolic inverse.

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 case ESCRIPT_CLASSIC_INTERPOLATION_WITH_FF_COUPLING:
157 return PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING;
158 case ESCRIPT_CLASSIC_INTERPOLATION:
159 return PASO_CLASSIC_INTERPOLATION;
160 case ESCRIPT_DIRECT_INTERPOLATION:
161 return PASO_DIRECT_INTERPOLATION;
162 case ESCRIPT_BOOMERAMG:
163 return PASO_BOOMERAMG;
164 case ESCRIPT_CIJP_FIXED_RANDOM_COARSENING:
165 return PASO_CIJP_FIXED_RANDOM_COARSENING;
166 case ESCRIPT_CIJP_COARSENING:
167 return PASO_CIJP_COARSENING;
168 case ESCRIPT_FALGOUT_COARSENING:
169 return PASO_FALGOUT_COARSENING;
170 case ESCRIPT_PMIS_COARSENING:
171 return PASO_PMIS_COARSENING;
172 case ESCRIPT_HMIS_COARSENING:
173 return PASO_HMIS_COARSENING;
174
175 default:
176 stringstream temp;
177 temp << "Error - Cannot map option value "<< option << " onto Paso";
178 throw DudleyAdapterException(temp.str());
179 }
180 }
181
182 void dudley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
183 {
184 Paso_SystemMatrix* mat=m_system_matrix.get();
185 int first_row_index = mat->row_distribution->first_component[mat->mpi_info->rank];
186 int last_row_index = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
187 int first_col_index = mat->col_distribution->first_component[mat->mpi_info->rank];
188 int last_col_index = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;
189
190 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
191
192 switch (mat->type) {
193 case MATRIX_FORMAT_DEFAULT: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n"); break;
194 case MATRIX_FORMAT_CSC: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n"); break;
195 case MATRIX_FORMAT_BLK1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n"); break;
196 case MATRIX_FORMAT_OFFSET1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n"); break;
197 case MATRIX_FORMAT_TRILINOS_CRS: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n"); break;
198 default: fprintf(stdout, "\tMatrix type unknown\n"); break;
199 }
200
201 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
202 fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
203 fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
204 fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
205 fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
206 fprintf(stdout, "\tcol_coupleBlock numRows %d\n", mat->col_coupleBlock->numRows);
207 fprintf(stdout, "\tcol_coupleBlock numCols %d\n", mat->col_coupleBlock->numCols);
208 fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", mat->col_coupleBlock->pattern->numOutput);
209 fprintf(stdout, "\trow_coupleBlock numRows %d\n", mat->row_coupleBlock->numRows);
210 fprintf(stdout, "\trow_coupleBlock numCols %d\n", mat->row_coupleBlock->numCols);
211 fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", mat->row_coupleBlock->pattern->numOutput);
212 fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
213 fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
214 fprintf(stdout, "\tblock_size %d\n", mat->block_size);
215 fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
216 fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
217
218 }
219
220 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, boost::python::object& options) const
221 {
222 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
223 Paso_Options paso_options;
224 options.attr("resetDiagnostics")();
225 escriptToPasoOptions(&paso_options,options);
226 if ( out.getDataPointSize() != getColumnBlockSize()) {
227 throw DudleyAdapterException("solve : column block size does not match the number of components of solution.");
228 } else if ( in.getDataPointSize() != getRowBlockSize()) {
229 throw DudleyAdapterException("solve : row block size does not match the number of components of right hand side.");
230 } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
231 throw DudleyAdapterException("solve : column function space and function space of solution don't match.");
232 } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
233 throw DudleyAdapterException("solve : row function space and function space of right hand side don't match.");
234 }
235 out.expand();
236 in.expand();
237 double* out_dp=out.getSampleDataRW(0);
238 double* in_dp=in.getSampleDataRW(0);
239 Paso_solve(mat,out_dp,in_dp,&paso_options);
240 pasoToEscriptOptions(&paso_options,options);
241 checkPasoError();
242 }
243
244 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
245 {
246 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
247 if ( col_q.getDataPointSize() != getColumnBlockSize()) {
248 throw DudleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
249 } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
250 throw DudleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
251 } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
252 throw DudleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
253 } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
254 throw DudleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
255 }
256 row_q.expand();
257 col_q.expand();
258 row_q.requireWrite();
259 col_q.requireWrite();
260 double* row_q_dp=row_q.getSampleDataRW(0);
261 double* col_q_dp=col_q.getSampleDataRW(0);
262 Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
263 checkPasoError();
264 }
265
266 void SystemMatrixAdapter::saveMM(const std::string& fileName) const
267 {
268 if( fileName.size() == 0 )
269 {
270 throw DudleyAdapterException("Null file name!");
271 }
272
273 char *fName = TMPMEMALLOC(fileName.size()+1,char);
274
275 strcpy(fName,fileName.c_str());
276 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
277 Paso_SystemMatrix_saveMM(mat,fName);
278 checkPasoError();
279 TMPMEMFREE(fName);
280
281 }
282
283 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
284 {
285 if( fileName.size() == 0 )
286 {
287 throw DudleyAdapterException("Null file name!");
288 }
289
290 char *fName = TMPMEMALLOC(fileName.size()+1,char);
291
292 strcpy(fName,fileName.c_str());
293 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
294 Paso_SystemMatrix_saveHB(mat,fName);
295 checkPasoError();
296 TMPMEMFREE(fName);
297
298 }
299
300 void SystemMatrixAdapter::resetValues() const
301 {
302 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
303 Paso_SystemMatrix_setValues(mat,0.);
304 Paso_solve_free(mat);
305 checkPasoError();
306 }
307
308 void SystemMatrixAdapter::pasoToEscriptOptions(const Paso_Options* paso_options,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 void SystemMatrixAdapter::escriptToPasoOptions(Paso_Options* paso_options, const boost::python::object& options)
325 {
326 #define EXTRACT(__key__,__val__,__type__) paso_options->__val__=boost::python::extract<__type__>(options.attr(__key__)())
327 #define EXTRACT_OPTION(__key__,__val__,__type__) paso_options->__val__=mapOptionToPaso(boost::python::extract<__type__>(options.attr(__key__)()))
328
329 Paso_Options_setDefaults(paso_options);
330 EXTRACT_OPTION("getSolverMethod",method,index_t);
331 EXTRACT_OPTION("getPackage",package,index_t);
332 EXTRACT("isVerbose",verbose,bool_t);
333 EXTRACT("isSymmetric",symmetric,bool_t);
334 EXTRACT("getTolerance",tolerance, double);
335 EXTRACT("getAbsoluteTolerance",absolute_tolerance, double);
336 EXTRACT("getInnerTolerance",inner_tolerance, double);
337 EXTRACT("adaptInnerTolerance",adapt_inner_tolerance, bool_t);
338 EXTRACT_OPTION("getReordering", reordering, index_t);
339 EXTRACT_OPTION("getPreconditioner", preconditioner, index_t);
340 EXTRACT("getIterMax", iter_max, dim_t);
341 EXTRACT("getInnerIterMax", inner_iter_max, dim_t);
342 EXTRACT("getDropTolerance", drop_tolerance, double);
343 EXTRACT("getDropStorage", drop_storage, double);
344 EXTRACT("getTruncation", truncation, dim_t);
345 EXTRACT("_getRestartForC", restart, index_t);
346 EXTRACT("getNumSweeps", sweeps, index_t);
347 EXTRACT("getNumPreSweeps", pre_sweeps, dim_t);
348 EXTRACT("getNumPostSweeps", post_sweeps, dim_t);
349 EXTRACT("getLevelMax", level_max, dim_t);
350 EXTRACT("getMinCoarseMatrixSize", min_coarse_matrix_size, dim_t);
351 EXTRACT("getCoarseningThreshold", coarsening_threshold, double);
352 EXTRACT("acceptConvergenceFailure", accept_failed_convergence, bool_t);
353 EXTRACT_OPTION("getCoarsening", coarsening_method, index_t);
354 EXTRACT_OPTION("getSmoother", smoother, index_t);
355 EXTRACT("getRelaxationFactor", relaxation_factor, double);
356 EXTRACT("useLocalPreconditioner", use_local_preconditioner, bool_t);
357 EXTRACT("getMinCoarseMatrixSparsity",min_coarse_sparsity, double);
358 EXTRACT("getNumRefinements",refinements, dim_t);
359 EXTRACT("getNumCoarseMatrixRefinements",coarse_matrix_refinements, dim_t);
360 EXTRACT("usePanel",usePanel, bool_t);
361 EXTRACT("getAMGInterpolation", interpolation_method, index_t);
362 EXTRACT("getDiagonalDominanceThreshold", diagonal_dominance_threshold, double);
363
364 #undef EXTRACT
365 #undef EXTRACT_OPTION
366 }
367
368
369 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26