src/CPPAdapter/SystemMatrixAdapter.cpp

/*
 ******************************************************************************
 *                                                                            *
 *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *
 *                                                                            *
 * This software is the property of ACcESS. No part of this code              *
 * may be copied in any form or by any means without the expressed written    *
 * consent of ACcESS.  Copying, use or modification of this software          *
 * by any unauthorised person is illegal unless that person has a software    *
 * license agreement with ACcESS.                                             *
 *                                                                            *
 ******************************************************************************
 $Id$
*/
extern "C" {
#include "paso/SystemMatrix.h"
#include "paso/Options.h"
}
#include "escript/Data/Data.h"
#include "escript/Data/UtilC.h"
#include "finley/CPPAdapter/SystemMatrixAdapter.h" 
#include "finley/CPPAdapter/FinleyAdapterException.h" 
#include "finley/CPPAdapter/FinleyError.h"
#include <boost/python/extract.hpp>

using namespace std;

namespace finley {

struct null_deleter
{
  void operator()(void const *ptr) const
  {
  }
};


SystemMatrixAdapter::SystemMatrixAdapter()
{
   throw FinleyAdapterException("Error - Illegal to generate default SystemMatrixAdapter.");
}

SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
                                         const int row_blocksize,
                                         const escript::FunctionSpace& row_functionspace,
                                         const int column_blocksize,
                                         const escript::FunctionSpace& column_functionspace):
AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
{
    m_system_matrix.reset(system_matrix,null_deleter());
}

SystemMatrixAdapter::~SystemMatrixAdapter()
{ 
    if (m_system_matrix.unique()) {
        Paso_SystemMatrix* mat=m_system_matrix.get();
        Paso_SystemMatrix_dealloc(mat);
    }
}

Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const 
{
   return m_system_matrix.get();
}

void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const 
{
   Paso_SystemMatrix* mat=getPaso_SystemMatrix();

  if ( x.getDataPointSize()  != getColumnBlockSize()) {
   throw FinleyAdapterException("matrix vector product : column block size does not match the number of components in input.");
  } else if (y.getDataPointSize() != getRowBlockSize()) {
   throw FinleyAdapterException("matrix vector product : row block size does not match the number of components in output.");
  } else if ( x.getFunctionSpace()  != getColumnFunctionSpace()) {
   throw FinleyAdapterException("matrix vector product : column function space and function space of input don't match.");
  } else if (y.getFunctionSpace() != getRowFunctionSpace()) {
   throw FinleyAdapterException("matrix vector product : row function space and function space of output don't match.");
  }
  x.expand();
  y.expand();
  double* x_dp=x.getSampleData(0);
  double* y_dp=y.getSampleData(0);
  Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
  checkPasoError();
}

int SystemMatrixAdapter::mapOptionToPaso(const int option)  {
   switch (option) {
       case  ESCRIPT_DEFAULT:
          return PASO_DEFAULT;
       case  ESCRIPT_DIRECT:
          return PASO_DIRECT;
       case  ESCRIPT_CHOLEVSKY:
          return PASO_CHOLEVSKY;
       case  ESCRIPT_PCG:
          return PASO_PCG;
       case  ESCRIPT_CR:
          return PASO_CR;
       case  ESCRIPT_CGS:
          return PASO_CGS;
       case  ESCRIPT_BICGSTAB:
          return PASO_BICGSTAB;
       case  ESCRIPT_SSOR:
          return PASO_SSOR;
       case  ESCRIPT_ILU0:
          return PASO_ILU0;
       case  ESCRIPT_ILUT:
          return PASO_ILUT;
       case  ESCRIPT_JACOBI:
          return PASO_JACOBI;
       case  ESCRIPT_GMRES:
          return PASO_GMRES;
       case  ESCRIPT_PRES20:
          return PASO_PRES20;
       case  ESCRIPT_NO_REORDERING:
          return PASO_NO_REORDERING;
       case  ESCRIPT_MINIMUM_FILL_IN:
          return PASO_MINIMUM_FILL_IN;
       case  ESCRIPT_NESTED_DISSECTION:
          return PASO_NESTED_DISSECTION;
       case  ESCRIPT_SCSL:
          return PASO_SCSL;
       case  ESCRIPT_MKL:
          return PASO_MKL;
       case  ESCRIPT_UMFPACK:
          return PASO_UMFPACK;
       case  ESCRIPT_ITERATIVE:
          return PASO_ITERATIVE;
       case  ESCRIPT_PASO:
          return PASO_PASO;
       case  ESCRIPT_LUMPING:
          return PASO_LUMPING;
       default:
           stringstream temp;
           temp << "Error - Cannot map option value "<< option << " onto Paso";
           throw FinleyAdapterException(temp.str());
    }
}

void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
{
    Paso_SystemMatrix* mat=getPaso_SystemMatrix();
    Paso_Options paso_options;
    Paso_Options_setDefaults(&paso_options);
    // extract options 
    #define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=boost::python::extract<__type__>(options.get(__key__))
    #define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
    EXTRACT("verbose",verbose,int);
    EXTRACT_OPTION("reordering",reordering,int);
    EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
    EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
    EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
    EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
    EXTRACT_OPTION("preconditioner",preconditioner,int);
    EXTRACT("iter_max",iter_max,int);
    EXTRACT("drop_tolerance",drop_tolerance,double);
    EXTRACT("drop_storage",drop_storage,double);
    EXTRACT("truncation",truncation,double);
    EXTRACT("restart",restart,double);
    #undef EXTRACT
    #undef EXTRACT_OPTION
    if ( out.getDataPointSize()  != getColumnBlockSize()) {
     throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
    } else if ( in.getDataPointSize() != getRowBlockSize()) {
     throw FinleyAdapterException("solve : row block size does not match the number of components of  right hand side.");
    } else if ( out.getFunctionSpace()  != getColumnFunctionSpace()) {
     throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
    } else if (in.getFunctionSpace() != getRowFunctionSpace()) {
     throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
    }
    out.expand();
    in.expand();
    double* out_dp=out.getSampleData(0);
    double* in_dp=in.getSampleData(0);
    Paso_solve(mat,out_dp,in_dp,&paso_options);
    checkPasoError();
}

void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
{
    Paso_SystemMatrix* mat = getPaso_SystemMatrix();
    if ( col_q.getDataPointSize()  != getColumnBlockSize()) {
     throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
    } else if ( row_q.getDataPointSize() != getRowBlockSize()) {
     throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
    } else if ( col_q.getFunctionSpace()  != getColumnFunctionSpace()) {
     throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
    } else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
     throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
    }
    row_q.expand();
    col_q.expand();
    double* row_q_dp=row_q.getSampleData(0);
    double* col_q_dp=col_q.getSampleData(0);
    Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
    checkPasoError();
}

void SystemMatrixAdapter::saveMM(const std::string& fileName) const
{
    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
    
    strcpy(fName,fileName.c_str());
    Paso_SystemMatrix* mat = getPaso_SystemMatrix();
    Paso_SystemMatrix_saveMM(mat,fName);
    checkPasoError();
    TMPMEMFREE(fName);

}

void SystemMatrixAdapter::saveHB(const std::string& fileName) const
{
    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;

    strcpy(fName,fileName.c_str());
    Paso_SystemMatrix* mat = getPaso_SystemMatrix();
    Paso_SystemMatrix_saveHB(mat,fName);
    checkPasoError();
    TMPMEMFREE(fName);

}

void SystemMatrixAdapter::resetValues() const
{
   Paso_SystemMatrix* mat = getPaso_SystemMatrix();
   Paso_SystemMatrix_setValues(mat,0.);
   Paso_solve_free(mat);
   checkPasoError();
}

}  // end of namespace
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	jgs	153	void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const
67	jgs	82	{
68	jgs	150	Paso_SystemMatrix* mat=getPaso_SystemMatrix();
69
70	jgs	153	if ( x.getDataPointSize() != getColumnBlockSize()) {
71	jgs	150	throw FinleyAdapterException("matrix vector product : column block size does not match the number of components in input.");
72			} else if (y.getDataPointSize() != getRowBlockSize()) {
73			throw FinleyAdapterException("matrix vector product : row block size does not match the number of components in output.");
74			} else if ( x.getFunctionSpace() != getColumnFunctionSpace()) {
75			throw FinleyAdapterException("matrix vector product : column function space and function space of input don't match.");
76			} else if (y.getFunctionSpace() != getRowFunctionSpace()) {
77			throw FinleyAdapterException("matrix vector product : row function space and function space of output don't match.");
78			}
79	jgs	153	x.expand();
80			y.expand();
81			double* x_dp=x.getSampleData(0);
82	jgs	150	double* y_dp=y.getSampleData(0);
83			Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
84			checkPasoError();
85	jgs	82	}
86
87	jgs	150	int SystemMatrixAdapter::mapOptionToPaso(const int option) {
88			switch (option) {
89			case ESCRIPT_DEFAULT:
90			return PASO_DEFAULT;
91			case ESCRIPT_DIRECT:
92			return PASO_DIRECT;
93			case ESCRIPT_CHOLEVSKY:
94			return PASO_CHOLEVSKY;
95			case ESCRIPT_PCG:
96			return PASO_PCG;
97			case ESCRIPT_CR:
98			return PASO_CR;
99			case ESCRIPT_CGS:
100			return PASO_CGS;
101			case ESCRIPT_BICGSTAB:
102			return PASO_BICGSTAB;
103			case ESCRIPT_SSOR:
104			return PASO_SSOR;
105			case ESCRIPT_ILU0:
106			return PASO_ILU0;
107			case ESCRIPT_ILUT:
108			return PASO_ILUT;
109			case ESCRIPT_JACOBI:
110			return PASO_JACOBI;
111			case ESCRIPT_GMRES:
112			return PASO_GMRES;
113			case ESCRIPT_PRES20:
114			return PASO_PRES20;
115			case ESCRIPT_NO_REORDERING:
116			return PASO_NO_REORDERING;
117			case ESCRIPT_MINIMUM_FILL_IN:
118			return PASO_MINIMUM_FILL_IN;
119			case ESCRIPT_NESTED_DISSECTION:
120			return PASO_NESTED_DISSECTION;
121			case ESCRIPT_SCSL:
122			return PASO_SCSL;
123			case ESCRIPT_MKL:
124			return PASO_MKL;
125			case ESCRIPT_UMFPACK:
126			return PASO_UMFPACK;
127			case ESCRIPT_ITERATIVE:
128			return PASO_ITERATIVE;
129			case ESCRIPT_PASO:
130			return PASO_PASO;
131			case ESCRIPT_LUMPING:
132			return PASO_LUMPING;
133			default:
134			stringstream temp;
135			temp << "Error - Cannot map option value "<< option << " onto Paso";
136			throw FinleyAdapterException(temp.str());
137			}
138			}
139
140	jgs	153	void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
141	jgs	82	{
142	jgs	150	Paso_SystemMatrix* mat=getPaso_SystemMatrix();
143			Paso_Options paso_options;
144			Paso_Options_setDefaults(&paso_options);
145	jgs	82	// extract options
146	jgs	150	#define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=boost::python::extract<__type__>(options.get(__key__))
147			#define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
148	jgs	82	EXTRACT("verbose",verbose,int);
149	jgs	150	EXTRACT_OPTION("reordering",reordering,int);
150	jgs	102	EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
151	jgs	150	EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
152	jgs	102	EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
153	jgs	150	EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
154			EXTRACT_OPTION("preconditioner",preconditioner,int);
155	jgs	82	EXTRACT("iter_max",iter_max,int);
156			EXTRACT("drop_tolerance",drop_tolerance,double);
157			EXTRACT("drop_storage",drop_storage,double);
158	jgs	102	EXTRACT("truncation",truncation,double);
159			EXTRACT("restart",restart,double);
160	jgs	82	#undef EXTRACT
161	jgs	150	#undef EXTRACT_OPTION
162	jgs	153	if ( out.getDataPointSize() != getColumnBlockSize()) {
163	jgs	150	throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
164			} else if ( in.getDataPointSize() != getRowBlockSize()) {
165			throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
166			} else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
167			throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
168			} else if (in.getFunctionSpace() != getRowFunctionSpace()) {
169			throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
170			}
171	jgs	153	out.expand();
172			in.expand();
173	jgs	150	double* out_dp=out.getSampleData(0);
174	jgs	153	double* in_dp=in.getSampleData(0);
175	jgs	150	Paso_solve(mat,out_dp,in_dp,&paso_options);
176			checkPasoError();
177	jgs	82	}
178
179	jgs	153	void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
180	jgs	82	{
181	jgs	150	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
182	jgs	153	if ( col_q.getDataPointSize() != getColumnBlockSize()) {
183	jgs	150	throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
184			} else if ( row_q.getDataPointSize() != getRowBlockSize()) {
185			throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
186			} else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
187			throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
188			} else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
189			throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
190			}
191	jgs	153	row_q.expand();
192			col_q.expand();
193			double* row_q_dp=row_q.getSampleData(0);
194			double* col_q_dp=col_q.getSampleData(0);
195	jgs	150	Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
196			checkPasoError();
197	jgs	82	}
198
199	jgs	102	void SystemMatrixAdapter::saveMM(const std::string& fileName) const
200			{
201	robwdcock	186	char fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char)NULL;
202
203			strcpy(fName,fileName.c_str());
204	jgs	150	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
205			Paso_SystemMatrix_saveMM(mat,fName);
206			checkPasoError();
207	robwdcock	186	TMPMEMFREE(fName);
208
209	jgs	102	}
210
211	jgs	123	void SystemMatrixAdapter::saveHB(const std::string& fileName) const
212			{
213	robwdcock	186	char fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char)NULL;
214
215	jgs	123	strcpy(fName,fileName.c_str());
216	jgs	150	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
217			Paso_SystemMatrix_saveHB(mat,fName);
218			checkPasoError();
219	robwdcock	186	TMPMEMFREE(fName);
220
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
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision