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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26