src/CPPAdapter/SystemMatrixAdapter.cpp


/*******************************************************
*
* Copyright (c) 2003-2009 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* 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();
   x.requireWrite();
   y.requireWrite();
   double* x_dp=x.getSampleDataRW(0);
   double* y_dp=y.getSampleDataRW(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_LUMPING:
          return PASO_LUMPING;
       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_MKL:
          return PASO_MKL;
       case  ESCRIPT_UMFPACK:
          return PASO_UMFPACK;
       case  ESCRIPT_ITERATIVE:
          return PASO_ITERATIVE;
       case  ESCRIPT_PASO:
          return PASO_PASO;
       case  ESCRIPT_AMG:
          return PASO_AMG;
       case  ESCRIPT_REC_ILU:
          return PASO_REC_ILU;
       case  ESCRIPT_TRILINOS:
          return PASO_TRILINOS;
       case  ESCRIPT_NONLINEAR_GMRES:
          return PASO_NONLINEAR_GMRES;
       case  ESCRIPT_TFQMR :
          return PASO_TFQMR;
       case  ESCRIPT_MINRES:
          return PASO_MINRES;
       case  ESCRIPT_GAUSS_SEIDEL:
          return PASO_GAUSS_SEIDEL;
       case  ESCRIPT_RILU:
          return PASO_RILU;
       case  ESCRIPT_DEFAULT_REORDERING:
          return PASO_DEFAULT_REORDERING;
       case  ESCRIPT_SUPER_LU:
          return PASO_SUPER_LU;
       case  ESCRIPT_PASTIX:
          return PASO_PASTIX;
       case  ESCRIPT_YAIR_SHAPIRA_COARSENING:
          return PASO_YAIR_SHAPIRA_COARSENING;
       case  ESCRIPT_RUGE_STUEBEN_COARSENING:
          return PASO_RUGE_STUEBEN_COARSENING;
       case  ESCRIPT_AGGREGATION_COARSENING:
          return PASO_AGGREGATION_COARSENING;
       case  ESCRIPT_NO_PRECONDITIONER:
          return PASO_NO_PRECONDITIONER;
       case  ESCRIPT_MIN_COARSE_MATRIX_SIZE:
          return PASO_MIN_COARSE_MATRIX_SIZE;
       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);

}

void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, boost::python::object& options) const
{
   Paso_SystemMatrix* mat=getPaso_SystemMatrix();
   Paso_Options paso_options;
   options.attr("resetDiagnostics")();
   escriptToPasoOptions(&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.getSampleDataRW(0);   
   double* in_dp=in.getSampleDataRW(0);     
   Paso_solve(mat,out_dp,in_dp,&paso_options);
   pasoToEscriptOptions(&paso_options,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();
   row_q.requireWrite();
   col_q.requireWrite();
   double* row_q_dp=row_q.getSampleDataRW(0);
   double* col_q_dp=col_q.getSampleDataRW(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::pasoToEscriptOptions(const Paso_Options* paso_options,boost::python::object& options) 
{
#define SET(__key__,__val__,__type__) options.attr("_updateDiagnostics")(__key__,(__type__)paso_options->__val__)
   SET("num_iter", num_iter, int);
   SET("num_level", num_level, int);
   SET("num_inner_iter", num_inner_iter, int);
   SET("time", time, double);
   SET("set_up_time", set_up_time, double);
   SET("net_time", net_time, double);
   SET("residual_norm", residual_norm, double);
   SET("converged",converged, bool);
#undef SET
}
void SystemMatrixAdapter::escriptToPasoOptions(Paso_Options* paso_options, const boost::python::object& options) 
{
#define EXTRACT(__key__,__val__,__type__) paso_options->__val__=boost::python::extract<__type__>(options.attr(__key__)())
#define EXTRACT_OPTION(__key__,__val__,__type__) paso_options->__val__=mapOptionToPaso(boost::python::extract<__type__>(options.attr(__key__)()))
 
   Paso_Options_setDefaults(paso_options);
   EXTRACT_OPTION("getSolverMethod",method,index_t);
   EXTRACT_OPTION("getPackage",package,index_t);
   EXTRACT("isVerbose",verbose,bool_t);
   EXTRACT("isSymmetric",symmetric,bool_t);
   EXTRACT("getTolerance",tolerance, double);
   EXTRACT("getAbsoluteTolerance",absolute_tolerance, double);
   EXTRACT("getInnerTolerance",inner_tolerance, double);
   EXTRACT("adaptInnerTolerance",adapt_inner_tolerance, bool_t);
   EXTRACT_OPTION("getReordering", reordering, index_t);
   EXTRACT_OPTION("getPreconditioner", preconditioner, index_t);
   EXTRACT("getIterMax", iter_max, dim_t);
   EXTRACT("getInnerIterMax", inner_iter_max, dim_t);
   EXTRACT("getDropTolerance", drop_tolerance, double);
   EXTRACT("getDropStorage", drop_storage, double);
   EXTRACT("getTruncation", truncation, dim_t);
   EXTRACT("_getRestartForC", restart, index_t);
   EXTRACT("getNumSweeps", sweeps, index_t);
   EXTRACT("getNumPreSweeps", pre_sweeps, dim_t);
   EXTRACT("getNumPostSweeps", post_sweeps, dim_t);
   EXTRACT("getLevelMax", level_max, dim_t);
   EXTRACT("getMinCoarseMatrixSize", min_coarse_matrix_size, dim_t);
   EXTRACT("getCoarseningThreshold", coarsening_threshold, double);
   EXTRACT("acceptConvergenceFailure", accept_failed_convergence, bool_t);
   EXTRACT_OPTION("getCoarsening", coarsening_method, index_t); 
   EXTRACT("getRelaxationFactor",  relaxation_factor,  double);  
  
#undef EXTRACT
#undef EXTRACT_OPTION
}
 

}  // end of namespace
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2009 by University of Queensland
5	* Earth Systems Science Computational Center (ESSCC)
6	* http://www.uq.edu.au/esscc
7	*
8	* Primary Business: Queensland, Australia
9	* Licensed under the Open Software License version 3.0
10	* http://www.opensource.org/licenses/osl-3.0.php
11	*
12	*******************************************************/
13
14
15	#ifdef PASO_MPI
16	#include <mpi.h>
17	#endif
18	#include "SystemMatrixAdapter.h"
19
20	using namespace std;
21
22	namespace finley {
23
24	struct null_deleter
25	{
26	void operator()(void const *ptr) const
27	{
28	}
29	};
30
31
32	SystemMatrixAdapter::SystemMatrixAdapter()
33	{
34	throw FinleyAdapterException("Error - Illegal to generate default SystemMatrixAdapter.");
35	}
36
37	SystemMatrixAdapter::SystemMatrixAdapter(Paso_SystemMatrix* system_matrix,
38	const int row_blocksize,
39	const escript::FunctionSpace& row_functionspace,
40	const int column_blocksize,
41	const escript::FunctionSpace& column_functionspace):
42	AbstractSystemMatrix(row_blocksize,row_functionspace,column_blocksize,column_functionspace)
43	{
44	m_system_matrix.reset(system_matrix,null_deleter());
45	}
46
47	SystemMatrixAdapter::~SystemMatrixAdapter()
48	{
49	if (m_system_matrix.unique()) {
50	Paso_SystemMatrix* mat=m_system_matrix.get();
51	Paso_SystemMatrix_free(mat);
52	}
53	}
54
55	Paso_SystemMatrix* SystemMatrixAdapter::getPaso_SystemMatrix() const
56	{
57	return m_system_matrix.get();
58	}
59
60	void SystemMatrixAdapter::ypAx(escript::Data& y,escript::Data& x) const
61	{
62	Paso_SystemMatrix* mat=getPaso_SystemMatrix();
63
64	if ( x.getDataPointSize() != getColumnBlockSize()) {
65	throw FinleyAdapterException("matrix vector product : column block size does not match the number of components in input.");
66	} else if (y.getDataPointSize() != getRowBlockSize()) {
67	throw FinleyAdapterException("matrix vector product : row block size does not match the number of components in output.");
68	} else if ( x.getFunctionSpace() != getColumnFunctionSpace()) {
69	throw FinleyAdapterException("matrix vector product : column function space and function space of input don't match.");
70	} else if (y.getFunctionSpace() != getRowFunctionSpace()) {
71	throw FinleyAdapterException("matrix vector product : row function space and function space of output don't match.");
72	}
73	x.expand();
74	y.expand();
75	x.requireWrite();
76	y.requireWrite();
77	double* x_dp=x.getSampleDataRW(0);
78	double* y_dp=y.getSampleDataRW(0);
79	Paso_SystemMatrix_MatrixVector(1., mat,x_dp, 1.,y_dp);
80	checkPasoError();
81	}
82
83	int SystemMatrixAdapter::mapOptionToPaso(const int option) {
84
85	switch (option) {
86	case ESCRIPT_DEFAULT:
87	return PASO_DEFAULT;
88	case ESCRIPT_DIRECT:
89	return PASO_DIRECT;
90	case ESCRIPT_CHOLEVSKY:
91	return PASO_CHOLEVSKY;
92	case ESCRIPT_PCG:
93	return PASO_PCG;
94	case ESCRIPT_CR:
95	return PASO_CR;
96	case ESCRIPT_CGS:
97	return PASO_CGS;
98	case ESCRIPT_BICGSTAB:
99	return PASO_BICGSTAB;
100	case ESCRIPT_SSOR:
101	return PASO_SSOR;
102	case ESCRIPT_ILU0:
103	return PASO_ILU0;
104	case ESCRIPT_ILUT:
105	return PASO_ILUT;
106	case ESCRIPT_JACOBI:
107	return PASO_JACOBI;
108	case ESCRIPT_GMRES:
109	return PASO_GMRES;
110	case ESCRIPT_PRES20:
111	return PASO_PRES20;
112	case ESCRIPT_LUMPING:
113	return PASO_LUMPING;
114	case ESCRIPT_NO_REORDERING:
115	return PASO_NO_REORDERING;
116	case ESCRIPT_MINIMUM_FILL_IN:
117	return PASO_MINIMUM_FILL_IN;
118	case ESCRIPT_NESTED_DISSECTION:
119	return PASO_NESTED_DISSECTION;
120	case ESCRIPT_MKL:
121	return PASO_MKL;
122	case ESCRIPT_UMFPACK:
123	return PASO_UMFPACK;
124	case ESCRIPT_ITERATIVE:
125	return PASO_ITERATIVE;
126	case ESCRIPT_PASO:
127	return PASO_PASO;
128	case ESCRIPT_AMG:
129	return PASO_AMG;
130	case ESCRIPT_REC_ILU:
131	return PASO_REC_ILU;
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	case ESCRIPT_MINRES:
139	return PASO_MINRES;
140	case ESCRIPT_GAUSS_SEIDEL:
141	return PASO_GAUSS_SEIDEL;
142	case ESCRIPT_RILU:
143	return PASO_RILU;
144	case ESCRIPT_DEFAULT_REORDERING:
145	return PASO_DEFAULT_REORDERING;
146	case ESCRIPT_SUPER_LU:
147	return PASO_SUPER_LU;
148	case ESCRIPT_PASTIX:
149	return PASO_PASTIX;
150	case ESCRIPT_YAIR_SHAPIRA_COARSENING:
151	return PASO_YAIR_SHAPIRA_COARSENING;
152	case ESCRIPT_RUGE_STUEBEN_COARSENING:
153	return PASO_RUGE_STUEBEN_COARSENING;
154	case ESCRIPT_AGGREGATION_COARSENING:
155	return PASO_AGGREGATION_COARSENING;
156	case ESCRIPT_NO_PRECONDITIONER:
157	return PASO_NO_PRECONDITIONER;
158	case ESCRIPT_MIN_COARSE_MATRIX_SIZE:
159	return PASO_MIN_COARSE_MATRIX_SIZE;
160	default:
161	stringstream temp;
162	temp << "Error - Cannot map option value "<< option << " onto Paso";
163	throw FinleyAdapterException(temp.str());
164	}
165	}
166
167	void finley::SystemMatrixAdapter::Print_Matrix_Info(const bool full=false) const
168	{
169	Paso_SystemMatrix* mat=m_system_matrix.get();
170	int first_row_index = mat->row_distribution->first_component[mat->mpi_info->rank];
171	int last_row_index = mat->row_distribution->first_component[mat->mpi_info->rank+1]-1;
172	int first_col_index = mat->col_distribution->first_component[mat->mpi_info->rank];
173	int last_col_index = mat->col_distribution->first_component[mat->mpi_info->rank+1]-1;
174
175	fprintf(stdout, "Print_Matrix_Info running on CPU %d of %d\n", mat->mpi_info->rank, mat->mpi_info->size);
176
177	switch (mat->type) {
178	case MATRIX_FORMAT_DEFAULT: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_DEFAULT\n"); break;
179	case MATRIX_FORMAT_CSC: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_CSC\n"); break;
180	case MATRIX_FORMAT_SYM: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_SYM\n"); break;
181	case MATRIX_FORMAT_BLK1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_BLK1\n"); break;
182	case MATRIX_FORMAT_OFFSET1: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_OFFSET1\n"); break;
183	case MATRIX_FORMAT_TRILINOS_CRS: fprintf(stdout, "\tMatrix type MATRIX_FORMAT_TRILINOS_CRS\n"); break;
184	default: fprintf(stdout, "\tMatrix type unknown\n"); break;
185	}
186
187	fprintf(stdout, "\trow indices run from %d to %d\n", first_row_index, last_row_index);
188	fprintf(stdout, "\tcol indices run from %d to %d\n", first_col_index, last_col_index);
189	fprintf(stdout, "\tmainBlock numRows %d\n", mat->mainBlock->numRows);
190	fprintf(stdout, "\tmainBlock numCols %d\n", mat->mainBlock->numCols);
191	fprintf(stdout, "\tmainBlock pattern numOutput %d\n", mat->mainBlock->pattern->numOutput);
192	fprintf(stdout, "\tcol_coupleBlock numRows %d\n", mat->col_coupleBlock->numRows);
193	fprintf(stdout, "\tcol_coupleBlock numCols %d\n", mat->col_coupleBlock->numCols);
194	fprintf(stdout, "\tcol_coupleBlock pattern numOutput %d\n", mat->col_coupleBlock->pattern->numOutput);
195	fprintf(stdout, "\trow_coupleBlock numRows %d\n", mat->row_coupleBlock->numRows);
196	fprintf(stdout, "\trow_coupleBlock numCols %d\n", mat->row_coupleBlock->numCols);
197	fprintf(stdout, "\trow_coupleBlock pattern numOutput %d\n", mat->row_coupleBlock->pattern->numOutput);
198	fprintf(stdout, "\trow_block_size %d\n", mat->row_block_size);
199	fprintf(stdout, "\tcol_block_size %d\n", mat->col_block_size);
200	fprintf(stdout, "\tblock_size %d\n", mat->block_size);
201	fprintf(stdout, "\tlogical_row_block_size %d\n", mat->logical_row_block_size);
202	fprintf(stdout, "\tlogical_col_block_size %d\n", mat->logical_col_block_size);
203
204	}
205
206	void SystemMatrixAdapter::setToSolution(escript::Data& out,escript::Data& in, boost::python::object& options) const
207	{
208	Paso_SystemMatrix* mat=getPaso_SystemMatrix();
209	Paso_Options paso_options;
210	options.attr("resetDiagnostics")();
211	escriptToPasoOptions(&paso_options,options);
212	if ( out.getDataPointSize() != getColumnBlockSize()) {
213	throw FinleyAdapterException("solve : column block size does not match the number of components of solution.");
214	} else if ( in.getDataPointSize() != getRowBlockSize()) {
215	throw FinleyAdapterException("solve : row block size does not match the number of components of right hand side.");
216	} else if ( out.getFunctionSpace() != getColumnFunctionSpace()) {
217	throw FinleyAdapterException("solve : column function space and function space of solution don't match.");
218	} else if (in.getFunctionSpace() != getRowFunctionSpace()) {
219	throw FinleyAdapterException("solve : row function space and function space of right hand side don't match.");
220	}
221	out.expand();
222	in.expand();
223	double* out_dp=out.getSampleDataRW(0);
224	double* in_dp=in.getSampleDataRW(0);
225	Paso_solve(mat,out_dp,in_dp,&paso_options);
226	pasoToEscriptOptions(&paso_options,options);
227	checkPasoError();
228	}
229
230	void SystemMatrixAdapter::nullifyRowsAndCols(escript::Data& row_q,escript::Data& col_q, const double mdv) const
231	{
232	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
233	if ( col_q.getDataPointSize() != getColumnBlockSize()) {
234	throw FinleyAdapterException("nullifyRowsAndCols : column block size does not match the number of components of column mask.");
235	} else if ( row_q.getDataPointSize() != getRowBlockSize()) {
236	throw FinleyAdapterException("nullifyRowsAndCols : row block size does not match the number of components of row mask.");
237	} else if ( col_q.getFunctionSpace() != getColumnFunctionSpace()) {
238	throw FinleyAdapterException("nullifyRowsAndCols : column function space and function space of column mask don't match.");
239	} else if (row_q.getFunctionSpace() != getRowFunctionSpace()) {
240	throw FinleyAdapterException("nullifyRowsAndCols : row function space and function space of row mask don't match.");
241	}
242	row_q.expand();
243	col_q.expand();
244	row_q.requireWrite();
245	col_q.requireWrite();
246	double* row_q_dp=row_q.getSampleDataRW(0);
247	double* col_q_dp=col_q.getSampleDataRW(0);
248	Paso_SystemMatrix_nullifyRowsAndCols(mat,row_q_dp,col_q_dp, mdv);
249	checkPasoError();
250	}
251
252	void SystemMatrixAdapter::saveMM(const std::string& fileName) const
253	{
254	if( fileName.size() == 0 )
255	{
256	throw FinleyAdapterException("Null file name!");
257	}
258
259	char *fName = TMPMEMALLOC(fileName.size()+1,char);
260
261	strcpy(fName,fileName.c_str());
262	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
263	Paso_SystemMatrix_saveMM(mat,fName);
264	checkPasoError();
265	TMPMEMFREE(fName);
266
267	}
268
269	void SystemMatrixAdapter::saveHB(const std::string& fileName) const
270	{
271	if( fileName.size() == 0 )
272	{
273	throw FinleyAdapterException("Null file name!");
274	}
275
276	char *fName = TMPMEMALLOC(fileName.size()+1,char);
277
278	strcpy(fName,fileName.c_str());
279	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
280	Paso_SystemMatrix_saveHB(mat,fName);
281	checkPasoError();
282	TMPMEMFREE(fName);
283
284	}
285
286	void SystemMatrixAdapter::resetValues() const
287	{
288	Paso_SystemMatrix* mat = getPaso_SystemMatrix();
289	Paso_SystemMatrix_setValues(mat,0.);
290	Paso_solve_free(mat);
291	checkPasoError();
292	}
293
294	void SystemMatrixAdapter::pasoToEscriptOptions(const Paso_Options* paso_options,boost::python::object& options)
295	{
296	#define SET(__key__,__val__,__type__) options.attr("_updateDiagnostics")(__key__,(__type__)paso_options->__val__)
297	SET("num_iter", num_iter, int);
298	SET("num_level", num_level, int);
299	SET("num_inner_iter", num_inner_iter, int);
300	SET("time", time, double);
301	SET("set_up_time", set_up_time, double);
302	SET("net_time", net_time, double);
303	SET("residual_norm", residual_norm, double);
304	SET("converged",converged, bool);
305	#undef SET
306	}
307	void SystemMatrixAdapter::escriptToPasoOptions(Paso_Options* paso_options, const boost::python::object& options)
308	{
309	#define EXTRACT(__key__,__val__,__type__) paso_options->__val__=boost::python::extract<__type__>(options.attr(__key__)())
310	#define EXTRACT_OPTION(__key__,__val__,__type__) paso_options->__val__=mapOptionToPaso(boost::python::extract<__type__>(options.attr(__key__)()))
311
312	Paso_Options_setDefaults(paso_options);
313	EXTRACT_OPTION("getSolverMethod",method,index_t);
314	EXTRACT_OPTION("getPackage",package,index_t);
315	EXTRACT("isVerbose",verbose,bool_t);
316	EXTRACT("isSymmetric",symmetric,bool_t);
317	EXTRACT("getTolerance",tolerance, double);
318	EXTRACT("getAbsoluteTolerance",absolute_tolerance, double);
319	EXTRACT("getInnerTolerance",inner_tolerance, double);
320	EXTRACT("adaptInnerTolerance",adapt_inner_tolerance, bool_t);
321	EXTRACT_OPTION("getReordering", reordering, index_t);
322	EXTRACT_OPTION("getPreconditioner", preconditioner, index_t);
323	EXTRACT("getIterMax", iter_max, dim_t);
324	EXTRACT("getInnerIterMax", inner_iter_max, dim_t);
325	EXTRACT("getDropTolerance", drop_tolerance, double);
326	EXTRACT("getDropStorage", drop_storage, double);
327	EXTRACT("getTruncation", truncation, dim_t);
328	EXTRACT("_getRestartForC", restart, index_t);
329	EXTRACT("getNumSweeps", sweeps, index_t);
330	EXTRACT("getNumPreSweeps", pre_sweeps, dim_t);
331	EXTRACT("getNumPostSweeps", post_sweeps, dim_t);
332	EXTRACT("getLevelMax", level_max, dim_t);
333	EXTRACT("getMinCoarseMatrixSize", min_coarse_matrix_size, dim_t);
334	EXTRACT("getCoarseningThreshold", coarsening_threshold, double);
335	EXTRACT("acceptConvergenceFailure", accept_failed_convergence, bool_t);
336	EXTRACT_OPTION("getCoarsening", coarsening_method, index_t);
337	EXTRACT("getRelaxationFactor", relaxation_factor, double);
338
339	#undef EXTRACT
340	#undef EXTRACT_OPTION
341	}
342
343
344	} // end of namespace
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision