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

Annotation of /branches/RW_WIN32/finley/src/CPPAdapter/SystemMatrixAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (15 years, 5 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/SystemMatrixAdapter.cpp
File size: 9765 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 jgs 82 /*
2     ******************************************************************************
3     * *
4     * COPYRIGHT ACcESS 2004 - All Rights Reserved *
5     * *
6     * This software is the property of ACcESS. No part of this code *
7     * may be copied in any form or by any means without the expressed written *
8     * consent of ACcESS. Copying, use or modification of this software *
9     * by any unauthorised person is illegal unless that person has a software *
10     * license agreement with ACcESS. *
11     * *
12     ******************************************************************************
13 jgs 102 $Id$
14 jgs 82 */
15     extern "C" {
16 jgs 150 #include "paso/SystemMatrix.h"
17     #include "paso/Options.h"
18 jgs 82 }
19     #include "escript/Data/Data.h"
20 jgs 150 #include "escript/Data/UtilC.h"
21 jgs 82 #include "finley/CPPAdapter/SystemMatrixAdapter.h"
22     #include "finley/CPPAdapter/FinleyAdapterException.h"
23     #include "finley/CPPAdapter/FinleyError.h"
24     #include <boost/python/extract.hpp>
25    
26     using namespace std;
27    
28     namespace finley {
29    
30 jgs 102 struct null_deleter
31     {
32     void operator()(void const *ptr) const
33     {
34     }
35     };
36    
37    
38 jgs 82 SystemMatrixAdapter::SystemMatrixAdapter()
39     {
40     throw FinleyAdapterException("Error - Illegal to generate default SystemMatrixAdapter.");
41     }
42    
43 jgs 150 SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
44 jgs 82 const int row_blocksize,
45     const escript::FunctionSpace& row_functionspace,
46     const int column_blocksize,
47     const escript::FunctionSpace& column_functionspace):
48 jgs 102 AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
49 jgs 82 {
50 jgs 102 m_system_matrix.reset(system_matrix,null_deleter());
51 jgs 82 }
52    
53     SystemMatrixAdapter::~SystemMatrixAdapter()
54     {
55     if (m_system_matrix.unique()) {
56 jgs 150 Paso_SystemMatrix* mat=m_system_matrix.get();
57     Paso_SystemMatrix_dealloc(mat);
58 jgs 82 }
59     }
60    
61 jgs 150 Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const
62 jgs 82 {
63     return m_system_matrix.get();
64     }
65    
66     void SystemMatrixAdapter::ypAx(escript::Data& y, const escript::Data& x) const
67     {
68 jgs 150 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
69    
70     if (!x.isExpanded()) {
71     throw FinleyAdapterException("matrix vector product : input Data object has to be expanded");
72     } else if (!y.isExpanded()) {
73     throw FinleyAdapterException("matrix vector product : output Data object has to be expanded");
74     } else if ( x.getDataPointSize() != getColumnBlockSize()) {
75     throw FinleyAdapterException("matrix vector product : column block size does not match the number of components in input.");
76     } else if (y.getDataPointSize() != getRowBlockSize()) {
77     throw FinleyAdapterException("matrix vector product : row block size does not match the number of components in output.");
78     } else if ( x.getFunctionSpace() != getColumnFunctionSpace()) {
79     throw FinleyAdapterException("matrix vector product : column function space and function space of input don't match.");
80     } else if (y.getFunctionSpace() != getRowFunctionSpace()) {
81     throw FinleyAdapterException("matrix vector product : row function space and function space of output don't match.");
82     }
83     double* x_dp=((escript::Data)x).getSampleData(0);
84     double* y_dp=y.getSampleData(0);
85     Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
86     checkPasoError();
87 jgs 82 }
88    
89 jgs 150 int SystemMatrixAdapter::mapOptionToPaso(const int option) {
90     switch (option) {
91     case ESCRIPT_DEFAULT:
92     return PASO_DEFAULT;
93     case ESCRIPT_DIRECT:
94     return PASO_DIRECT;
95     case ESCRIPT_CHOLEVSKY:
96     return PASO_CHOLEVSKY;
97     case ESCRIPT_PCG:
98     return PASO_PCG;
99     case ESCRIPT_CR:
100     return PASO_CR;
101     case ESCRIPT_CGS:
102     return PASO_CGS;
103     case ESCRIPT_BICGSTAB:
104     return PASO_BICGSTAB;
105     case ESCRIPT_SSOR:
106     return PASO_SSOR;
107     case ESCRIPT_ILU0:
108     return PASO_ILU0;
109     case ESCRIPT_ILUT:
110     return PASO_ILUT;
111     case ESCRIPT_JACOBI:
112     return PASO_JACOBI;
113     case ESCRIPT_GMRES:
114     return PASO_GMRES;
115     case ESCRIPT_PRES20:
116     return PASO_PRES20;
117     case ESCRIPT_NO_REORDERING:
118     return PASO_NO_REORDERING;
119     case ESCRIPT_MINIMUM_FILL_IN:
120     return PASO_MINIMUM_FILL_IN;
121     case ESCRIPT_NESTED_DISSECTION:
122     return PASO_NESTED_DISSECTION;
123     case ESCRIPT_SCSL:
124     return PASO_SCSL;
125     case ESCRIPT_MKL:
126     return PASO_MKL;
127     case ESCRIPT_UMFPACK:
128     return PASO_UMFPACK;
129     case ESCRIPT_ITERATIVE:
130     return PASO_ITERATIVE;
131     case ESCRIPT_PASO:
132     return PASO_PASO;
133     case ESCRIPT_LUMPING:
134     return PASO_LUMPING;
135     default:
136     stringstream temp;
137     temp << "Error - Cannot map option value "<< option << " onto Paso";
138     throw FinleyAdapterException(temp.str());
139     }
140     }
141    
142 jgs 82 void SystemMatrixAdapter::setToSolution(escript::Data& out, const escript::Data& in, const boost::python::dict& options) const
143     {
144 jgs 150 Paso_SystemMatrix* mat=getPaso_SystemMatrix();
145     Paso_Options paso_options;
146     Paso_Options_setDefaults(&paso_options);
147 jgs 82 // extract options
148 jgs 150 #define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=boost::python::extract<__type__>(options.get(__key__))
149     #define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
150 jgs 82 EXTRACT("verbose",verbose,int);
151 jgs 150 EXTRACT_OPTION("reordering",reordering,int);
152 jgs 102 EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
153 jgs 150 EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
154 jgs 102 EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
155 jgs 150 EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
156     EXTRACT_OPTION("preconditioner",preconditioner,int);
157 jgs 82 EXTRACT("iter_max",iter_max,int);
158     EXTRACT("drop_tolerance",drop_tolerance,double);
159     EXTRACT("drop_storage",drop_storage,double);
160 jgs 102 EXTRACT("truncation",truncation,double);
161     EXTRACT("restart",restart,double);
162 jgs 82 #undef EXTRACT
163 jgs 150 #undef EXTRACT_OPTION
164     if (! out.isExpanded()) {
165     throw FinleyAdapterException("solve : result Data object has to be expanded");
166     } else if (!in.isExpanded()) {
167     throw FinleyAdapterException("solve : right hand side Data object has to be expanded");
168     } else if ( out.getDataPointSize() != getColumnBlockSize()) {
169     throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
170     } else if ( in.getDataPointSize() != getRowBlockSize()) {
171     throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
172     } else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
173     throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
174     } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
175     throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
176     }
177     double* out_dp=out.getSampleData(0);
178     double* in_dp=((escript::Data) in).getSampleData(0);
179     Paso_solve(mat,out_dp,in_dp,&paso_options);
180     checkPasoError();
181 jgs 82 }
182    
183     void SystemMatrixAdapter::nullifyRowsAndCols(const escript::Data& row_q, const escript::Data& col_q, const double mdv) const
184     {
185 jgs 150 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
186     if (! col_q.isExpanded()) {
187     throw FinleyAdapterException("nullifyRowsAndCols : column mask has to be expanded");
188     } else if (!row_q.isExpanded()) {
189     throw FinleyAdapterException("nullifyRowsAndCols : row mask has to be expanded");
190     } else if ( col_q.getDataPointSize() != getColumnBlockSize()) {
191     throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
192     } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
193     throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
194     } else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
195     throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
196     } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
197     throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
198     }
199     double* row_q_dp=((escript::Data) row_q).getSampleData(0);
200     double* col_q_dp=((escript::Data) col_q).getSampleData(0);
201     Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
202     checkPasoError();
203 jgs 82 }
204    
205 jgs 102 void SystemMatrixAdapter::saveMM(const std::string& fileName) const
206     {
207     char fName[fileName.size()+1];
208     strcpy(fName,fileName.c_str());
209 jgs 150 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
210     Paso_SystemMatrix_saveMM(mat,fName);
211     checkPasoError();
212 jgs 102 }
213    
214 jgs 123 void SystemMatrixAdapter::saveHB(const std::string& fileName) const
215     {
216     char fName[fileName.size()+1];
217     strcpy(fName,fileName.c_str());
218 jgs 150 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
219     Paso_SystemMatrix_saveHB(mat,fName);
220     checkPasoError();
221 jgs 123 }
222    
223 jgs 149 void SystemMatrixAdapter::resetValues() const
224 jgs 102 {
225 jgs 150 Paso_SystemMatrix* mat = getPaso_SystemMatrix();
226     Paso_SystemMatrix_setValues(mat,0.);
227     Paso_solve_free(mat);
228     checkPasoError();
229 jgs 108 }
230 jgs 102
231 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