finley/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 "SystemMatrix.h"
#include "Options.h"
}
#include "escript/Data/Data.h"
#include "escript/Data/UtilC.h"
#include "SystemMatrixAdapter.h" 
#include "FinleyAdapterException.h" 
#include "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;
       case  ESCRIPT_AMG:
          return PASO_AMG;
       case  ESCRIPT_RILU:
          return PASO_RILU;
       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];
    strcpy(fName,fileName.c_str());
    Paso_SystemMatrix* mat = getPaso_SystemMatrix();
    Paso_SystemMatrix_saveMM(mat,fName);
    checkPasoError();
}

void SystemMatrixAdapter::saveHB(const std::string& fileName) const
{
    char fName[fileName.size()+1];
    strcpy(fName,fileName.c_str());
    Paso_SystemMatrix* mat = getPaso_SystemMatrix();
    Paso_SystemMatrix_saveHB(mat,fName);
    checkPasoError();
}

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

}  // end of namespace
1	/*
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	$Id$
14	*/
15	extern "C" {
16	#include "SystemMatrix.h"
17	#include "Options.h"
18	}
19	#include "escript/Data/Data.h"
20	#include "escript/Data/UtilC.h"
21	#include "SystemMatrixAdapter.h"
22	#include "FinleyAdapterException.h"
23	#include "FinleyError.h"
24	#include <boost/python/extract.hpp>
25
26	using namespace std;
27
28	namespace finley {
29
30	struct null_deleter
31	{
32	void operator()(void const *ptr) const
33	{
34	}
35	};
36
37
38	SystemMatrixAdapter::SystemMatrixAdapter()
39	{
40	throw FinleyAdapterException("Error - Illegal to generate default SystemMatrixAdapter.");
41	}
42
43	SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
44	const int row_blocksize,
45	const escript::FunctionSpace& row_functionspace,
46	const int column_blocksize,
47	const escript::FunctionSpace& column_functionspace):
48	AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
49	{
50	m_system_matrix.reset(system_matrix,null_deleter());
51	}
52
53	SystemMatrixAdapter::~SystemMatrixAdapter()
54	{
55	if (m_system_matrix.unique()) {
56	Paso_SystemMatrix* mat=m_system_matrix.get();
57	Paso_SystemMatrix_dealloc(mat);
58	}
59	}
60
61	Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const
62	{
63	return m_system_matrix.get();
64	}
65
66	void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const
67	{
68	Paso_SystemMatrix* mat=getPaso_SystemMatrix();
69
70	if ( x.getDataPointSize() != getColumnBlockSize()) {
71	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	x.expand();
80	y.expand();
81	double* x_dp=x.getSampleData(0);
82	double* y_dp=y.getSampleData(0);
83	Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
84	checkPasoError();
85	}
86
87	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	case ESCRIPT_AMG:
134	return PASO_AMG;
135	case ESCRIPT_RILU:
136	return PASO_RILU;
137	default:
138	stringstream temp;
139	temp << "Error - Cannot map option value "<< option << " onto Paso";
140	throw FinleyAdapterException(temp.str());
141	}
142	}
143
144	void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
145	{
146	Paso_SystemMatrix* mat=getPaso_SystemMatrix();
147	Paso_Options paso_options;
148	Paso_Options_setDefaults(&paso_options);
149	// extract options
150	#define EXTRACT(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=boost::python::extract<__type__>(options.get(__key__))
151	#define EXTRACT_OPTION(__key__,__val__,__type__) if ( options.has_key(__key__)) paso_options.__val__=mapOptionToPaso(boost::python::extract<__type__>(options.get(__key__)))
152	EXTRACT("verbose",verbose,int);
153	EXTRACT_OPTION("reordering",reordering,int);
154	EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
155	EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
156	EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
157	EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
158	EXTRACT_OPTION("preconditioner",preconditioner,int);
159	EXTRACT("iter_max",iter_max,int);
160	EXTRACT("drop_tolerance",drop_tolerance,double);
161	EXTRACT("drop_storage",drop_storage,double);
162	EXTRACT("truncation",truncation,double);
163	EXTRACT("restart",restart,double);
164	#undef EXTRACT
165	#undef EXTRACT_OPTION
166	if ( out.getDataPointSize() != getColumnBlockSize()) {
167	throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
168	} else if ( in.getDataPointSize() != getRowBlockSize()) {
169	throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
170	} else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
171	throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
172	} else if (in.getFunctionSpace() != getRowFunctionSpace()) {
173	throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
174	}
175	out.expand();
176	in.expand();
177	double* out_dp=out.getSampleData(0);
178	double* in_dp=in.getSampleData(0);
179	Paso_solve(mat,out_dp,in_dp,&paso_options);
180	checkPasoError();
181	}
182
183	void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
184	{
185	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
186	if ( col_q.getDataPointSize() != getColumnBlockSize()) {
187	throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
188	} else if ( row_q.getDataPointSize() != getRowBlockSize()) {
189	throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
190	} else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
191	throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
192	} else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
193	throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
194	}
195	row_q.expand();
196	col_q.expand();
197	double* row_q_dp=row_q.getSampleData(0);
198	double* col_q_dp=col_q.getSampleData(0);
199	Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
200	checkPasoError();
201	}
202
203	void SystemMatrixAdapter::saveMM(const std::string& fileName) const
204	{
205	char fName[fileName.size()+1];
206	strcpy(fName,fileName.c_str());
207	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
208	Paso_SystemMatrix_saveMM(mat,fName);
209	checkPasoError();
210	}
211
212	void SystemMatrixAdapter::saveHB(const std::string& fileName) const
213	{
214	char fName[fileName.size()+1];
215	strcpy(fName,fileName.c_str());
216	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
217	Paso_SystemMatrix_saveHB(mat,fName);
218	checkPasoError();
219	}
220
221	void SystemMatrixAdapter::resetValues() const
222	{
223	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
224	Paso_SystemMatrix_setValues(mat,0.);
225	Paso_solve_free(mat);
226	checkPasoError();
227	}
228
229	} // end of namespace
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision