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, const escript::Data& x) const 
{
   Paso_SystemMatrix* mat=getPaso_SystemMatrix();

  if (!x.isExpanded()) {
   throw FinleyAdapterException("matrix vector product : input Data object has to be expanded");
  } else if (!y.isExpanded()) {
   throw FinleyAdapterException("matrix vector product : output Data object has to be expanded");
  } else 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.");
  }
  double* x_dp=((escript::Data)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, const 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.isExpanded()) {
     throw FinleyAdapterException("solve : result Data object has to be expanded");
    } else if (!in.isExpanded()) {
     throw FinleyAdapterException("solve : right hand side Data object has to be expanded");
    } else 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.");
    }
    double* out_dp=out.getSampleData(0);
    double* in_dp=((escript::Data) in).getSampleData(0);
    Paso_solve(mat,out_dp,in_dp,&paso_options);
    checkPasoError();
}

void SystemMatrixAdapter::nullifyRowsAndCols(const escript::Data& row_q, const escript::Data& col_q, const double mdv) const
{
    Paso_SystemMatrix* mat = getPaso_SystemMatrix();
    if (! col_q.isExpanded()) {
     throw FinleyAdapterException("nullifyRowsAndCols : column mask has to be expanded");
    } else if (!row_q.isExpanded()) {
     throw FinleyAdapterException("nullifyRowsAndCols : row mask has to be expanded");
    } else 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.");
    }
    double* row_q_dp=((escript::Data) row_q).getSampleData(0);
    double* col_q_dp=((escript::Data) 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 "paso/SystemMatrix.h"
17	#include "paso/Options.h"
18	}
19	#include "escript/Data/Data.h"
20	#include "escript/Data/UtilC.h"
21	#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	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, const escript::Data& x) const
67	{
68	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	}
88
89	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	void SystemMatrixAdapter::setToSolution(escript::Data& out, const escript::Data& in, const boost::python::dict& options) const
143	{
144	Paso_SystemMatrix* mat=getPaso_SystemMatrix();
145	Paso_Options paso_options;
146	Paso_Options_setDefaults(&paso_options);
147	// extract options
148	#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	EXTRACT("verbose",verbose,int);
151	EXTRACT_OPTION("reordering",reordering,int);
152	EXTRACT(ESCRIPT_TOLERANCE_KEY,tolerance,double);
153	EXTRACT_OPTION(ESCRIPT_METHOD_KEY,method,int);
154	EXTRACT(ESCRIPT_SYMMETRY_KEY,symmetric,int);
155	EXTRACT_OPTION(ESCRIPT_PACKAGE_KEY,package,int);
156	EXTRACT_OPTION("preconditioner",preconditioner,int);
157	EXTRACT("iter_max",iter_max,int);
158	EXTRACT("drop_tolerance",drop_tolerance,double);
159	EXTRACT("drop_storage",drop_storage,double);
160	EXTRACT("truncation",truncation,double);
161	EXTRACT("restart",restart,double);
162	#undef EXTRACT
163	#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	}
182
183	void SystemMatrixAdapter::nullifyRowsAndCols(const escript::Data& row_q, const escript::Data& col_q, const double mdv) const
184	{
185	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	}
204
205	void SystemMatrixAdapter::saveMM(const std::string& fileName) const
206	{
207	char fName[fileName.size()+1];
208	strcpy(fName,fileName.c_str());
209	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
210	Paso_SystemMatrix_saveMM(mat,fName);
211	checkPasoError();
212	}
213
214	void SystemMatrixAdapter::saveHB(const std::string& fileName) const
215	{
216	char fName[fileName.size()+1];
217	strcpy(fName,fileName.c_str());
218	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
219	Paso_SystemMatrix_saveHB(mat,fName);
220	checkPasoError();
221	}
222
223	void SystemMatrixAdapter::resetValues() const
224	{
225	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
226	Paso_SystemMatrix_setValues(mat,0.);
227	Paso_solve_free(mat);
228	checkPasoError();
229	}
230
231	} // end of namespace
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision