/[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 1819 - (hide annotations)
Tue Sep 30 05:58:06 2008 UTC (12 years, 11 months ago) by artak
File size: 12098 byte(s)
Firs version of symmetric Gauss-Seidel preconditioner with coloring
1 jgs 472
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
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 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     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 jgs 82 }
80    
81 jgs 150 int SystemMatrixAdapter::mapOptionToPaso(const int option) {
82     switch (option) {
83 gross 1639 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 artak 1703 case ESCRIPT_TFQMR:
136     return PASO_TFQMR;
137 artak 1787 case ESCRIPT_MINRES:
138     return PASO_MINRES;
139 artak 1819 case ESCRIPT_GS:
140     return PASO_GS;
141 gross 1639 default:
142     stringstream temp;
143     temp << "Error - Cannot map option value "<< option << " onto Paso";
144     throw FinleyAdapterException(temp.str());
145     }
146 jgs 150 }
147    
148 ksteube 1339 void finley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
149     {
150 phornby 1628 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 ksteube 1339
156 phornby 1628 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
157 ksteube 1339
158 phornby 1628 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 ksteube 1339
168 phornby 1628 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 ksteube 1339
186 phornby 1628 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 ksteube 1339
195     }
196    
197 jgs 153 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
198 jgs 82 {
199 phornby 1628 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 jgs 82 }
218    
219 jgs 153 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
220 jgs 82 {
221 phornby 1628 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 jgs 82 }
238    
239 jgs 102 void SystemMatrixAdapter::saveMM(const std::string& fileName) const
240     {
241 phornby 1628 if( fileName.size() == 0 )
242     {
243     throw FinleyAdapterException("Null file name!");
244     }
245    
246     char *fName = TMPMEMALLOC(fileName.size()+1,char);
247 woo409 757
248 phornby 1628 strcpy(fName,fileName.c_str());
249     Paso_SystemMatrix* mat = getPaso_SystemMatrix();
250     Paso_SystemMatrix_saveMM(mat,fName);
251     checkPasoError();
252     TMPMEMFREE(fName);
253 woo409 757
254 jgs 102 }
255    
256 jgs 123 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
257     {
258 phornby 1628 if( fileName.size() == 0 )
259     {
260     throw FinleyAdapterException("Null file name!");
261     }
262 woo409 757
263 phornby 1628 char *fName = TMPMEMALLOC(fileName.size()+1,char);
264 woo409 757
265 phornby 1628 strcpy(fName,fileName.c_str());
266     Paso_SystemMatrix* mat = getPaso_SystemMatrix();
267     Paso_SystemMatrix_saveHB(mat,fName);
268     checkPasoError();
269     TMPMEMFREE(fName);
270    
271 jgs 123 }
272    
273 jgs 149 void SystemMatrixAdapter::resetValues() const
274 jgs 102 {
275 jgs 150 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
276     Paso_SystemMatrix_setValues(mat,0.);
277     Paso_solve_free(mat);
278     checkPasoError();
279 jgs 108 }
280 jgs 102
281 gross 1364 void SystemMatrixAdapter::dictToPasoOptions(Paso_Options* paso_options, const boost::python::dict& options)
282     {
283 phornby 1628 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 artak 1819 EXTRACT("precNumSteps",precNumSteps,int);
299 phornby 1628 #undef EXTRACT
300     #undef EXTRACT_OPTION
301 gross 1364 }
302 phornby 1628
303 gross 1364
304 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