/[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 1823 - (show annotations)
Wed Oct 1 05:56:05 2008 UTC (11 years, 11 months ago) by artak
File size: 12086 byte(s)
GS preconditioner now works for block_size 2 and 3 as well. Sweep parameter introdiced, but not working yet
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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 double* x_dp=x.getSampleData(0);
76 double* y_dp=y.getSampleData(0);
77 Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
78 checkPasoError();
79 }
80
81 int SystemMatrixAdapter::mapOptionToPaso(const int option) {
82 switch (option) {
83 case ESCRIPT_DEFAULT:
84 return PASO_DEFAULT;
85 case ESCRIPT_DIRECT:
86 return PASO_DIRECT;
87 case ESCRIPT_CHOLEVSKY:
88 return PASO_CHOLEVSKY;
89 case ESCRIPT_PCG:
90 return PASO_PCG;
91 case ESCRIPT_CR:
92 return PASO_CR;
93 case ESCRIPT_CGS:
94 return PASO_CGS;
95 case ESCRIPT_BICGSTAB:
96 return PASO_BICGSTAB;
97 case ESCRIPT_SSOR:
98 return PASO_SSOR;
99 case ESCRIPT_ILU0:
100 return PASO_ILU0;
101 case ESCRIPT_ILUT:
102 return PASO_ILUT;
103 case ESCRIPT_JACOBI:
104 return PASO_JACOBI;
105 case ESCRIPT_GMRES:
106 return PASO_GMRES;
107 case ESCRIPT_PRES20:
108 return PASO_PRES20;
109 case ESCRIPT_NO_REORDERING:
110 return PASO_NO_REORDERING;
111 case ESCRIPT_MINIMUM_FILL_IN:
112 return PASO_MINIMUM_FILL_IN;
113 case ESCRIPT_NESTED_DISSECTION:
114 return PASO_NESTED_DISSECTION;
115 case ESCRIPT_SCSL:
116 return PASO_SCSL;
117 case ESCRIPT_MKL:
118 return PASO_MKL;
119 case ESCRIPT_UMFPACK:
120 return PASO_UMFPACK;
121 case ESCRIPT_ITERATIVE:
122 return PASO_ITERATIVE;
123 case ESCRIPT_PASO:
124 return PASO_PASO;
125 case ESCRIPT_LUMPING:
126 return PASO_LUMPING;
127 case ESCRIPT_AMG:
128 return PASO_AMG;
129 case ESCRIPT_RILU:
130 return PASO_RILU;
131 case ESCRIPT_TRILINOS:
132 return PASO_TRILINOS;
133 case ESCRIPT_NONLINEAR_GMRES:
134 return PASO_NONLINEAR_GMRES;
135 case ESCRIPT_TFQMR:
136 return PASO_TFQMR;
137 case ESCRIPT_MINRES:
138 return PASO_MINRES;
139 case ESCRIPT_GS:
140 return PASO_GS;
141 default:
142 stringstream temp;
143 temp << "Error - Cannot map option value "<< option << " onto Paso";
144 throw FinleyAdapterException(temp.str());
145 }
146 }
147
148 void finley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
149 {
150 Paso_SystemMatrix* mat=m_system_matrix.get();
151 int first_row_index = mat->row_distribution->first_component[mat->mpi_info->rank];
152 int last_row_index = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
153 int first_col_index = mat->col_distribution->first_component[mat->mpi_info->rank];
154 int last_col_index = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;
155
156 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
157
158 switch (mat->type) {
159 case MATRIX_FORMAT_DEFAULT: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n"); break;
160 case MATRIX_FORMAT_CSC: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n"); break;
161 case MATRIX_FORMAT_SYM: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_SYM\n"); break;
162 case MATRIX_FORMAT_BLK1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n"); break;
163 case MATRIX_FORMAT_OFFSET1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n"); break;
164 case MATRIX_FORMAT_TRILINOS_CRS: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n"); break;
165 default: fprintf(stdout, "\tMatrix type unknown\n"); break;
166 }
167
168 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
169 fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
170 fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
171 fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
172 fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
173 fprintf(stdout, "\tcol_coupleBlock numRows %d\n", mat->col_coupleBlock->numRows);
174 fprintf(stdout, "\tcol_coupleBlock numCols %d\n", mat->col_coupleBlock->numCols);
175 fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", mat->col_coupleBlock->pattern->numOutput);
176 fprintf(stdout, "\trow_coupleBlock numRows %d\n", mat->row_coupleBlock->numRows);
177 fprintf(stdout, "\trow_coupleBlock numCols %d\n", mat->row_coupleBlock->numCols);
178 fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", mat->row_coupleBlock->pattern->numOutput);
179 fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
180 fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
181 fprintf(stdout, "\tblock_size %d\n", mat->block_size);
182 fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
183 fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
184 fprintf(stdout, "\tlogical_block_size %d\n", mat->logical_block_size);
185
186 if (full) {
187 printf("\trow_distribution: ");
188 for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->row_distribution[i]);
189 printf("\n");
190 printf("\tcol_distribution: ");
191 for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->col_distribution[i]);
192 printf("\n");
193 }
194
195 }
196
197 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
198 {
199 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
200 Paso_Options paso_options;
201 dictToPasoOptions(&paso_options,options);
202 if ( out.getDataPointSize() != getColumnBlockSize()) {
203 throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
204 } else if ( in.getDataPointSize() != getRowBlockSize()) {
205 throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
206 } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
207 throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
208 } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
209 throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
210 }
211 out.expand();
212 in.expand();
213 double* out_dp=out.getSampleData(0);
214 double* in_dp=in.getSampleData(0);
215 Paso_solve(mat,out_dp,in_dp,&paso_options);
216 checkPasoError();
217 }
218
219 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
220 {
221 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
222 if ( col_q.getDataPointSize() != getColumnBlockSize()) {
223 throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
224 } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
225 throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
226 } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
227 throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
228 } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
229 throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
230 }
231 row_q.expand();
232 col_q.expand();
233 double* row_q_dp=row_q.getSampleData(0);
234 double* col_q_dp=col_q.getSampleData(0);
235 Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
236 checkPasoError();
237 }
238
239 void SystemMatrixAdapter::saveMM(const std::string& fileName) const
240 {
241 if( fileName.size() == 0 )
242 {
243 throw FinleyAdapterException("Null file name!");
244 }
245
246 char *fName = TMPMEMALLOC(fileName.size()+1,char);
247
248 strcpy(fName,fileName.c_str());
249 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
250 Paso_SystemMatrix_saveMM(mat,fName);
251 checkPasoError();
252 TMPMEMFREE(fName);
253
254 }
255
256 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
257 {
258 if( fileName.size() == 0 )
259 {
260 throw FinleyAdapterException("Null file name!");
261 }
262
263 char *fName = TMPMEMALLOC(fileName.size()+1,char);
264
265 strcpy(fName,fileName.c_str());
266 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
267 Paso_SystemMatrix_saveHB(mat,fName);
268 checkPasoError();
269 TMPMEMFREE(fName);
270
271 }
272
273 void SystemMatrixAdapter::resetValues() const
274 {
275 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
276 Paso_SystemMatrix_setValues(mat,0.);
277 Paso_solve_free(mat);
278 checkPasoError();
279 }
280
281 void SystemMatrixAdapter::dictToPasoOptions(Paso_Options* paso_options, const boost::python::dict& options)
282 {
283 Paso_Options_setDefaults(paso_options);
284 #define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options->__val__=boost::python::extract<__type__>(options.get(__key__))
285 #define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options->__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
286 EXTRACT("verbose",verbose,int);
287 EXTRACT_OPTION("reordering",reordering,int);
288 EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
289 EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
290 EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
291 EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
292 EXTRACT_OPTION("preconditioner",preconditioner,int);
293 EXTRACT("iter_max",iter_max,int);
294 EXTRACT("drop_tolerance",drop_tolerance,double);
295 EXTRACT("drop_storage",drop_storage,double);
296 EXTRACT("truncation",truncation,int);
297 EXTRACT("restart",restart,int);
298 EXTRACT("sweeps",sweeps,int);
299 #undef EXTRACT
300 #undef EXTRACT_OPTION
301 }
302
303
304 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26