/[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 1787 - (hide annotations)
Mon Sep 15 01:36:34 2008 UTC (11 years, 6 months ago) by artak
File size: 12037 byte(s)
MINRES solver is added to escript. Additional 16 tests are added to run_simplesolve for MINRES and TFQMR solvers
1 jgs 472
2 ksteube 1312 /* $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 ksteube 817 #ifdef PASO_MPI
17     #include <mpi.h>
18     #endif
19 jgs 203 #include "SystemMatrixAdapter.h"
20 jgs 82
21     using namespace std;
22    
23     namespace finley {
24    
25 jgs 102 struct null_deleter
26     {
27 phornby 1628 void operator()(void const *ptr) const
28     {
29     }
30 jgs 102 };
31    
32    
33 jgs 82 SystemMatrixAdapter::SystemMatrixAdapter()
34     {
35     throw FinleyAdapterException("Error - Illegal to generate default SystemMatrixAdapter.");
36     }
37    
38 jgs 150 SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
39 jgs 82 const int row_blocksize,
40     const escript::FunctionSpace& row_functionspace,
41     const int column_blocksize,
42     const escript::FunctionSpace& column_functionspace):
43 jgs 102 AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
44 jgs 82 {
45 phornby 1628 m_system_matrix.reset(system_matrix,null_deleter());
46 jgs 82 }
47    
48     SystemMatrixAdapter::~SystemMatrixAdapter()
49     {
50 phornby 1628 if (m_system_matrix.unique()) {
51     Paso_SystemMatrix* mat=m_system_matrix.get();
52     Paso_SystemMatrix_free(mat);
53     }
54 jgs 82 }
55    
56 jgs 150 Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const
57 jgs 82 {
58     return m_system_matrix.get();
59     }
60    
61 jgs 153 void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const
62 jgs 82 {
63 jgs 150 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
64    
65 phornby 1628 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 jgs 82 }
81    
82 jgs 150 int SystemMatrixAdapter::mapOptionToPaso(const int option) {
83     switch (option) {
84 gross 1639 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     case ESCRIPT_NONLINEAR_GMRES:
135     return PASO_NONLINEAR_GMRES;
136 artak 1703 case ESCRIPT_TFQMR:
137     return PASO_TFQMR;
138 artak 1787 case ESCRIPT_MINRES:
139     return PASO_MINRES;
140 gross 1639 default:
141     stringstream temp;
142     temp << "Error - Cannot map option value "<< option << " onto Paso";
143     throw FinleyAdapterException(temp.str());
144     }
145 jgs 150 }
146    
147 ksteube 1339 void finley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
148     {
149 phornby 1628 Paso_SystemMatrix* mat=m_system_matrix.get();
150     int first_row_index = mat->row_distribution->first_component[mat->mpi_info->rank];
151     int last_row_index = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
152     int first_col_index = mat->col_distribution->first_component[mat->mpi_info->rank];
153     int last_col_index = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;
154 ksteube 1339
155 phornby 1628 fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
156 ksteube 1339
157 phornby 1628 switch (mat->type) {
158     case MATRIX_FORMAT_DEFAULT: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n"); break;
159     case MATRIX_FORMAT_CSC: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n"); break;
160     case MATRIX_FORMAT_SYM: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_SYM\n"); break;
161     case MATRIX_FORMAT_BLK1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n"); break;
162     case MATRIX_FORMAT_OFFSET1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n"); break;
163     case MATRIX_FORMAT_TRILINOS_CRS: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n"); break;
164     default: fprintf(stdout, "\tMatrix type unknown\n"); break;
165     }
166 ksteube 1339
167 phornby 1628 fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
168     fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
169     fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
170     fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
171     fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
172     fprintf(stdout, "\tcol_coupleBlock numRows %d\n", mat->col_coupleBlock->numRows);
173     fprintf(stdout, "\tcol_coupleBlock numCols %d\n", mat->col_coupleBlock->numCols);
174     fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", mat->col_coupleBlock->pattern->numOutput);
175     fprintf(stdout, "\trow_coupleBlock numRows %d\n", mat->row_coupleBlock->numRows);
176     fprintf(stdout, "\trow_coupleBlock numCols %d\n", mat->row_coupleBlock->numCols);
177     fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", mat->row_coupleBlock->pattern->numOutput);
178     fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
179     fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
180     fprintf(stdout, "\tblock_size %d\n", mat->block_size);
181     fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
182     fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
183     fprintf(stdout, "\tlogical_block_size %d\n", mat->logical_block_size);
184 ksteube 1339
185 phornby 1628 if (full) {
186     printf("\trow_distribution: ");
187     for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->row_distribution[i]);
188     printf("\n");
189     printf("\tcol_distribution: ");
190     for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->col_distribution[i]);
191     printf("\n");
192     }
193 ksteube 1339
194     }
195    
196 jgs 153 void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
197 jgs 82 {
198 phornby 1628 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
199     Paso_Options paso_options;
200     dictToPasoOptions(&paso_options,options);
201     if ( out.getDataPointSize() != getColumnBlockSize()) {
202     throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
203     } else if ( in.getDataPointSize() != getRowBlockSize()) {
204     throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
205     } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
206     throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
207     } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
208     throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
209     }
210     out.expand();
211     in.expand();
212     double* out_dp=out.getSampleData(0);
213     double* in_dp=in.getSampleData(0);
214     Paso_solve(mat,out_dp,in_dp,&paso_options);
215     checkPasoError();
216 jgs 82 }
217    
218 jgs 153 void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
219 jgs 82 {
220 phornby 1628 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
221     if ( col_q.getDataPointSize() != getColumnBlockSize()) {
222     throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
223     } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
224     throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
225     } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
226     throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
227     } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
228     throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
229     }
230     row_q.expand();
231     col_q.expand();
232     double* row_q_dp=row_q.getSampleData(0);
233     double* col_q_dp=col_q.getSampleData(0);
234     Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
235     checkPasoError();
236 jgs 82 }
237    
238 jgs 102 void SystemMatrixAdapter::saveMM(const std::string& fileName) const
239     {
240 phornby 1628 if( fileName.size() == 0 )
241     {
242     throw FinleyAdapterException("Null file name!");
243     }
244    
245     char *fName = TMPMEMALLOC(fileName.size()+1,char);
246 woo409 757
247 phornby 1628 strcpy(fName,fileName.c_str());
248     Paso_SystemMatrix* mat = getPaso_SystemMatrix();
249     Paso_SystemMatrix_saveMM(mat,fName);
250     checkPasoError();
251     TMPMEMFREE(fName);
252 woo409 757
253 jgs 102 }
254    
255 jgs 123 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
256     {
257 phornby 1628 if( fileName.size() == 0 )
258     {
259     throw FinleyAdapterException("Null file name!");
260     }
261 woo409 757
262 phornby 1628 char *fName = TMPMEMALLOC(fileName.size()+1,char);
263 woo409 757
264 phornby 1628 strcpy(fName,fileName.c_str());
265     Paso_SystemMatrix* mat = getPaso_SystemMatrix();
266     Paso_SystemMatrix_saveHB(mat,fName);
267     checkPasoError();
268     TMPMEMFREE(fName);
269    
270 jgs 123 }
271    
272 jgs 149 void SystemMatrixAdapter::resetValues() const
273 jgs 102 {
274 jgs 150 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
275     Paso_SystemMatrix_setValues(mat,0.);
276     Paso_solve_free(mat);
277     checkPasoError();
278 jgs 108 }
279 jgs 102
280 gross 1364 void SystemMatrixAdapter::dictToPasoOptions(Paso_Options* paso_options, const boost::python::dict& options)
281     {
282 phornby 1628 Paso_Options_setDefaults(paso_options);
283     #define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options->__val__=boost::python::extract<__type__>(options.get(__key__))
284     #define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options->__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
285     EXTRACT("verbose",verbose,int);
286     EXTRACT_OPTION("reordering",reordering,int);
287     EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
288     EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
289     EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
290     EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
291     EXTRACT_OPTION("preconditioner",preconditioner,int);
292     EXTRACT("iter_max",iter_max,int);
293     EXTRACT("drop_tolerance",drop_tolerance,double);
294     EXTRACT("drop_storage",drop_storage,double);
295     EXTRACT("truncation",truncation,int);
296     EXTRACT("restart",restart,int);
297     #undef EXTRACT
298     #undef EXTRACT_OPTION
299 gross 1364 }
300 phornby 1628
301 gross 1364
302 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