/[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 1364 - (show annotations)
Mon Dec 17 07:22:45 2007 UTC (13 years, 6 months ago) by gross
File size: 12482 byte(s)
finley interface to paso's transport solver added.
1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #ifdef PASO_MPI
17 #include <mpi.h>
18 #endif
19 #include "SystemMatrixAdapter.h"
20
21 using namespace std;
22
23 namespace finley {
24
25 struct null_deleter
26 {
27 void operator()(void const *ptr) const
28 {
29 }
30 };
31
32
33 SystemMatrixAdapter::SystemMatrixAdapter()
34 {
35 throw FinleyAdapterException("Error - Illegal to generate default SystemMatrixAdapter.");
36 }
37
38 SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
39 const int row_blocksize,
40 const escript::FunctionSpace& row_functionspace,
41 const int column_blocksize,
42 const escript::FunctionSpace& column_functionspace):
43 AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
44 {
45 m_system_matrix.reset(system_matrix,null_deleter());
46 }
47
48 SystemMatrixAdapter::~SystemMatrixAdapter()
49 {
50 if (m_system_matrix.unique()) {
51 Paso_SystemMatrix* mat=m_system_matrix.get();
52 Paso_SystemMatrix_free(mat);
53 }
54 }
55
56 Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const
57 {
58 return m_system_matrix.get();
59 }
60
61 void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const
62 {
63 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
64
65 if ( x.getDataPointSize() != getColumnBlockSize()) {
66 throw FinleyAdapterException("matrix vector product : column block size does not match the number of components in input.");
67 } else if (y.getDataPointSize() != getRowBlockSize()) {
68 throw FinleyAdapterException("matrix vector product : row block size does not match the number of components in output.");
69 } else if ( x.getFunctionSpace() != getColumnFunctionSpace()) {
70 throw FinleyAdapterException("matrix vector product : column function space and function space of input don't match.");
71 } else if (y.getFunctionSpace() != getRowFunctionSpace()) {
72 throw FinleyAdapterException("matrix vector product : row function space and function space of output don't match.");
73 }
74 x.expand();
75 y.expand();
76 double* x_dp=x.getSampleData(0);
77 double* y_dp=y.getSampleData(0);
78 Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
79 checkPasoError();
80 }
81
82 int SystemMatrixAdapter::mapOptionToPaso(const int option) {
83 switch (option) {
84 case ESCRIPT_DEFAULT:
85 return PASO_DEFAULT;
86 case ESCRIPT_DIRECT:
87 return PASO_DIRECT;
88 case ESCRIPT_CHOLEVSKY:
89 return PASO_CHOLEVSKY;
90 case ESCRIPT_PCG:
91 return PASO_PCG;
92 case ESCRIPT_CR:
93 return PASO_CR;
94 case ESCRIPT_CGS:
95 return PASO_CGS;
96 case ESCRIPT_BICGSTAB:
97 return PASO_BICGSTAB;
98 case ESCRIPT_SSOR:
99 return PASO_SSOR;
100 case ESCRIPT_ILU0:
101 return PASO_ILU0;
102 case ESCRIPT_ILUT:
103 return PASO_ILUT;
104 case ESCRIPT_JACOBI:
105 return PASO_JACOBI;
106 case ESCRIPT_GMRES:
107 return PASO_GMRES;
108 case ESCRIPT_PRES20:
109 return PASO_PRES20;
110 case ESCRIPT_NO_REORDERING:
111 return PASO_NO_REORDERING;
112 case ESCRIPT_MINIMUM_FILL_IN:
113 return PASO_MINIMUM_FILL_IN;
114 case ESCRIPT_NESTED_DISSECTION:
115 return PASO_NESTED_DISSECTION;
116 case ESCRIPT_SCSL:
117 return PASO_SCSL;
118 case ESCRIPT_MKL:
119 return PASO_MKL;
120 case ESCRIPT_UMFPACK:
121 return PASO_UMFPACK;
122 case ESCRIPT_ITERATIVE:
123 return PASO_ITERATIVE;
124 case ESCRIPT_PASO:
125 return PASO_PASO;
126 case ESCRIPT_LUMPING:
127 return PASO_LUMPING;
128 case ESCRIPT_AMG:
129 return PASO_AMG;
130 case ESCRIPT_RILU:
131 return PASO_RILU;
132 case ESCRIPT_TRILINOS:
133 return PASO_TRILINOS;
134 default:
135 stringstream temp;
136 temp << "Error - Cannot map option value "<< option << " onto Paso";
137 throw FinleyAdapterException(temp.str());
138 }
139 }
140
141 void finley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
142 {
143 Paso_SystemMatrix* mat=m_system_matrix.get();
144 int first_row_index = mat->row_distribution->first_component[mat->mpi_info->rank];
145 int last_row_index = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
146 int first_col_index = mat->col_distribution->first_component[mat->mpi_info->rank];
147 int last_col_index = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;
148
149 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
150
151 switch (mat->type) {
152 case MATRIX_FORMAT_DEFAULT: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n"); break;
153 case MATRIX_FORMAT_CSC: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n"); break;
154 case MATRIX_FORMAT_SYM: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_SYM\n"); break;
155 case MATRIX_FORMAT_BLK1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n"); break;
156 case MATRIX_FORMAT_OFFSET1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n"); break;
157 case MATRIX_FORMAT_TRILINOS_CRS: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n"); break;
158 default: fprintf(stdout, "\tMatrix type unknown\n"); break;
159 }
160
161 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
162 fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
163 fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
164 fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
165 fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
166 fprintf(stdout, "\tcoupleBlock numRows %d\n", mat->coupleBlock->numRows);
167 fprintf(stdout, "\tcoupleBlock numCols %d\n", mat->coupleBlock->numCols);
168 fprintf(stdout, "\tcoupleBlock pattern numOutput %d\n", mat->coupleBlock->pattern->numOutput);
169 fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
170 fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
171 fprintf(stdout, "\tblock_size %d\n", mat->block_size);
172 fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
173 fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
174 fprintf(stdout, "\tlogical_block_size %d\n", mat->logical_block_size);
175
176 if (full) {
177 printf("\trow_distribution: ");
178 for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->row_distribution[i]);
179 printf("\n");
180 printf("\tcol_distribution: ");
181 for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->col_distribution[i]);
182 printf("\n");
183 }
184
185 }
186
187 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
188 {
189 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
190 Paso_Options paso_options;
191 dictToPasoOptions(&paso_options,options);
192 // Paso_Options_setDefaults(&paso_options);
193 // extract options
194 // #define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=boost::python::extract<__type__>(options.get(__key__))
195 // #define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
196 // EXTRACT("verbose",verbose,int);
197 // EXTRACT_OPTION("reordering",reordering,int);
198 // EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
199 // EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
200 // EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
201 // EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
202 // EXTRACT_OPTION("preconditioner",preconditioner,int);
203 // EXTRACT("iter_max",iter_max,int);
204 // EXTRACT("drop_tolerance",drop_tolerance,double);
205 // EXTRACT("drop_storage",drop_storage,double);
206 // EXTRACT("truncation",truncation,int);
207 // EXTRACT("restart",restart,int);
208 // #undef EXTRACT
209 // #undef EXTRACT_OPTION
210 if ( out.getDataPointSize() != getColumnBlockSize()) {
211 throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
212 } else if ( in.getDataPointSize() != getRowBlockSize()) {
213 throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
214 } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
215 throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
216 } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
217 throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
218 }
219 out.expand();
220 in.expand();
221 double* out_dp=out.getSampleData(0);
222 double* in_dp=in.getSampleData(0);
223 Paso_solve(mat,out_dp,in_dp,&paso_options);
224 checkPasoError();
225 }
226
227 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
228 {
229 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
230 if ( col_q.getDataPointSize() != getColumnBlockSize()) {
231 throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
232 } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
233 throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
234 } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
235 throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
236 } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
237 throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
238 }
239 row_q.expand();
240 col_q.expand();
241 double* row_q_dp=row_q.getSampleData(0);
242 double* col_q_dp=col_q.getSampleData(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 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
250
251 strcpy(fName,fileName.c_str());
252 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
253 Paso_SystemMatrix_saveMM(mat,fName);
254 checkPasoError();
255 TMPMEMFREE(fName);
256
257 }
258
259 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
260 {
261 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
262
263 strcpy(fName,fileName.c_str());
264 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
265 Paso_SystemMatrix_saveHB(mat,fName);
266 checkPasoError();
267 TMPMEMFREE(fName);
268
269 }
270
271 void SystemMatrixAdapter::resetValues() const
272 {
273 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
274 Paso_SystemMatrix_setValues(mat,0.);
275 Paso_solve_free(mat);
276 checkPasoError();
277 }
278
279 void SystemMatrixAdapter::dictToPasoOptions(Paso_Options* paso_options, const boost::python::dict& options)
280 {
281 Paso_Options_setDefaults(paso_options);
282 #define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options->__val__=boost::python::extract<__type__>(options.get(__key__))
283 #define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options->__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
284 EXTRACT("verbose",verbose,int);
285 EXTRACT_OPTION("reordering",reordering,int);
286 EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
287 EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
288 EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
289 EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
290 EXTRACT_OPTION("preconditioner",preconditioner,int);
291 EXTRACT("iter_max",iter_max,int);
292 EXTRACT("drop_tolerance",drop_tolerance,double);
293 EXTRACT("drop_storage",drop_storage,double);
294 EXTRACT("truncation",truncation,int);
295 EXTRACT("restart",restart,int);
296 #undef EXTRACT
297 #undef EXTRACT_OPTION
298 }
299
300
301 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26