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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3051 - (hide annotations)
Tue Jun 29 00:45:49 2010 UTC (9 years, 5 months ago) by lgao
File size: 14146 byte(s)
Hybrid MPI/OpenMP versioned Gauss_Seidel preconditioner is added. 
To use it, use "SolverOptions.GAUSS_SEIDEL_MPI" in python script. 


1 jgs 472
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 ksteube 817 #ifdef PASO_MPI
16     #include <mpi.h>
17     #endif
18 jgs 203 #include "SystemMatrixAdapter.h"
19 jgs 82
20     using namespace std;
21    
22     namespace finley {
23    
24 jgs 102 struct null_deleter
25     {
26 phornby 1628 void operator()(void const *ptr) const
27     {
28     }
29 jgs 102 };
30    
31    
32 jgs 82 SystemMatrixAdapter::SystemMatrixAdapter()
33     {
34     throw FinleyAdapterException("Error - Illegal to generate default SystemMatrixAdapter.");
35     }
36    
37 jgs 150 SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
38 jgs 82 const int row_blocksize,
39     const escript::FunctionSpace& row_functionspace,
40     const int column_blocksize,
41     const escript::FunctionSpace& column_functionspace):
42 jgs 102 AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
43 jgs 82 {
44 phornby 1628 m_system_matrix.reset(system_matrix,null_deleter());
45 jgs 82 }
46    
47     SystemMatrixAdapter::~SystemMatrixAdapter()
48     {
49 phornby 1628 if (m_system_matrix.unique()) {
50     Paso_SystemMatrix* mat=m_system_matrix.get();
51     Paso_SystemMatrix_free(mat);
52     }
53 jgs 82 }
54    
55 jgs 150 Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const
56 jgs 82 {
57     return m_system_matrix.get();
58     }
59    
60 jgs 153 void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const
61 jgs 82 {
62 jgs 150 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
63    
64 phornby 1628 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 jfenwick 2271 x.requireWrite();
76     y.requireWrite();
77     double* x_dp=x.getSampleDataRW(0);
78     double* y_dp=y.getSampleDataRW(0);
79 phornby 1628 Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
80     checkPasoError();
81 jgs 82 }
82    
83 jgs 150 int SystemMatrixAdapter::mapOptionToPaso(const int option) {
84 gross 2474
85 jgs 150 switch (option) {
86 gross 1639 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 gross 2474 case ESCRIPT_LUMPING:
113     return PASO_LUMPING;
114 gross 1639 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 artak 2760 case ESCRIPT_AMLI:
131     return PASO_AMLI;
132 gross 2474 case ESCRIPT_REC_ILU:
133     return PASO_REC_ILU;
134 gross 1639 case ESCRIPT_TRILINOS:
135     return PASO_TRILINOS;
136     case ESCRIPT_NONLINEAR_GMRES:
137     return PASO_NONLINEAR_GMRES;
138 gross 2474 case ESCRIPT_TFQMR :
139 artak 1703 return PASO_TFQMR;
140 artak 1787 case ESCRIPT_MINRES:
141     return PASO_MINRES;
142 gross 2474 case ESCRIPT_GAUSS_SEIDEL:
143     return PASO_GAUSS_SEIDEL;
144 lgao 3051 case ESCRIPT_GAUSS_SEIDEL_MPI:
145     return PASO_GAUSS_SEIDEL_MPI;
146 gross 2474 case ESCRIPT_RILU:
147     return PASO_RILU;
148     case ESCRIPT_DEFAULT_REORDERING:
149     return PASO_DEFAULT_REORDERING;
150     case ESCRIPT_SUPER_LU:
151     return PASO_SUPER_LU;
152     case ESCRIPT_PASTIX:
153     return PASO_PASTIX;
154     case ESCRIPT_YAIR_SHAPIRA_COARSENING:
155     return PASO_YAIR_SHAPIRA_COARSENING;
156     case ESCRIPT_RUGE_STUEBEN_COARSENING:
157     return PASO_RUGE_STUEBEN_COARSENING;
158 artak 2816 case ESCRIPT_STANDARD_COARSENING:
159     return PASO_STANDARD_COARSENING;
160 gross 2474 case ESCRIPT_AGGREGATION_COARSENING:
161     return PASO_AGGREGATION_COARSENING;
162     case ESCRIPT_NO_PRECONDITIONER:
163     return PASO_NO_PRECONDITIONER;
164 artak 2524 case ESCRIPT_MIN_COARSE_MATRIX_SIZE:
165     return PASO_MIN_COARSE_MATRIX_SIZE;
166 gross 1639 default:
167     stringstream temp;
168     temp << "Error - Cannot map option value "<< option << " onto Paso";
169     throw FinleyAdapterException(temp.str());
170     }
171 jgs 150 }
172    
173 ksteube 1339 void finley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
174     {
175 phornby 1628 Paso_SystemMatrix* mat=m_system_matrix.get();
176     int first_row_index = mat->row_distribution->first_component[mat->mpi_info->rank];
177     int last_row_index = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
178     int first_col_index = mat->col_distribution->first_component[mat->mpi_info->rank];
179     int last_col_index = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;
180 ksteube 1339
181 phornby 1628 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
182 ksteube 1339
183 phornby 1628 switch (mat->type) {
184     case MATRIX_FORMAT_DEFAULT: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n"); break;
185     case MATRIX_FORMAT_CSC: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n"); break;
186     case MATRIX_FORMAT_SYM: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_SYM\n"); break;
187     case MATRIX_FORMAT_BLK1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n"); break;
188     case MATRIX_FORMAT_OFFSET1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n"); break;
189     case MATRIX_FORMAT_TRILINOS_CRS: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n"); break;
190     default: fprintf(stdout, "\tMatrix type unknown\n"); break;
191     }
192 ksteube 1339
193 phornby 1628 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
194     fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
195     fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
196     fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
197     fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
198     fprintf(stdout, "\tcol_coupleBlock numRows %d\n", mat->col_coupleBlock->numRows);
199     fprintf(stdout, "\tcol_coupleBlock numCols %d\n", mat->col_coupleBlock->numCols);
200     fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", mat->col_coupleBlock->pattern->numOutput);
201     fprintf(stdout, "\trow_coupleBlock numRows %d\n", mat->row_coupleBlock->numRows);
202     fprintf(stdout, "\trow_coupleBlock numCols %d\n", mat->row_coupleBlock->numCols);
203     fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", mat->row_coupleBlock->pattern->numOutput);
204     fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
205     fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
206     fprintf(stdout, "\tblock_size %d\n", mat->block_size);
207     fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
208     fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
209 ksteube 1339
210     }
211    
212 gross 2474 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, boost::python::object& options) const
213 jgs 82 {
214 phornby 1628 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
215     Paso_Options paso_options;
216 gross 2474 options.attr("resetDiagnostics")();
217     escriptToPasoOptions(&paso_options,options);
218 phornby 1628 if ( out.getDataPointSize() != getColumnBlockSize()) {
219     throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
220     } else if ( in.getDataPointSize() != getRowBlockSize()) {
221     throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
222     } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
223     throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
224     } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
225     throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
226     }
227     out.expand();
228     in.expand();
229 jfenwick 2271 double* out_dp=out.getSampleDataRW(0);
230     double* in_dp=in.getSampleDataRW(0);
231 phornby 1628 Paso_solve(mat,out_dp,in_dp,&paso_options);
232 gross 2474 pasoToEscriptOptions(&paso_options,options);
233 phornby 1628 checkPasoError();
234 jgs 82 }
235    
236 jgs 153 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
237 jgs 82 {
238 phornby 1628 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
239     if ( col_q.getDataPointSize() != getColumnBlockSize()) {
240     throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
241     } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
242     throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
243     } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
244     throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
245     } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
246     throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
247     }
248     row_q.expand();
249     col_q.expand();
250 jfenwick 2271 row_q.requireWrite();
251     col_q.requireWrite();
252     double* row_q_dp=row_q.getSampleDataRW(0);
253     double* col_q_dp=col_q.getSampleDataRW(0);
254 phornby 1628 Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
255     checkPasoError();
256 jgs 82 }
257    
258 jgs 102 void SystemMatrixAdapter::saveMM(const std::string& fileName) const
259     {
260 phornby 1628 if( fileName.size() == 0 )
261     {
262     throw FinleyAdapterException("Null file name!");
263     }
264    
265     char *fName = TMPMEMALLOC(fileName.size()+1,char);
266 woo409 757
267 phornby 1628 strcpy(fName,fileName.c_str());
268     Paso_SystemMatrix* mat = getPaso_SystemMatrix();
269     Paso_SystemMatrix_saveMM(mat,fName);
270     checkPasoError();
271     TMPMEMFREE(fName);
272 woo409 757
273 jgs 102 }
274    
275 jgs 123 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
276     {
277 phornby 1628 if( fileName.size() == 0 )
278     {
279     throw FinleyAdapterException("Null file name!");
280     }
281 woo409 757
282 phornby 1628 char *fName = TMPMEMALLOC(fileName.size()+1,char);
283 woo409 757
284 phornby 1628 strcpy(fName,fileName.c_str());
285     Paso_SystemMatrix* mat = getPaso_SystemMatrix();
286     Paso_SystemMatrix_saveHB(mat,fName);
287     checkPasoError();
288     TMPMEMFREE(fName);
289    
290 jgs 123 }
291    
292 jgs 149 void SystemMatrixAdapter::resetValues() const
293 jgs 102 {
294 jgs 150 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
295     Paso_SystemMatrix_setValues(mat,0.);
296     Paso_solve_free(mat);
297     checkPasoError();
298 jgs 108 }
299 jgs 102
300 gross 2474 void SystemMatrixAdapter::pasoToEscriptOptions(const Paso_Options* paso_options,boost::python::object& options)
301 gross 1364 {
302 gross 2474 #define SET(__key__,__val__,__type__) options.attr("_updateDiagnostics")(__key__,(__type__)paso_options->__val__)
303     SET("num_iter", num_iter, int);
304     SET("num_level", num_level, int);
305     SET("num_inner_iter", num_inner_iter, int);
306     SET("time", time, double);
307     SET("set_up_time", set_up_time, double);
308 lgao 2596 SET("net_time", net_time, double);
309 gross 2474 SET("residual_norm", residual_norm, double);
310     SET("converged",converged, bool);
311     #undef SET
312     }
313     void SystemMatrixAdapter::escriptToPasoOptions(Paso_Options* paso_options, const boost::python::object& options)
314     {
315     #define EXTRACT(__key__,__val__,__type__) paso_options->__val__=boost::python::extract<__type__>(options.attr(__key__)())
316     #define EXTRACT_OPTION(__key__,__val__,__type__) paso_options->__val__=mapOptionToPaso(boost::python::extract<__type__>(options.attr(__key__)()))
317    
318 phornby 1628 Paso_Options_setDefaults(paso_options);
319 gross 2474 EXTRACT_OPTION("getSolverMethod",method,index_t);
320     EXTRACT_OPTION("getPackage",package,index_t);
321     EXTRACT("isVerbose",verbose,bool_t);
322     EXTRACT("isSymmetric",symmetric,bool_t);
323     EXTRACT("getTolerance",tolerance, double);
324     EXTRACT("getAbsoluteTolerance",absolute_tolerance, double);
325     EXTRACT("getInnerTolerance",inner_tolerance, double);
326     EXTRACT("adaptInnerTolerance",adapt_inner_tolerance, bool_t);
327     EXTRACT_OPTION("getReordering", reordering, index_t);
328     EXTRACT_OPTION("getPreconditioner", preconditioner, index_t);
329     EXTRACT("getIterMax", iter_max, dim_t);
330     EXTRACT("getInnerIterMax", inner_iter_max, dim_t);
331     EXTRACT("getDropTolerance", drop_tolerance, double);
332     EXTRACT("getDropStorage", drop_storage, double);
333     EXTRACT("getTruncation", truncation, dim_t);
334     EXTRACT("_getRestartForC", restart, index_t);
335     EXTRACT("getNumSweeps", sweeps, index_t);
336     EXTRACT("getNumPreSweeps", pre_sweeps, dim_t);
337     EXTRACT("getNumPostSweeps", post_sweeps, dim_t);
338     EXTRACT("getLevelMax", level_max, dim_t);
339 artak 2524 EXTRACT("getMinCoarseMatrixSize", min_coarse_matrix_size, dim_t);
340 gross 2474 EXTRACT("getCoarseningThreshold", coarsening_threshold, double);
341     EXTRACT("acceptConvergenceFailure", accept_failed_convergence, bool_t);
342 artak 2828 EXTRACT_OPTION("getCoarsening", coarsening_method, index_t);
343     EXTRACT_OPTION("getSmoother", smoother, index_t);
344 gross 2474 EXTRACT("getRelaxationFactor", relaxation_factor, double);
345    
346 phornby 1628 #undef EXTRACT
347     #undef EXTRACT_OPTION
348 gross 1364 }
349 phornby 1628
350 gross 1364
351 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26