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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3676 - (hide annotations)
Thu Nov 17 03:52:25 2011 UTC (7 years, 10 months ago) by jfenwick
File size: 15017 byte(s)
Removed some debug text
Corrected some dudley where they should not have been.
Set savanna to build shared for now

1 jfenwick 3661
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 paso {
19    
20     struct null_deleter
21     {
22     void operator()(void const *ptr) const
23     {
24 jfenwick 3676 }
25 jfenwick 3661 };
26    
27    
28     SystemMatrixAdapter::SystemMatrixAdapter()
29     {
30     throw PasoException("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 PasoException("matrix vector product : column block size does not match the number of components in input.");
62     } else if (y.getDataPointSize() != getRowBlockSize()) {
63     throw PasoException("matrix vector product : row block size does not match the number of components in output.");
64     } else if ( x.getFunctionSpace() != getColumnFunctionSpace()) {
65     throw PasoException("matrix vector product : column function space and function space of input don't match.");
66     } else if (y.getFunctionSpace() != getRowFunctionSpace()) {
67     throw PasoException("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 PasoException(temp.str());
179     }
180     }
181    
182     void paso::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 PasoException("solve : column block size does not match the number of components of solution.");
228     } else if ( in.getDataPointSize() != getRowBlockSize()) {
229     throw PasoException("solve : row block size does not match the number of components of right hand side.");
230     } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
231     throw PasoException("solve : column function space and function space of solution don't match.");
232     } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
233     throw PasoException("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 PasoException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
249     } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
250     throw PasoException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
251     } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
252     throw PasoException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
253     } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
254     throw PasoException("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 PasoException("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 PasoException("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

  ViewVC Help
Powered by ViewVC 1.1.26