src/CPPAdapter/SystemMatrixAdapter.cpp


/* $Id$ */

/*******************************************************
 *
 *           Copyright 2003-2007 by ACceSS MNRF
 *       Copyright 2007 by University of Queensland
 *
 *                http://esscc.uq.edu.au
 *        Primary Business: Queensland, Australia
 *  Licensed under the Open Software License version 3.0
 *     http://www.opensource.org/licenses/osl-3.0.php
 *
 *******************************************************/

#ifdef PASO_MPI
#include <mpi.h>
#endif
#include "SystemMatrixAdapter.h" 

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_free(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;
       case  ESCRIPT_TRILINOS:
          return PASO_TRILINOS;
       case  ESCRIPT_NONLINEAR_GMRES:
          return PASO_NONLINEAR_GMRES;
       case  ESCRIPT_TFQMR:
          return PASO_TFQMR;
       default:
           stringstream temp;
           temp << "Error - Cannot map option value "<< option << " onto Paso";
           throw FinleyAdapterException(temp.str());
    }
}

void finley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
{
   Paso_SystemMatrix* mat=m_system_matrix.get();
   int first_row_index  = mat->row_distribution->first_component[mat->mpi_info->rank];
   int last_row_index   = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
   int first_col_index  = mat->col_distribution->first_component[mat->mpi_info->rank];
   int last_col_index   = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;

   fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);

   switch (mat->type) {
   case MATRIX_FORMAT_DEFAULT:      fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n"); break;
   case MATRIX_FORMAT_CSC:      fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n"); break;
   case MATRIX_FORMAT_SYM:      fprintf(stdout, "\tMatrix type MATRIX_FORMAT_SYM\n"); break;
   case MATRIX_FORMAT_BLK1:     fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n"); break;
   case MATRIX_FORMAT_OFFSET1:      fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n"); break;
   case MATRIX_FORMAT_TRILINOS_CRS: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n"); break;
   default:             fprintf(stdout, "\tMatrix type unknown\n"); break;
   }

   fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
   fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
   fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
   fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
   fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
   fprintf(stdout, "\tcol_coupleBlock numRows %d\n", mat->col_coupleBlock->numRows);
   fprintf(stdout, "\tcol_coupleBlock numCols %d\n", mat->col_coupleBlock->numCols);
   fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", mat->col_coupleBlock->pattern->numOutput);
   fprintf(stdout, "\trow_coupleBlock numRows %d\n", mat->row_coupleBlock->numRows);
   fprintf(stdout, "\trow_coupleBlock numCols %d\n", mat->row_coupleBlock->numCols);
   fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", mat->row_coupleBlock->pattern->numOutput);
   fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
   fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
   fprintf(stdout, "\tblock_size %d\n", mat->block_size);
   fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
   fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
   fprintf(stdout, "\tlogical_block_size %d\n", mat->logical_block_size);

   if (full) {
      printf("\trow_distribution: ");
      for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->row_distribution[i]);
      printf("\n");
      printf("\tcol_distribution: ");
      for(int i=0; i<=mat->mpi_info->size; i++) printf("%3d ", mat->col_distribution[i]);
      printf("\n");
   }

}

void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, const boost::python::dict& options) const
{
   Paso_SystemMatrix* mat=getPaso_SystemMatrix();
   Paso_Options paso_options;
   dictToPasoOptions(&paso_options,options);
   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
{
   if( fileName.size() == 0 )
   {
      throw FinleyAdapterException("Null file name!");
   }

   char *fName = TMPMEMALLOC(fileName.size()+1,char);
    
   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
{
   if( fileName.size() == 0 )
   {
      throw FinleyAdapterException("Null file name!");
   }

   char *fName = TMPMEMALLOC(fileName.size()+1,char);

   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();
}

void SystemMatrixAdapter::dictToPasoOptions(Paso_Options* paso_options, const boost::python::dict& options) 
{
   Paso_Options_setDefaults(paso_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,int);
   EXTRACT("restart",restart,int);
#undef EXTRACT
#undef EXTRACT_OPTION
}
 

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