[Scipy-svn] r3613 - trunk/scipy/sandbox/multigrid/multigridtools

scipy-svn at scipy.org scipy-svn at scipy.org
Mon Dec 3 18:30:08 EST 2007


Author: wnbell
Date: 2007-12-03 17:30:01 -0600 (Mon, 03 Dec 2007)
New Revision: 3613

Modified:
   trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i
   trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.py
   trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx
   trunk/scipy/sandbox/multigrid/multigridtools/ruge_stuben.h
   trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
Log:
templated multigridtools index types
removed 64-bit indices


Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i	2007-12-03 22:01:43 UTC (rev 3612)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i	2007-12-03 23:30:01 UTC (rev 3613)
@@ -61,13 +61,7 @@
 %enddef
 
 
-I_IN_ARRAY1( int         )
-I_IN_ARRAY1( long long   )
-T_IN_ARRAY1( float       )
-T_IN_ARRAY1( double      )
 
-
-
  /*
   * OUT types
   */
@@ -103,13 +97,8 @@
 };
 %enddef
 
-I_ARRAY_ARGOUT( int    )
-I_ARRAY_ARGOUT( long long    )
-T_ARRAY_ARGOUT( float  )
-T_ARRAY_ARGOUT( double )
 
 
-
  /*
   * INPLACE types
   */
@@ -126,12 +115,32 @@
 };
 %enddef
 
-I_INPLACE_ARRAY1( int    )
-I_INPLACE_ARRAY1( long long    )
-T_INPLACE_ARRAY1( float  )
-T_INPLACE_ARRAY1( double )
 
 
+/*
+ * Macros to instantiate index types and data types
+ */
+%define DECLARE_INDEX_TYPE( ctype )
+I_IN_ARRAY1( ctype )
+I_ARRAY_ARGOUT( ctype )
+I_INPLACE_ARRAY1( ctype )
+%enddef
+
+%define DECLARE_DATA_TYPE( ctype )
+T_IN_ARRAY1( ctype )
+T_ARRAY_ARGOUT( ctype )
+T_INPLACE_ARRAY1( ctype )
+%enddef
+
+/*
+ * Create all desired index and data types here
+ */
+DECLARE_INDEX_TYPE( int       )
+
+DECLARE_DATA_TYPE( float               )
+DECLARE_DATA_TYPE( double              )
+
+
 %include "ruge_stuben.h"
 %include "smoothed_aggregation.h"
 %include "relaxation.h"
@@ -143,14 +152,11 @@
 %define INSTANTIATE_BOTH( f_name )
 %template(f_name)   f_name<int,float>;
 %template(f_name)   f_name<int,double>;
-%template(f_name)   f_name<long long,float>;
-%template(f_name)   f_name<long long,double>;
 /* 64-bit indices would go here */
 %enddef
  
 %define INSTANTIATE_INDEX( f_name )
 %template(f_name)   f_name<int>;
-%template(f_name)   f_name<long long>;
 %enddef
 
 %define INSTANTIATE_DATA( f_name )
@@ -159,14 +165,11 @@
 %enddef
  
  
+INSTANTIATE_INDEX(sa_get_aggregates)
 
-INSTANTIATE_DATA(rs_strong_connections)
-INSTANTIATE_DATA(rs_interpolation)
-
-INSTANTIATE_DATA(sa_strong_connections)
-INSTANTIATE_DATA(sa_smoother)
-/*INSTANTIATE_INDEX(sa_get_aggregates)*/
-
+INSTANTIATE_BOTH(rs_strong_connections)
+INSTANTIATE_BOTH(rs_interpolation)
+INSTANTIATE_BOTH(sa_strong_connections)
 INSTANTIATE_BOTH(gauss_seidel)
 INSTANTIATE_BOTH(jacobi)
 

Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.py	2007-12-03 22:01:43 UTC (rev 3612)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.py	2007-12-03 23:30:01 UTC (rev 3613)
@@ -1,5 +1,5 @@
 # This file was automatically generated by SWIG (http://www.swig.org).
-# Version 1.3.32
+# Version 1.3.34
 #
 # Don't modify this file, modify the SWIG interface instead.
 # This file is compatible with both classic and new-style classes.
@@ -53,84 +53,60 @@
 F_NODE = _multigridtools.F_NODE
 
 def sa_get_aggregates(*args):
-    """sa_get_aggregates(int n_row, int Ap, int Aj, std::vector<(int)> Bj)"""
-    return _multigridtools.sa_get_aggregates(*args)
+  """sa_get_aggregates(int n_row, int Ap, int Aj, std::vector<(int)> Bj)"""
+  return _multigridtools.sa_get_aggregates(*args)
 
 
 def rs_strong_connections(*args):
+  """
+    rs_strong_connections(int n_row, float theta, int Ap, int Aj, float Ax, std::vector<(int)> Sp, 
+        std::vector<(int)> Sj, 
+        std::vector<(float)> Sx)
+    rs_strong_connections(int n_row, double theta, int Ap, int Aj, double Ax, 
+        std::vector<(int)> Sp, std::vector<(int)> Sj, 
+        std::vector<(double)> Sx)
     """
-      rs_strong_connections(int n_row, float theta, int Ap, int Aj, float Ax, std::vector<(int)> Sp,
-          std::vector<(int)> Sj,
-          std::vector<(float)> Sx)
-      rs_strong_connections(int n_row, double theta, int Ap, int Aj, double Ax,
-          std::vector<(int)> Sp, std::vector<(int)> Sj,
-          std::vector<(double)> Sx)
-      """
-    return _multigridtools.rs_strong_connections(*args)
+  return _multigridtools.rs_strong_connections(*args)
 
 def rs_interpolation(*args):
+  """
+    rs_interpolation(int n_nodes, int Ap, int Aj, float Ax, int Sp, int Sj, 
+        float Sx, int Tp, int Tj, float Tx, std::vector<(int)> Bp, 
+        std::vector<(int)> Bj, std::vector<(float)> Bx)
+    rs_interpolation(int n_nodes, int Ap, int Aj, double Ax, int Sp, int Sj, 
+        double Sx, int Tp, int Tj, double Tx, std::vector<(int)> Bp, 
+        std::vector<(int)> Bj, std::vector<(double)> Bx)
     """
-      rs_interpolation(int n_nodes, int Ap, int Aj, float Ax, int Sp, int Sj,
-          float Sx, int Tp, int Tj, float Tx, std::vector<(int)> Bp,
-          std::vector<(int)> Bj, std::vector<(float)> Bx)
-      rs_interpolation(int n_nodes, int Ap, int Aj, double Ax, int Sp, int Sj,
-          double Sx, int Tp, int Tj, double Tx, std::vector<(int)> Bp,
-          std::vector<(int)> Bj, std::vector<(double)> Bx)
-      """
-    return _multigridtools.rs_interpolation(*args)
+  return _multigridtools.rs_interpolation(*args)
 
 def sa_strong_connections(*args):
+  """
+    sa_strong_connections(int n_row, float epsilon, int Ap, int Aj, float Ax, 
+        std::vector<(int)> Sp, std::vector<(int)> Sj, 
+        std::vector<(float)> Sx)
+    sa_strong_connections(int n_row, double epsilon, int Ap, int Aj, double Ax, 
+        std::vector<(int)> Sp, std::vector<(int)> Sj, 
+        std::vector<(double)> Sx)
     """
-      sa_strong_connections(int n_row, float epsilon, int Ap, int Aj, float Ax,
-          std::vector<(int)> Sp, std::vector<(int)> Sj,
-          std::vector<(float)> Sx)
-      sa_strong_connections(int n_row, double epsilon, int Ap, int Aj, double Ax,
-          std::vector<(int)> Sp, std::vector<(int)> Sj,
-          std::vector<(double)> Sx)
-      """
-    return _multigridtools.sa_strong_connections(*args)
+  return _multigridtools.sa_strong_connections(*args)
 
-def sa_smoother(*args):
-    """
-      sa_smoother(int n_row, float omega, int Ap, int Aj, float Ax, int Sp,
-          int Sj, float Sx, std::vector<(int)> Bp,
-          std::vector<(int)> Bj, std::vector<(float)> Bx)
-      sa_smoother(int n_row, double omega, int Ap, int Aj, double Ax,
-          int Sp, int Sj, double Sx, std::vector<(int)> Bp,
-          std::vector<(int)> Bj, std::vector<(double)> Bx)
-      """
-    return _multigridtools.sa_smoother(*args)
-
 def gauss_seidel(*args):
+  """
+    gauss_seidel(int n_row, int Ap, int Aj, float Ax, float x, float b, 
+        int row_start, int row_stop, int row_step)
+    gauss_seidel(int n_row, int Ap, int Aj, double Ax, double x, double b, 
+        int row_start, int row_stop, int row_step)
     """
-      gauss_seidel(int n_row, int Ap, int Aj, float Ax, float x, float b,
-          int row_start, int row_stop, int row_step)
-      gauss_seidel(int n_row, int Ap, int Aj, double Ax, double x, double b,
-          int row_start, int row_stop, int row_step)
-      gauss_seidel(long long n_row, long long Ap, long long Aj, float Ax,
-          float x, float b, long long row_start, long long row_stop,
-          long long row_step)
-      gauss_seidel(long long n_row, long long Ap, long long Aj, double Ax,
-          double x, double b, long long row_start,
-          long long row_stop, long long row_step)
-      """
-    return _multigridtools.gauss_seidel(*args)
+  return _multigridtools.gauss_seidel(*args)
 
 def jacobi(*args):
+  """
+    jacobi(int n_row, int Ap, int Aj, float Ax, float x, float b, 
+        float temp, int row_start, int row_stop, 
+        int row_step, float omega)
+    jacobi(int n_row, int Ap, int Aj, double Ax, double x, double b, 
+        double temp, int row_start, int row_stop, 
+        int row_step, double omega)
     """
-      jacobi(int n_row, int Ap, int Aj, float Ax, float x, float b,
-          float temp, int row_start, int row_stop,
-          int row_step, float omega)
-      jacobi(int n_row, int Ap, int Aj, double Ax, double x, double b,
-          double temp, int row_start, int row_stop,
-          int row_step, double omega)
-      jacobi(long long n_row, long long Ap, long long Aj, float Ax,
-          float x, float b, float temp, long long row_start,
-          long long row_stop, long long row_step,
-          float omega)
-      jacobi(long long n_row, long long Ap, long long Aj, double Ax,
-          double x, double b, double temp, long long row_start,
-          long long row_stop, long long row_step,
-          double omega)
-      """
-    return _multigridtools.jacobi(*args)
+  return _multigridtools.jacobi(*args)
+

Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx	2007-12-03 22:01:43 UTC (rev 3612)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx	2007-12-03 23:30:01 UTC (rev 3613)
@@ -1,6 +1,6 @@
 /* ----------------------------------------------------------------------------
  * This file was automatically generated by SWIG (http://www.swig.org).
- * Version 1.3.32
+ * Version 1.3.34
  * 
  * This file is not intended to be easily readable and contains a number of 
  * coding conventions designed to improve portability and efficiency. Do not make
@@ -1124,7 +1124,7 @@
       return 1;
     } else {
       PyErr_Format(PyExc_TypeError, "%s expected %s%d arguments, got none", 
-		   name, (min == max ? "" : "at least "), min);
+		   name, (min == max ? "" : "at least "), (int)min);
       return 0;
     }
   }  
@@ -1135,11 +1135,11 @@
     register Py_ssize_t l = PyTuple_GET_SIZE(args);
     if (l < min) {
       PyErr_Format(PyExc_TypeError, "%s expected %s%d arguments, got %d", 
-		   name, (min == max ? "" : "at least "), min, l);
+		   name, (min == max ? "" : "at least "), (int)min, (int)l);
       return 0;
     } else if (l > max) {
       PyErr_Format(PyExc_TypeError, "%s expected %s%d arguments, got %d", 
-		   name, (min == max ? "" : "at most "), max, l);
+		   name, (min == max ? "" : "at most "), (int)max, (int)l);
       return 0;
     } else {
       register int i;
@@ -2501,7 +2501,7 @@
 
 #define SWIG_name    "_multigridtools"
 
-#define SWIGVERSION 0x010332 
+#define SWIGVERSION 0x010334 
 #define SWIG_VERSION SWIGVERSION
 
 
@@ -3099,43 +3099,6 @@
   return res;
 }
 
-
-SWIGINTERN int
-SWIG_AsVal_long_SS_long (PyObject *obj, long long *val)
-{
-  int res = SWIG_TypeError;
-  if (PyLong_Check(obj)) {
-    long long v = PyLong_AsLongLong(obj);
-    if (!PyErr_Occurred()) {
-      if (val) *val = v;
-      return SWIG_OK;
-    } else {
-      PyErr_Clear();
-    }
-  } else {
-    long v;
-    res = SWIG_AsVal_long (obj,&v);
-    if (SWIG_IsOK(res)) {
-      if (val) *val = v;
-      return res;
-    }
-  }
-#ifdef SWIG_PYTHON_CAST_MODE
-  {
-    const double mant_max = 1LL << DBL_MANT_DIG;
-    const double mant_min = -mant_max;
-    double d;
-    res = SWIG_AsVal_double (obj,&d);
-    if (SWIG_IsOK(res) && SWIG_CanCastAsInteger(&d, mant_min, mant_max)) {
-      if (val) *val = (long long)(d);
-      return SWIG_AddCast(res);
-    }
-    res = SWIG_TypeError;
-  }
-#endif
-  return res;
-}
-
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -3186,7 +3149,7 @@
     
     arg3 = (int*) array3->data;
   }
-  sa_get_aggregates(arg1,(int const (*))arg2,(int const (*))arg3,arg4);
+  sa_get_aggregates<int >(arg1,(int const (*))arg2,(int const (*))arg3,arg4);
   resultobj = SWIG_Py_Void();
   {
     int length = (arg4)->size(); 
@@ -3295,7 +3258,7 @@
     
     arg5 = (float*) array5->data;
   }
-  rs_strong_connections<float >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(float const (*))arg5,arg6,arg7,arg8);
+  rs_strong_connections<int,float >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(float const (*))arg5,arg6,arg7,arg8);
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
@@ -3424,7 +3387,7 @@
     
     arg5 = (double*) array5->data;
   }
-  rs_strong_connections<double >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(double const (*))arg5,arg6,arg7,arg8);
+  rs_strong_connections<int,double >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(double const (*))arg5,arg6,arg7,arg8);
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
@@ -3545,7 +3508,7 @@
   }
   
 fail:
-  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'rs_strong_connections'.\n  Possible C/C++ prototypes are:\n""    rs_strong_connections<(float)>(int const,float const,int const [],int const [],float const [],std::vector<int > *,std::vector<int > *,std::vector<float > *)\n""    rs_strong_connections<(double)>(int const,double const,int const [],int const [],double const [],std::vector<int > *,std::vector<int > *,std::vector<double > *)\n");
+  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'rs_strong_connections'.\n  Possible C/C++ prototypes are:\n""    rs_strong_connections<(int,float)>(int const,float const,int const [],int const [],float const [],std::vector<int > *,std::vector<int > *,std::vector<float > *)\n""    rs_strong_connections<(int,double)>(int const,double const,int const [],int const [],double const [],std::vector<int > *,std::vector<int > *,std::vector<double > *)\n");
   return NULL;
 }
 
@@ -3707,7 +3670,7 @@
     
     arg10 = (float*) array10->data;
   }
-  rs_interpolation<float >(arg1,(int const (*))arg2,(int const (*))arg3,(float const (*))arg4,(int const (*))arg5,(int const (*))arg6,(float const (*))arg7,(int const (*))arg8,(int const (*))arg9,(float const (*))arg10,arg11,arg12,arg13);
+  rs_interpolation<int,float >(arg1,(int const (*))arg2,(int const (*))arg3,(float const (*))arg4,(int const (*))arg5,(int const (*))arg6,(float const (*))arg7,(int const (*))arg8,(int const (*))arg9,(float const (*))arg10,arg11,arg12,arg13);
   resultobj = SWIG_Py_Void();
   {
     int length = (arg11)->size(); 
@@ -3947,7 +3910,7 @@
     
     arg10 = (double*) array10->data;
   }
-  rs_interpolation<double >(arg1,(int const (*))arg2,(int const (*))arg3,(double const (*))arg4,(int const (*))arg5,(int const (*))arg6,(double const (*))arg7,(int const (*))arg8,(int const (*))arg9,(double const (*))arg10,arg11,arg12,arg13);
+  rs_interpolation<int,double >(arg1,(int const (*))arg2,(int const (*))arg3,(double const (*))arg4,(int const (*))arg5,(int const (*))arg6,(double const (*))arg7,(int const (*))arg8,(int const (*))arg9,(double const (*))arg10,arg11,arg12,arg13);
   resultobj = SWIG_Py_Void();
   {
     int length = (arg11)->size(); 
@@ -4152,7 +4115,7 @@
   }
   
 fail:
-  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'rs_interpolation'.\n  Possible C/C++ prototypes are:\n""    rs_interpolation<(float)>(int const,int const [],int const [],float const [],int const [],int const [],float const [],int const [],int const [],float const [],std::vector<int > *,std::vector<int > *,std::vector<float > *)\n""    rs_interpolation<(double)>(int const,int const [],int const [],double const [],int const [],int const [],double const [],int const [],int const [],double const [],std::vector<int > *,std::vector<int > *,std::vector<double > *)\n");
+  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'rs_interpolation'.\n  Possible C/C++ prototypes are:\n""    rs_interpolation<(int,float)>(int const,int const [],int const [],float const [],int const [],int const [],float const [],int const [],int const [],float const [],std::vector<int > *,std::vector<int > *,std::vector<float > *)\n""    rs_interpolation<(int,double)>(int const,int const [],int const [],double const [],int const [],int const [],double const [],int const [],int const [],double const [],std::vector<int > *,std::vector<int > *,std::vector<double > *)\n");
   return NULL;
 }
 
@@ -4239,7 +4202,7 @@
     
     arg5 = (float*) array5->data;
   }
-  sa_strong_connections<float >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(float const (*))arg5,arg6,arg7,arg8);
+  sa_strong_connections<int,float >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(float const (*))arg5,arg6,arg7,arg8);
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
@@ -4368,7 +4331,7 @@
     
     arg5 = (double*) array5->data;
   }
-  sa_strong_connections<double >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(double const (*))arg5,arg6,arg7,arg8);
+  sa_strong_connections<int,double >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(double const (*))arg5,arg6,arg7,arg8);
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
@@ -4489,498 +4452,11 @@
   }
   
 fail:
-  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'sa_strong_connections'.\n  Possible C/C++ prototypes are:\n""    sa_strong_connections<(float)>(int const,float const,int const [],int const [],float const [],std::vector<int > *,std::vector<int > *,std::vector<float > *)\n""    sa_strong_connections<(double)>(int const,double const,int const [],int const [],double const [],std::vector<int > *,std::vector<int > *,std::vector<double > *)\n");
+  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'sa_strong_connections'.\n  Possible C/C++ prototypes are:\n""    sa_strong_connections<(int,float)>(int const,float const,int const [],int const [],float const [],std::vector<int > *,std::vector<int > *,std::vector<float > *)\n""    sa_strong_connections<(int,double)>(int const,double const,int const [],int const [],double const [],std::vector<int > *,std::vector<int > *,std::vector<double > *)\n");
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_sa_smoother__SWIG_1(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  int arg1 ;
-  float arg2 ;
-  int *arg3 ;
-  int *arg4 ;
-  float *arg5 ;
-  int *arg6 ;
-  int *arg7 ;
-  float *arg8 ;
-  std::vector<int > *arg9 = (std::vector<int > *) 0 ;
-  std::vector<int > *arg10 = (std::vector<int > *) 0 ;
-  std::vector<float > *arg11 = (std::vector<float > *) 0 ;
-  int val1 ;
-  int ecode1 = 0 ;
-  float val2 ;
-  int ecode2 = 0 ;
-  PyArrayObject *array3 = NULL ;
-  int is_new_object3 ;
-  PyArrayObject *array4 = NULL ;
-  int is_new_object4 ;
-  PyArrayObject *array5 = NULL ;
-  int is_new_object5 ;
-  PyArrayObject *array6 = NULL ;
-  int is_new_object6 ;
-  PyArrayObject *array7 = NULL ;
-  int is_new_object7 ;
-  PyArrayObject *array8 = NULL ;
-  int is_new_object8 ;
-  std::vector<int > *tmp9 ;
-  std::vector<int > *tmp10 ;
-  std::vector<float > *tmp11 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  PyObject * obj2 = 0 ;
-  PyObject * obj3 = 0 ;
-  PyObject * obj4 = 0 ;
-  PyObject * obj5 = 0 ;
-  PyObject * obj6 = 0 ;
-  PyObject * obj7 = 0 ;
-  
-  {
-    tmp9 = new std::vector<int>(); 
-    arg9 = tmp9; 
-  }
-  {
-    tmp10 = new std::vector<int>(); 
-    arg10 = tmp10; 
-  }
-  {
-    tmp11 = new std::vector<float>(); 
-    arg11 = tmp11; 
-  }
-  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOO:sa_smoother",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7)) SWIG_fail;
-  ecode1 = SWIG_AsVal_int(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "sa_smoother" "', argument " "1"" of type '" "int""'");
-  } 
-  arg1 = static_cast< int >(val1);
-  ecode2 = SWIG_AsVal_float(obj1, &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "sa_smoother" "', argument " "2"" of type '" "float""'");
-  } 
-  arg2 = static_cast< float >(val2);
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_INT, &is_new_object3);
-    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
-      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
-    
-    arg3 = (int*) array3->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array4 = obj_to_array_contiguous_allow_conversion(obj3, PyArray_INT, &is_new_object4);
-    if (!array4 || !require_dimensions(array4,1) || !require_size(array4,size,1)
-      || !require_contiguous(array4)   || !require_native(array4)) SWIG_fail;
-    
-    arg4 = (int*) array4->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array5 = obj_to_array_contiguous_allow_conversion(obj4, PyArray_FLOAT, &is_new_object5);
-    if (!array5 || !require_dimensions(array5,1) || !require_size(array5,size,1)
-      || !require_contiguous(array5)   || !require_native(array5)) SWIG_fail;
-    
-    arg5 = (float*) array5->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array6 = obj_to_array_contiguous_allow_conversion(obj5, PyArray_INT, &is_new_object6);
-    if (!array6 || !require_dimensions(array6,1) || !require_size(array6,size,1)
-      || !require_contiguous(array6)   || !require_native(array6)) SWIG_fail;
-    
-    arg6 = (int*) array6->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array7 = obj_to_array_contiguous_allow_conversion(obj6, PyArray_INT, &is_new_object7);
-    if (!array7 || !require_dimensions(array7,1) || !require_size(array7,size,1)
-      || !require_contiguous(array7)   || !require_native(array7)) SWIG_fail;
-    
-    arg7 = (int*) array7->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array8 = obj_to_array_contiguous_allow_conversion(obj7, PyArray_FLOAT, &is_new_object8);
-    if (!array8 || !require_dimensions(array8,1) || !require_size(array8,size,1)
-      || !require_contiguous(array8)   || !require_native(array8)) SWIG_fail;
-    
-    arg8 = (float*) array8->data;
-  }
-  sa_smoother<float >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(float const (*))arg5,(int const (*))arg6,(int const (*))arg7,(float const (*))arg8,arg9,arg10,arg11);
-  resultobj = SWIG_Py_Void();
-  {
-    int length = (arg9)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
-    memcpy(PyArray_DATA(obj),&((*(arg9))[0]),sizeof(int)*length);	 
-    delete arg9; 
-    resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
-  }
-  {
-    int length = (arg10)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
-    memcpy(PyArray_DATA(obj),&((*(arg10))[0]),sizeof(int)*length);	 
-    delete arg10; 
-    resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
-  }
-  {
-    int length = (arg11)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length,PyArray_FLOAT); 
-    memcpy(PyArray_DATA(obj),&((*(arg11))[0]),sizeof(float)*length);	 
-    delete arg11; 
-    resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object5 && array5) Py_DECREF(array5);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  {
-    if (is_new_object7 && array7) Py_DECREF(array7);
-  }
-  {
-    if (is_new_object8 && array8) Py_DECREF(array8);
-  }
-  return resultobj;
-fail:
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object5 && array5) Py_DECREF(array5);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  {
-    if (is_new_object7 && array7) Py_DECREF(array7);
-  }
-  {
-    if (is_new_object8 && array8) Py_DECREF(array8);
-  }
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_sa_smoother__SWIG_2(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  int arg1 ;
-  double arg2 ;
-  int *arg3 ;
-  int *arg4 ;
-  double *arg5 ;
-  int *arg6 ;
-  int *arg7 ;
-  double *arg8 ;
-  std::vector<int > *arg9 = (std::vector<int > *) 0 ;
-  std::vector<int > *arg10 = (std::vector<int > *) 0 ;
-  std::vector<double > *arg11 = (std::vector<double > *) 0 ;
-  int val1 ;
-  int ecode1 = 0 ;
-  double val2 ;
-  int ecode2 = 0 ;
-  PyArrayObject *array3 = NULL ;
-  int is_new_object3 ;
-  PyArrayObject *array4 = NULL ;
-  int is_new_object4 ;
-  PyArrayObject *array5 = NULL ;
-  int is_new_object5 ;
-  PyArrayObject *array6 = NULL ;
-  int is_new_object6 ;
-  PyArrayObject *array7 = NULL ;
-  int is_new_object7 ;
-  PyArrayObject *array8 = NULL ;
-  int is_new_object8 ;
-  std::vector<int > *tmp9 ;
-  std::vector<int > *tmp10 ;
-  std::vector<double > *tmp11 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  PyObject * obj2 = 0 ;
-  PyObject * obj3 = 0 ;
-  PyObject * obj4 = 0 ;
-  PyObject * obj5 = 0 ;
-  PyObject * obj6 = 0 ;
-  PyObject * obj7 = 0 ;
-  
-  {
-    tmp9 = new std::vector<int>(); 
-    arg9 = tmp9; 
-  }
-  {
-    tmp10 = new std::vector<int>(); 
-    arg10 = tmp10; 
-  }
-  {
-    tmp11 = new std::vector<double>(); 
-    arg11 = tmp11; 
-  }
-  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOO:sa_smoother",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7)) SWIG_fail;
-  ecode1 = SWIG_AsVal_int(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "sa_smoother" "', argument " "1"" of type '" "int""'");
-  } 
-  arg1 = static_cast< int >(val1);
-  ecode2 = SWIG_AsVal_double(obj1, &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "sa_smoother" "', argument " "2"" of type '" "double""'");
-  } 
-  arg2 = static_cast< double >(val2);
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_INT, &is_new_object3);
-    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
-      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
-    
-    arg3 = (int*) array3->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array4 = obj_to_array_contiguous_allow_conversion(obj3, PyArray_INT, &is_new_object4);
-    if (!array4 || !require_dimensions(array4,1) || !require_size(array4,size,1)
-      || !require_contiguous(array4)   || !require_native(array4)) SWIG_fail;
-    
-    arg4 = (int*) array4->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array5 = obj_to_array_contiguous_allow_conversion(obj4, PyArray_DOUBLE, &is_new_object5);
-    if (!array5 || !require_dimensions(array5,1) || !require_size(array5,size,1)
-      || !require_contiguous(array5)   || !require_native(array5)) SWIG_fail;
-    
-    arg5 = (double*) array5->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array6 = obj_to_array_contiguous_allow_conversion(obj5, PyArray_INT, &is_new_object6);
-    if (!array6 || !require_dimensions(array6,1) || !require_size(array6,size,1)
-      || !require_contiguous(array6)   || !require_native(array6)) SWIG_fail;
-    
-    arg6 = (int*) array6->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array7 = obj_to_array_contiguous_allow_conversion(obj6, PyArray_INT, &is_new_object7);
-    if (!array7 || !require_dimensions(array7,1) || !require_size(array7,size,1)
-      || !require_contiguous(array7)   || !require_native(array7)) SWIG_fail;
-    
-    arg7 = (int*) array7->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array8 = obj_to_array_contiguous_allow_conversion(obj7, PyArray_DOUBLE, &is_new_object8);
-    if (!array8 || !require_dimensions(array8,1) || !require_size(array8,size,1)
-      || !require_contiguous(array8)   || !require_native(array8)) SWIG_fail;
-    
-    arg8 = (double*) array8->data;
-  }
-  sa_smoother<double >(arg1,arg2,(int const (*))arg3,(int const (*))arg4,(double const (*))arg5,(int const (*))arg6,(int const (*))arg7,(double const (*))arg8,arg9,arg10,arg11);
-  resultobj = SWIG_Py_Void();
-  {
-    int length = (arg9)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
-    memcpy(PyArray_DATA(obj),&((*(arg9))[0]),sizeof(int)*length);	 
-    delete arg9; 
-    resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
-  }
-  {
-    int length = (arg10)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
-    memcpy(PyArray_DATA(obj),&((*(arg10))[0]),sizeof(int)*length);	 
-    delete arg10; 
-    resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
-  }
-  {
-    int length = (arg11)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length,PyArray_DOUBLE); 
-    memcpy(PyArray_DATA(obj),&((*(arg11))[0]),sizeof(double)*length);	 
-    delete arg11; 
-    resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object5 && array5) Py_DECREF(array5);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  {
-    if (is_new_object7 && array7) Py_DECREF(array7);
-  }
-  {
-    if (is_new_object8 && array8) Py_DECREF(array8);
-  }
-  return resultobj;
-fail:
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object5 && array5) Py_DECREF(array5);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  {
-    if (is_new_object7 && array7) Py_DECREF(array7);
-  }
-  {
-    if (is_new_object8 && array8) Py_DECREF(array8);
-  }
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_sa_smoother(PyObject *self, PyObject *args) {
-  int argc;
-  PyObject *argv[9];
-  int ii;
-  
-  if (!PyTuple_Check(args)) SWIG_fail;
-  argc = (int)PyObject_Length(args);
-  for (ii = 0; (ii < argc) && (ii < 8); ii++) {
-    argv[ii] = PyTuple_GET_ITEM(args,ii);
-  }
-  if (argc == 8) {
-    int _v;
-    {
-      int res = SWIG_AsVal_int(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      {
-        int res = SWIG_AsVal_float(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_INT)) ? 1 : 0;
-        }
-        if (_v) {
-          {
-            _v = (is_array(argv[3]) && PyArray_CanCastSafely(PyArray_TYPE(argv[3]),PyArray_INT)) ? 1 : 0;
-          }
-          if (_v) {
-            {
-              _v = (is_array(argv[4]) && PyArray_CanCastSafely(PyArray_TYPE(argv[4]),PyArray_FLOAT)) ? 1 : 0;
-            }
-            if (_v) {
-              {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_INT)) ? 1 : 0;
-              }
-              if (_v) {
-                {
-                  _v = (is_array(argv[6]) && PyArray_CanCastSafely(PyArray_TYPE(argv[6]),PyArray_INT)) ? 1 : 0;
-                }
-                if (_v) {
-                  {
-                    _v = (is_array(argv[7]) && PyArray_CanCastSafely(PyArray_TYPE(argv[7]),PyArray_FLOAT)) ? 1 : 0;
-                  }
-                  if (_v) {
-                    return _wrap_sa_smoother__SWIG_1(self, args);
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-  if (argc == 8) {
-    int _v;
-    {
-      int res = SWIG_AsVal_int(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      {
-        int res = SWIG_AsVal_double(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_INT)) ? 1 : 0;
-        }
-        if (_v) {
-          {
-            _v = (is_array(argv[3]) && PyArray_CanCastSafely(PyArray_TYPE(argv[3]),PyArray_INT)) ? 1 : 0;
-          }
-          if (_v) {
-            {
-              _v = (is_array(argv[4]) && PyArray_CanCastSafely(PyArray_TYPE(argv[4]),PyArray_DOUBLE)) ? 1 : 0;
-            }
-            if (_v) {
-              {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_INT)) ? 1 : 0;
-              }
-              if (_v) {
-                {
-                  _v = (is_array(argv[6]) && PyArray_CanCastSafely(PyArray_TYPE(argv[6]),PyArray_INT)) ? 1 : 0;
-                }
-                if (_v) {
-                  {
-                    _v = (is_array(argv[7]) && PyArray_CanCastSafely(PyArray_TYPE(argv[7]),PyArray_DOUBLE)) ? 1 : 0;
-                  }
-                  if (_v) {
-                    return _wrap_sa_smoother__SWIG_2(self, args);
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-  
-fail:
-  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'sa_smoother'.\n  Possible C/C++ prototypes are:\n""    sa_smoother<(float)>(int const,float const,int const [],int const [],float const [],int const [],int const [],float const [],std::vector<int > *,std::vector<int > *,std::vector<float > *)\n""    sa_smoother<(double)>(int const,double const,int const [],int const [],double const [],int const [],int const [],double const [],std::vector<int > *,std::vector<int > *,std::vector<double > *)\n");
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_gauss_seidel__SWIG_1(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   int arg1 ;
@@ -5253,278 +4729,6 @@
 }
 
 
-SWIGINTERN PyObject *_wrap_gauss_seidel__SWIG_3(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  long long arg1 ;
-  long long *arg2 ;
-  long long *arg3 ;
-  float *arg4 ;
-  float *arg5 ;
-  float *arg6 ;
-  long long arg7 ;
-  long long arg8 ;
-  long long arg9 ;
-  long long val1 ;
-  int ecode1 = 0 ;
-  PyArrayObject *array2 = NULL ;
-  int is_new_object2 ;
-  PyArrayObject *array3 = NULL ;
-  int is_new_object3 ;
-  PyArrayObject *array4 = NULL ;
-  int is_new_object4 ;
-  PyArrayObject *temp5 = NULL ;
-  PyArrayObject *array6 = NULL ;
-  int is_new_object6 ;
-  long long val7 ;
-  int ecode7 = 0 ;
-  long long val8 ;
-  int ecode8 = 0 ;
-  long long val9 ;
-  int ecode9 = 0 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  PyObject * obj2 = 0 ;
-  PyObject * obj3 = 0 ;
-  PyObject * obj4 = 0 ;
-  PyObject * obj5 = 0 ;
-  PyObject * obj6 = 0 ;
-  PyObject * obj7 = 0 ;
-  PyObject * obj8 = 0 ;
-  
-  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOO:gauss_seidel",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8)) SWIG_fail;
-  ecode1 = SWIG_AsVal_long_SS_long(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "gauss_seidel" "', argument " "1"" of type '" "long long""'");
-  } 
-  arg1 = static_cast< long long >(val1);
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_LONGLONG, &is_new_object2);
-    if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
-      || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
-    
-    arg2 = (long long*) array2->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_LONGLONG, &is_new_object3);
-    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
-      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
-    
-    arg3 = (long long*) array3->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array4 = obj_to_array_contiguous_allow_conversion(obj3, PyArray_FLOAT, &is_new_object4);
-    if (!array4 || !require_dimensions(array4,1) || !require_size(array4,size,1)
-      || !require_contiguous(array4)   || !require_native(array4)) SWIG_fail;
-    
-    arg4 = (float*) array4->data;
-  }
-  {
-    temp5 = obj_to_array_no_conversion(obj4,PyArray_FLOAT);
-    if (!temp5  || !require_contiguous(temp5) || !require_native(temp5)) SWIG_fail;
-    arg5 = (float*) array_data(temp5);
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array6 = obj_to_array_contiguous_allow_conversion(obj5, PyArray_FLOAT, &is_new_object6);
-    if (!array6 || !require_dimensions(array6,1) || !require_size(array6,size,1)
-      || !require_contiguous(array6)   || !require_native(array6)) SWIG_fail;
-    
-    arg6 = (float*) array6->data;
-  }
-  ecode7 = SWIG_AsVal_long_SS_long(obj6, &val7);
-  if (!SWIG_IsOK(ecode7)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "gauss_seidel" "', argument " "7"" of type '" "long long""'");
-  } 
-  arg7 = static_cast< long long >(val7);
-  ecode8 = SWIG_AsVal_long_SS_long(obj7, &val8);
-  if (!SWIG_IsOK(ecode8)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "gauss_seidel" "', argument " "8"" of type '" "long long""'");
-  } 
-  arg8 = static_cast< long long >(val8);
-  ecode9 = SWIG_AsVal_long_SS_long(obj8, &val9);
-  if (!SWIG_IsOK(ecode9)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "gauss_seidel" "', argument " "9"" of type '" "long long""'");
-  } 
-  arg9 = static_cast< long long >(val9);
-  gauss_seidel<long long,float >(arg1,(long long const (*))arg2,(long long const (*))arg3,(float const (*))arg4,arg5,(float const (*))arg6,arg7,arg8,arg9);
-  resultobj = SWIG_Py_Void();
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return resultobj;
-fail:
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_gauss_seidel__SWIG_4(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  long long arg1 ;
-  long long *arg2 ;
-  long long *arg3 ;
-  double *arg4 ;
-  double *arg5 ;
-  double *arg6 ;
-  long long arg7 ;
-  long long arg8 ;
-  long long arg9 ;
-  long long val1 ;
-  int ecode1 = 0 ;
-  PyArrayObject *array2 = NULL ;
-  int is_new_object2 ;
-  PyArrayObject *array3 = NULL ;
-  int is_new_object3 ;
-  PyArrayObject *array4 = NULL ;
-  int is_new_object4 ;
-  PyArrayObject *temp5 = NULL ;
-  PyArrayObject *array6 = NULL ;
-  int is_new_object6 ;
-  long long val7 ;
-  int ecode7 = 0 ;
-  long long val8 ;
-  int ecode8 = 0 ;
-  long long val9 ;
-  int ecode9 = 0 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  PyObject * obj2 = 0 ;
-  PyObject * obj3 = 0 ;
-  PyObject * obj4 = 0 ;
-  PyObject * obj5 = 0 ;
-  PyObject * obj6 = 0 ;
-  PyObject * obj7 = 0 ;
-  PyObject * obj8 = 0 ;
-  
-  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOO:gauss_seidel",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8)) SWIG_fail;
-  ecode1 = SWIG_AsVal_long_SS_long(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "gauss_seidel" "', argument " "1"" of type '" "long long""'");
-  } 
-  arg1 = static_cast< long long >(val1);
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_LONGLONG, &is_new_object2);
-    if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
-      || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
-    
-    arg2 = (long long*) array2->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_LONGLONG, &is_new_object3);
-    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
-      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
-    
-    arg3 = (long long*) array3->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array4 = obj_to_array_contiguous_allow_conversion(obj3, PyArray_DOUBLE, &is_new_object4);
-    if (!array4 || !require_dimensions(array4,1) || !require_size(array4,size,1)
-      || !require_contiguous(array4)   || !require_native(array4)) SWIG_fail;
-    
-    arg4 = (double*) array4->data;
-  }
-  {
-    temp5 = obj_to_array_no_conversion(obj4,PyArray_DOUBLE);
-    if (!temp5  || !require_contiguous(temp5) || !require_native(temp5)) SWIG_fail;
-    arg5 = (double*) array_data(temp5);
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array6 = obj_to_array_contiguous_allow_conversion(obj5, PyArray_DOUBLE, &is_new_object6);
-    if (!array6 || !require_dimensions(array6,1) || !require_size(array6,size,1)
-      || !require_contiguous(array6)   || !require_native(array6)) SWIG_fail;
-    
-    arg6 = (double*) array6->data;
-  }
-  ecode7 = SWIG_AsVal_long_SS_long(obj6, &val7);
-  if (!SWIG_IsOK(ecode7)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "gauss_seidel" "', argument " "7"" of type '" "long long""'");
-  } 
-  arg7 = static_cast< long long >(val7);
-  ecode8 = SWIG_AsVal_long_SS_long(obj7, &val8);
-  if (!SWIG_IsOK(ecode8)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "gauss_seidel" "', argument " "8"" of type '" "long long""'");
-  } 
-  arg8 = static_cast< long long >(val8);
-  ecode9 = SWIG_AsVal_long_SS_long(obj8, &val9);
-  if (!SWIG_IsOK(ecode9)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "gauss_seidel" "', argument " "9"" of type '" "long long""'");
-  } 
-  arg9 = static_cast< long long >(val9);
-  gauss_seidel<long long,double >(arg1,(long long const (*))arg2,(long long const (*))arg3,(double const (*))arg4,arg5,(double const (*))arg6,arg7,arg8,arg9);
-  resultobj = SWIG_Py_Void();
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return resultobj;
-fail:
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_gauss_seidel(PyObject *self, PyObject *args) {
   int argc;
   PyObject *argv[10];
@@ -5641,115 +4845,9 @@
       }
     }
   }
-  if (argc == 9) {
-    int _v;
-    {
-      int res = SWIG_AsVal_long_SS_long(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      {
-        _v = (is_array(argv[1]) && PyArray_CanCastSafely(PyArray_TYPE(argv[1]),PyArray_LONGLONG)) ? 1 : 0;
-      }
-      if (_v) {
-        {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_LONGLONG)) ? 1 : 0;
-        }
-        if (_v) {
-          {
-            _v = (is_array(argv[3]) && PyArray_CanCastSafely(PyArray_TYPE(argv[3]),PyArray_FLOAT)) ? 1 : 0;
-          }
-          if (_v) {
-            {
-              _v = (is_array(argv[4]) && PyArray_CanCastSafely(PyArray_TYPE(argv[4]),PyArray_FLOAT)) ? 1 : 0;
-            }
-            if (_v) {
-              {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_FLOAT)) ? 1 : 0;
-              }
-              if (_v) {
-                {
-                  int res = SWIG_AsVal_long_SS_long(argv[6], NULL);
-                  _v = SWIG_CheckState(res);
-                }
-                if (_v) {
-                  {
-                    int res = SWIG_AsVal_long_SS_long(argv[7], NULL);
-                    _v = SWIG_CheckState(res);
-                  }
-                  if (_v) {
-                    {
-                      int res = SWIG_AsVal_long_SS_long(argv[8], NULL);
-                      _v = SWIG_CheckState(res);
-                    }
-                    if (_v) {
-                      return _wrap_gauss_seidel__SWIG_3(self, args);
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-  if (argc == 9) {
-    int _v;
-    {
-      int res = SWIG_AsVal_long_SS_long(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      {
-        _v = (is_array(argv[1]) && PyArray_CanCastSafely(PyArray_TYPE(argv[1]),PyArray_LONGLONG)) ? 1 : 0;
-      }
-      if (_v) {
-        {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_LONGLONG)) ? 1 : 0;
-        }
-        if (_v) {
-          {
-            _v = (is_array(argv[3]) && PyArray_CanCastSafely(PyArray_TYPE(argv[3]),PyArray_DOUBLE)) ? 1 : 0;
-          }
-          if (_v) {
-            {
-              _v = (is_array(argv[4]) && PyArray_CanCastSafely(PyArray_TYPE(argv[4]),PyArray_DOUBLE)) ? 1 : 0;
-            }
-            if (_v) {
-              {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_DOUBLE)) ? 1 : 0;
-              }
-              if (_v) {
-                {
-                  int res = SWIG_AsVal_long_SS_long(argv[6], NULL);
-                  _v = SWIG_CheckState(res);
-                }
-                if (_v) {
-                  {
-                    int res = SWIG_AsVal_long_SS_long(argv[7], NULL);
-                    _v = SWIG_CheckState(res);
-                  }
-                  if (_v) {
-                    {
-                      int res = SWIG_AsVal_long_SS_long(argv[8], NULL);
-                      _v = SWIG_CheckState(res);
-                    }
-                    if (_v) {
-                      return _wrap_gauss_seidel__SWIG_4(self, args);
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-    }
-  }
   
 fail:
-  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'gauss_seidel'.\n  Possible C/C++ prototypes are:\n""    gauss_seidel<(int,float)>(int const,int const [],int const [],float const [],float [],float const [],int const,int const,int const)\n""    gauss_seidel<(int,double)>(int const,int const [],int const [],double const [],double [],double const [],int const,int const,int const)\n""    gauss_seidel<(long long,float)>(long long const,long long const [],long long const [],float const [],float [],float const [],long long const,long long const,long long const)\n""    gauss_seidel<(long long,double)>(long long const,long long const [],long long const [],double const [],double [],double const [],long long const,long long const,long long const)\n");
+  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'gauss_seidel'.\n  Possible C/C++ prototypes are:\n""    gauss_seidel<(int,float)>(int const,int const [],int const [],float const [],float [],float const [],int const,int const,int const)\n""    gauss_seidel<(int,double)>(int const,int const [],int const [],double const [],double [],double const [],int const,int const,int const)\n");
   return NULL;
 }
 
@@ -6060,312 +5158,6 @@
 }
 
 
-SWIGINTERN PyObject *_wrap_jacobi__SWIG_3(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  long long arg1 ;
-  long long *arg2 ;
-  long long *arg3 ;
-  float *arg4 ;
-  float *arg5 ;
-  float *arg6 ;
-  float *arg7 ;
-  long long arg8 ;
-  long long arg9 ;
-  long long arg10 ;
-  float arg11 ;
-  long long val1 ;
-  int ecode1 = 0 ;
-  PyArrayObject *array2 = NULL ;
-  int is_new_object2 ;
-  PyArrayObject *array3 = NULL ;
-  int is_new_object3 ;
-  PyArrayObject *array4 = NULL ;
-  int is_new_object4 ;
-  PyArrayObject *temp5 = NULL ;
-  PyArrayObject *array6 = NULL ;
-  int is_new_object6 ;
-  PyArrayObject *temp7 = NULL ;
-  long long val8 ;
-  int ecode8 = 0 ;
-  long long val9 ;
-  int ecode9 = 0 ;
-  long long val10 ;
-  int ecode10 = 0 ;
-  float val11 ;
-  int ecode11 = 0 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  PyObject * obj2 = 0 ;
-  PyObject * obj3 = 0 ;
-  PyObject * obj4 = 0 ;
-  PyObject * obj5 = 0 ;
-  PyObject * obj6 = 0 ;
-  PyObject * obj7 = 0 ;
-  PyObject * obj8 = 0 ;
-  PyObject * obj9 = 0 ;
-  PyObject * obj10 = 0 ;
-  
-  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOOOO:jacobi",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8,&obj9,&obj10)) SWIG_fail;
-  ecode1 = SWIG_AsVal_long_SS_long(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "jacobi" "', argument " "1"" of type '" "long long""'");
-  } 
-  arg1 = static_cast< long long >(val1);
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_LONGLONG, &is_new_object2);
-    if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
-      || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
-    
-    arg2 = (long long*) array2->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_LONGLONG, &is_new_object3);
-    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
-      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
-    
-    arg3 = (long long*) array3->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array4 = obj_to_array_contiguous_allow_conversion(obj3, PyArray_FLOAT, &is_new_object4);
-    if (!array4 || !require_dimensions(array4,1) || !require_size(array4,size,1)
-      || !require_contiguous(array4)   || !require_native(array4)) SWIG_fail;
-    
-    arg4 = (float*) array4->data;
-  }
-  {
-    temp5 = obj_to_array_no_conversion(obj4,PyArray_FLOAT);
-    if (!temp5  || !require_contiguous(temp5) || !require_native(temp5)) SWIG_fail;
-    arg5 = (float*) array_data(temp5);
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array6 = obj_to_array_contiguous_allow_conversion(obj5, PyArray_FLOAT, &is_new_object6);
-    if (!array6 || !require_dimensions(array6,1) || !require_size(array6,size,1)
-      || !require_contiguous(array6)   || !require_native(array6)) SWIG_fail;
-    
-    arg6 = (float*) array6->data;
-  }
-  {
-    temp7 = obj_to_array_no_conversion(obj6,PyArray_FLOAT);
-    if (!temp7  || !require_contiguous(temp7) || !require_native(temp7)) SWIG_fail;
-    arg7 = (float*) array_data(temp7);
-  }
-  ecode8 = SWIG_AsVal_long_SS_long(obj7, &val8);
-  if (!SWIG_IsOK(ecode8)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "jacobi" "', argument " "8"" of type '" "long long""'");
-  } 
-  arg8 = static_cast< long long >(val8);
-  ecode9 = SWIG_AsVal_long_SS_long(obj8, &val9);
-  if (!SWIG_IsOK(ecode9)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "jacobi" "', argument " "9"" of type '" "long long""'");
-  } 
-  arg9 = static_cast< long long >(val9);
-  ecode10 = SWIG_AsVal_long_SS_long(obj9, &val10);
-  if (!SWIG_IsOK(ecode10)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode10), "in method '" "jacobi" "', argument " "10"" of type '" "long long""'");
-  } 
-  arg10 = static_cast< long long >(val10);
-  ecode11 = SWIG_AsVal_float(obj10, &val11);
-  if (!SWIG_IsOK(ecode11)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode11), "in method '" "jacobi" "', argument " "11"" of type '" "float""'");
-  } 
-  arg11 = static_cast< float >(val11);
-  jacobi<long long,float >(arg1,(long long const (*))arg2,(long long const (*))arg3,(float const (*))arg4,arg5,(float const (*))arg6,arg7,arg8,arg9,arg10,arg11);
-  resultobj = SWIG_Py_Void();
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return resultobj;
-fail:
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_jacobi__SWIG_4(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  long long arg1 ;
-  long long *arg2 ;
-  long long *arg3 ;
-  double *arg4 ;
-  double *arg5 ;
-  double *arg6 ;
-  double *arg7 ;
-  long long arg8 ;
-  long long arg9 ;
-  long long arg10 ;
-  double arg11 ;
-  long long val1 ;
-  int ecode1 = 0 ;
-  PyArrayObject *array2 = NULL ;
-  int is_new_object2 ;
-  PyArrayObject *array3 = NULL ;
-  int is_new_object3 ;
-  PyArrayObject *array4 = NULL ;
-  int is_new_object4 ;
-  PyArrayObject *temp5 = NULL ;
-  PyArrayObject *array6 = NULL ;
-  int is_new_object6 ;
-  PyArrayObject *temp7 = NULL ;
-  long long val8 ;
-  int ecode8 = 0 ;
-  long long val9 ;
-  int ecode9 = 0 ;
-  long long val10 ;
-  int ecode10 = 0 ;
-  double val11 ;
-  int ecode11 = 0 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  PyObject * obj2 = 0 ;
-  PyObject * obj3 = 0 ;
-  PyObject * obj4 = 0 ;
-  PyObject * obj5 = 0 ;
-  PyObject * obj6 = 0 ;
-  PyObject * obj7 = 0 ;
-  PyObject * obj8 = 0 ;
-  PyObject * obj9 = 0 ;
-  PyObject * obj10 = 0 ;
-  
-  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOOOO:jacobi",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8,&obj9,&obj10)) SWIG_fail;
-  ecode1 = SWIG_AsVal_long_SS_long(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "jacobi" "', argument " "1"" of type '" "long long""'");
-  } 
-  arg1 = static_cast< long long >(val1);
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_LONGLONG, &is_new_object2);
-    if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
-      || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
-    
-    arg2 = (long long*) array2->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_LONGLONG, &is_new_object3);
-    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
-      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
-    
-    arg3 = (long long*) array3->data;
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array4 = obj_to_array_contiguous_allow_conversion(obj3, PyArray_DOUBLE, &is_new_object4);
-    if (!array4 || !require_dimensions(array4,1) || !require_size(array4,size,1)
-      || !require_contiguous(array4)   || !require_native(array4)) SWIG_fail;
-    
-    arg4 = (double*) array4->data;
-  }
-  {
-    temp5 = obj_to_array_no_conversion(obj4,PyArray_DOUBLE);
-    if (!temp5  || !require_contiguous(temp5) || !require_native(temp5)) SWIG_fail;
-    arg5 = (double*) array_data(temp5);
-  }
-  {
-    npy_intp size[1] = {
-      -1
-    };
-    array6 = obj_to_array_contiguous_allow_conversion(obj5, PyArray_DOUBLE, &is_new_object6);
-    if (!array6 || !require_dimensions(array6,1) || !require_size(array6,size,1)
-      || !require_contiguous(array6)   || !require_native(array6)) SWIG_fail;
-    
-    arg6 = (double*) array6->data;
-  }
-  {
-    temp7 = obj_to_array_no_conversion(obj6,PyArray_DOUBLE);
-    if (!temp7  || !require_contiguous(temp7) || !require_native(temp7)) SWIG_fail;
-    arg7 = (double*) array_data(temp7);
-  }
-  ecode8 = SWIG_AsVal_long_SS_long(obj7, &val8);
-  if (!SWIG_IsOK(ecode8)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "jacobi" "', argument " "8"" of type '" "long long""'");
-  } 
-  arg8 = static_cast< long long >(val8);
-  ecode9 = SWIG_AsVal_long_SS_long(obj8, &val9);
-  if (!SWIG_IsOK(ecode9)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "jacobi" "', argument " "9"" of type '" "long long""'");
-  } 
-  arg9 = static_cast< long long >(val9);
-  ecode10 = SWIG_AsVal_long_SS_long(obj9, &val10);
-  if (!SWIG_IsOK(ecode10)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode10), "in method '" "jacobi" "', argument " "10"" of type '" "long long""'");
-  } 
-  arg10 = static_cast< long long >(val10);
-  ecode11 = SWIG_AsVal_double(obj10, &val11);
-  if (!SWIG_IsOK(ecode11)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode11), "in method '" "jacobi" "', argument " "11"" of type '" "double""'");
-  } 
-  arg11 = static_cast< double >(val11);
-  jacobi<long long,double >(arg1,(long long const (*))arg2,(long long const (*))arg3,(double const (*))arg4,arg5,(double const (*))arg6,arg7,arg8,arg9,arg10,arg11);
-  resultobj = SWIG_Py_Void();
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return resultobj;
-fail:
-  {
-    if (is_new_object2 && array2) Py_DECREF(array2);
-  }
-  {
-    if (is_new_object3 && array3) Py_DECREF(array3);
-  }
-  {
-    if (is_new_object4 && array4) Py_DECREF(array4);
-  }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_jacobi(PyObject *self, PyObject *args) {
   int argc;
   PyObject *argv[12];
@@ -6504,137 +5296,9 @@
       }
     }
   }
-  if (argc == 11) {
-    int _v;
-    {
-      int res = SWIG_AsVal_long_SS_long(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      {
-        _v = (is_array(argv[1]) && PyArray_CanCastSafely(PyArray_TYPE(argv[1]),PyArray_LONGLONG)) ? 1 : 0;
-      }
-      if (_v) {
-        {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_LONGLONG)) ? 1 : 0;
-        }
-        if (_v) {
-          {
-            _v = (is_array(argv[3]) && PyArray_CanCastSafely(PyArray_TYPE(argv[3]),PyArray_FLOAT)) ? 1 : 0;
-          }
-          if (_v) {
-            {
-              _v = (is_array(argv[4]) && PyArray_CanCastSafely(PyArray_TYPE(argv[4]),PyArray_FLOAT)) ? 1 : 0;
-            }
-            if (_v) {
-              {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_FLOAT)) ? 1 : 0;
-              }
-              if (_v) {
-                {
-                  _v = (is_array(argv[6]) && PyArray_CanCastSafely(PyArray_TYPE(argv[6]),PyArray_FLOAT)) ? 1 : 0;
-                }
-                if (_v) {
-                  {
-                    int res = SWIG_AsVal_long_SS_long(argv[7], NULL);
-                    _v = SWIG_CheckState(res);
-                  }
-                  if (_v) {
-                    {
-                      int res = SWIG_AsVal_long_SS_long(argv[8], NULL);
-                      _v = SWIG_CheckState(res);
-                    }
-                    if (_v) {
-                      {
-                        int res = SWIG_AsVal_long_SS_long(argv[9], NULL);
-                        _v = SWIG_CheckState(res);
-                      }
-                      if (_v) {
-                        {
-                          int res = SWIG_AsVal_float(argv[10], NULL);
-                          _v = SWIG_CheckState(res);
-                        }
-                        if (_v) {
-                          return _wrap_jacobi__SWIG_3(self, args);
-                        }
-                      }
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-  if (argc == 11) {
-    int _v;
-    {
-      int res = SWIG_AsVal_long_SS_long(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      {
-        _v = (is_array(argv[1]) && PyArray_CanCastSafely(PyArray_TYPE(argv[1]),PyArray_LONGLONG)) ? 1 : 0;
-      }
-      if (_v) {
-        {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_LONGLONG)) ? 1 : 0;
-        }
-        if (_v) {
-          {
-            _v = (is_array(argv[3]) && PyArray_CanCastSafely(PyArray_TYPE(argv[3]),PyArray_DOUBLE)) ? 1 : 0;
-          }
-          if (_v) {
-            {
-              _v = (is_array(argv[4]) && PyArray_CanCastSafely(PyArray_TYPE(argv[4]),PyArray_DOUBLE)) ? 1 : 0;
-            }
-            if (_v) {
-              {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_DOUBLE)) ? 1 : 0;
-              }
-              if (_v) {
-                {
-                  _v = (is_array(argv[6]) && PyArray_CanCastSafely(PyArray_TYPE(argv[6]),PyArray_DOUBLE)) ? 1 : 0;
-                }
-                if (_v) {
-                  {
-                    int res = SWIG_AsVal_long_SS_long(argv[7], NULL);
-                    _v = SWIG_CheckState(res);
-                  }
-                  if (_v) {
-                    {
-                      int res = SWIG_AsVal_long_SS_long(argv[8], NULL);
-                      _v = SWIG_CheckState(res);
-                    }
-                    if (_v) {
-                      {
-                        int res = SWIG_AsVal_long_SS_long(argv[9], NULL);
-                        _v = SWIG_CheckState(res);
-                      }
-                      if (_v) {
-                        {
-                          int res = SWIG_AsVal_double(argv[10], NULL);
-                          _v = SWIG_CheckState(res);
-                        }
-                        if (_v) {
-                          return _wrap_jacobi__SWIG_4(self, args);
-                        }
-                      }
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-    }
-  }
   
 fail:
-  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'jacobi'.\n  Possible C/C++ prototypes are:\n""    jacobi<(int,float)>(int const,int const [],int const [],float const [],float [],float const [],float [],int const,int const,int const,float const)\n""    jacobi<(int,double)>(int const,int const [],int const [],double const [],double [],double const [],double [],int const,int const,int const,double const)\n""    jacobi<(long long,float)>(long long const,long long const [],long long const [],float const [],float [],float const [],float [],long long const,long long const,long long const,float const)\n""    jacobi<(long long,double)>(long long const,long long const [],long long const [],double const [],double [],double const [],double [],long long const,long long const,long long const,double const)\n");
+  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'jacobi'.\n  Possible C/C++ prototypes are:\n""    jacobi<(int,float)>(int const,int const [],int const [],float const [],float [],float const [],float [],int const,int const,int const,float const)\n""    jacobi<(int,double)>(int const,int const [],int const [],double const [],double [],double const [],double [],int const,int const,int const,double const)\n");
   return NULL;
 }
 
@@ -6665,25 +5329,11 @@
 		"    std::vector<(int)> Sp, std::vector<(int)> Sj, \n"
 		"    std::vector<(double)> Sx)\n"
 		""},
-	 { (char *)"sa_smoother", _wrap_sa_smoother, METH_VARARGS, (char *)"\n"
-		"sa_smoother(int n_row, float omega, int Ap, int Aj, float Ax, int Sp, \n"
-		"    int Sj, float Sx, std::vector<(int)> Bp, \n"
-		"    std::vector<(int)> Bj, std::vector<(float)> Bx)\n"
-		"sa_smoother(int n_row, double omega, int Ap, int Aj, double Ax, \n"
-		"    int Sp, int Sj, double Sx, std::vector<(int)> Bp, \n"
-		"    std::vector<(int)> Bj, std::vector<(double)> Bx)\n"
-		""},
 	 { (char *)"gauss_seidel", _wrap_gauss_seidel, METH_VARARGS, (char *)"\n"
 		"gauss_seidel(int n_row, int Ap, int Aj, float Ax, float x, float b, \n"
 		"    int row_start, int row_stop, int row_step)\n"
 		"gauss_seidel(int n_row, int Ap, int Aj, double Ax, double x, double b, \n"
 		"    int row_start, int row_stop, int row_step)\n"
-		"gauss_seidel(long long n_row, long long Ap, long long Aj, float Ax, \n"
-		"    float x, float b, long long row_start, long long row_stop, \n"
-		"    long long row_step)\n"
-		"gauss_seidel(long long n_row, long long Ap, long long Aj, double Ax, \n"
-		"    double x, double b, long long row_start, \n"
-		"    long long row_stop, long long row_step)\n"
 		""},
 	 { (char *)"jacobi", _wrap_jacobi, METH_VARARGS, (char *)"\n"
 		"jacobi(int n_row, int Ap, int Aj, float Ax, float x, float b, \n"
@@ -6692,14 +5342,6 @@
 		"jacobi(int n_row, int Ap, int Aj, double Ax, double x, double b, \n"
 		"    double temp, int row_start, int row_stop, \n"
 		"    int row_step, double omega)\n"
-		"jacobi(long long n_row, long long Ap, long long Aj, float Ax, \n"
-		"    float x, float b, float temp, long long row_start, \n"
-		"    long long row_stop, long long row_step, \n"
-		"    float omega)\n"
-		"jacobi(long long n_row, long long Ap, long long Aj, double Ax, \n"
-		"    double x, double b, double temp, long long row_start, \n"
-		"    long long row_stop, long long row_step, \n"
-		"    double omega)\n"
 		""},
 	 { NULL, NULL, 0, NULL }
 };

Modified: trunk/scipy/sandbox/multigrid/multigridtools/ruge_stuben.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/ruge_stuben.h	2007-12-03 22:01:43 UTC (rev 3612)
+++ trunk/scipy/sandbox/multigrid/multigridtools/ruge_stuben.h	2007-12-03 23:30:01 UTC (rev 3613)
@@ -13,30 +13,31 @@
 enum NodeType {U_NODE,  C_NODE, F_NODE};
 
 
-template<class T>
-void rs_strong_connections(const int n_row,
-			   const T theta,
-			   const int Ap[], const int Aj[], const T Ax[],
-			   std::vector<int> * Sp, std::vector<int> * Sj, std::vector<T> * Sx){
+template<class I, class T>
+void rs_strong_connections(const I n_row,
+                           const T theta,
+                           const I Ap[], const I Aj[], const T Ax[],
+                           std::vector<I> * Sp, std::vector<I> * Sj, std::vector<T> * Sx){
+
     //Sp,Sj form a CSR representation where the i-th row contains
     //the indices of all the strong connections from node i
     Sp->push_back(0);
 
     //Compute lambdas for each node
-    for(int i = 0; i < n_row; i++){
+    for(I i = 0; i < n_row; i++){
         T min_offdiagonal = 0.0;
 
-        int row_start = Ap[i];
-        int row_end   = Ap[i+1];
-        for(int jj = row_start; jj < row_end; jj++){
+        I row_start = Ap[i];
+        I row_end   = Ap[i+1];
+        for(I jj = row_start; jj < row_end; jj++){
             min_offdiagonal = std::min(min_offdiagonal,Ax[jj]); //assumes diagonal is positive!
         }
 
         T threshold = theta*min_offdiagonal;
-        for(int jj = row_start; jj < row_end; jj++){
+        for(I jj = row_start; jj < row_end; jj++){
             if(Ax[jj] < threshold){
-	            Sj->push_back(Aj[jj]);
-	            Sx->push_back(Ax[jj]);
+                Sj->push_back(Aj[jj]);
+                Sx->push_back(Ax[jj]);
             }
         }
 
@@ -47,341 +48,339 @@
 
 
 
-template<class T>
-void rs_interpolation(const int n_nodes,
-		      const int Ap[], const int Aj[], const T Ax[],
-		      const int Sp[], const int Sj[], const T Sx[],
-		      const int Tp[], const int Tj[], const T Tx[],
-		      std::vector<int> * Bp, std::vector<int> * Bj, std::vector<T> * Bx){
-  
-  std::vector<int> lambda(n_nodes,0);
+template<class I, class T>
+void rs_interpolation(const I n_nodes,
+        const I Ap[], const I Aj[], const T Ax[],
+        const I Sp[], const I Sj[], const T Sx[],
+        const I Tp[], const I Tj[], const T Tx[],
+        std::vector<I> * Bp, std::vector<I> * Bj, std::vector<T> * Bx){
 
-  //compute lambdas
-  for(int i = 0; i < n_nodes; i++){
-    lambda[i] = Tp[i+1] - Tp[i];
-  }
+    std::vector<I> lambda(n_nodes,0);
 
+    //compute lambdas
+    for(I i = 0; i < n_nodes; i++){
+        lambda[i] = Tp[i+1] - Tp[i];
+    }
 
-  //for each value of lambda, create an interval of nodes with that value
-  // ptr - is the first index of the interval
-  // count - is the number of indices in that interval
-  // index to node - the node located at a given index
-  // node to index - the index of a given node
-  std::vector<int> interval_ptr(n_nodes,0);
-  std::vector<int> interval_count(n_nodes,0);
-  std::vector<int> index_to_node(n_nodes);
-  std::vector<int> node_to_index(n_nodes);
 
-  for(int i = 0; i < n_nodes; i++){
-    interval_count[lambda[i]]++;
-  }
-  for(int i = 0, cumsum = 0; i < n_nodes; i++){
-    interval_ptr[i] = cumsum;
-    cumsum += interval_count[i];
-    interval_count[i] = 0;
-  }
-  for(int i = 0; i < n_nodes; i++){
-    int lambda_i = lambda[i];
-    int index    = interval_ptr[lambda_i]+interval_count[lambda_i];
-    index_to_node[index] = i;
-    node_to_index[i]     = index;
-    interval_count[lambda_i]++;
-  }
+    //for each value of lambda, create an interval of nodes with that value
+    // ptr - is the first index of the interval
+    // count - is the number of indices in that interval
+    // index to node - the node located at a given index
+    // node to index - the index of a given node
+    std::vector<I> interval_ptr(n_nodes,0);
+    std::vector<I> interval_count(n_nodes,0);
+    std::vector<I> index_to_node(n_nodes);
+    std::vector<I> node_to_index(n_nodes);
 
+    for(I i = 0; i < n_nodes; i++){
+        interval_count[lambda[i]]++;
+    }
+    for(I i = 0, cumsum = 0; i < n_nodes; i++){
+        interval_ptr[i] = cumsum;
+        cumsum += interval_count[i];
+        interval_count[i] = 0;
+    }
+    for(I i = 0; i < n_nodes; i++){
+        I lambda_i = lambda[i];
+        I index    = interval_ptr[lambda_i]+interval_count[lambda_i];
+        index_to_node[index] = i;
+        node_to_index[i]     = index;
+        interval_count[lambda_i]++;
+    }
 
 
-  
 
-  std::vector<NodeType> NodeSets(n_nodes,U_NODE);
 
-  //Now add elements to C and F, in decending order of lambda
-  for(int top_index = n_nodes - 1; top_index > -1; top_index--){
-    int i        = index_to_node[top_index];
-    int lambda_i = lambda[i];
+
+    std::vector<NodeType> NodeSets(n_nodes,U_NODE);
+
+    //Now add elements to C and F, in decending order of lambda
+    for(I top_index = n_nodes - 1; top_index > -1; top_index--){
+        I i        = index_to_node[top_index];
+        I lambda_i = lambda[i];
 #ifdef DEBUG
-    {
+        {
 #ifdef DEBUG_PRINT
-      std::cout << "top_index " << top_index << std::endl;
-      std::cout << "i         " << i << std::endl;
-      std::cout << "lambda_i  " << lambda_i << std::endl;
+            std::cout << "top_index " << top_index << std::endl;
+            std::cout << "i         " << i << std::endl;
+            std::cout << "lambda_i  " << lambda_i << std::endl;
 
-      for(int i = 0; i < n_nodes; i++){
-	std::cout << i << "=";
-	if(NodeSets[i] == U_NODE)
-	  std::cout << "U";
-	else if(NodeSets[i] == F_NODE)
-	  std::cout << "F";
-	else
-	  std::cout << "C";
-	std::cout << " ";
-      }
-      std::cout << std::endl;
+            for(I i = 0; i < n_nodes; i++){
+                std::cout << i << "=";
+                if(NodeSets[i] == U_NODE)
+                    std::cout << "U";
+                else if(NodeSets[i] == F_NODE)
+                    std::cout << "F";
+                else
+                    std::cout << "C";
+                std::cout << " ";
+            }
+            std::cout << std::endl;
 
-      std::cout << "node_to_index" << std::endl;
-      for(int i = 0; i < n_nodes; i++){
-	std::cout << i << "->" << node_to_index[i] << "  ";
-      }
-      std::cout << std::endl;
-      std::cout << "index_to_node" << std::endl;
-      for(int i = 0; i < n_nodes; i++){
-	std::cout << i << "->" << index_to_node[i] << "  ";
-      }
-      std::cout << std::endl;
+            std::cout << "node_to_index" << std::endl;
+            for(I i = 0; i < n_nodes; i++){
+                std::cout << i << "->" << node_to_index[i] << "  ";
+            }
+            std::cout << std::endl;
+            std::cout << "index_to_node" << std::endl;
+            for(I i = 0; i < n_nodes; i++){
+                std::cout << i << "->" << index_to_node[i] << "  ";
+            }
+            std::cout << std::endl;
 
-      std::cout << "interval_count ";
-      for(int i = 0; i < n_nodes; i++){
-	std::cout << interval_count[i] << " ";
-      }
-      std::cout << std::endl;
+            std::cout << "interval_count ";
+            for(I i = 0; i < n_nodes; i++){
+                std::cout << interval_count[i] << " ";
+            }
+            std::cout << std::endl;
 #endif
 
-      //make sure arrays are correct
-      for(int n = 0; n < n_nodes; n++){
-	assert(index_to_node[node_to_index[n]] == n);
-      }
+            //make sure arrays are correct
+            for(I n = 0; n < n_nodes; n++){
+                assert(index_to_node[node_to_index[n]] == n);
+            }
 
-      //make sure intervals are reasonable
-      int sum_intervals = 0;
-      for(int n = 0; n < n_nodes; n++){
-	assert(interval_count[n] >= 0);
-	if(interval_count[n] > 0){
-	  assert(interval_ptr[n] == sum_intervals);
-	}
-	sum_intervals += interval_count[n];
-      }
-      assert(sum_intervals == top_index+1);
+            //make sure intervals are reasonable
+            I sum_intervals = 0;
+            for(I n = 0; n < n_nodes; n++){
+                assert(interval_count[n] >= 0);
+                if(interval_count[n] > 0){
+                    assert(interval_ptr[n] == sum_intervals);
+                }
+                sum_intervals += interval_count[n];
+            }
+            assert(sum_intervals == top_index+1);
 
-      
-      if(interval_count[lambda_i] <= 0){
-	std::cout << "top_index " << top_index << std::endl;
-	std::cout << "lambda_i " << lambda_i << std::endl;
-	std::cout << "interval_count[lambda_i] " << interval_count[lambda_i] << std::endl;
-	std::cout << "top_index " << top_index << std::endl;
-	std::cout << "i         " << i << std::endl;
-	std::cout << "lambda_i  " << lambda_i << std::endl;
-      }
-      
-      
-      for(int n = 0; n <= top_index; n++){
-	assert(NodeSets[index_to_node[n]] != C_NODE);
-      }
-    }
-    assert(node_to_index[i] == top_index);
-    assert(interval_ptr[lambda_i] + interval_count[lambda_i] - 1 == top_index);
-    //max interval should have at least one element
-    assert(interval_count[lambda_i] > 0);    
+
+            if(interval_count[lambda_i] <= 0){
+                std::cout << "top_index " << top_index << std::endl;
+                std::cout << "lambda_i " << lambda_i << std::endl;
+                std::cout << "interval_count[lambda_i] " << interval_count[lambda_i] << std::endl;
+                std::cout << "top_index " << top_index << std::endl;
+                std::cout << "i         " << i << std::endl;
+                std::cout << "lambda_i  " << lambda_i << std::endl;
+            }
+
+
+            for(I n = 0; n <= top_index; n++){
+                assert(NodeSets[index_to_node[n]] != C_NODE);
+            }
+        }
+        assert(node_to_index[i] == top_index);
+        assert(interval_ptr[lambda_i] + interval_count[lambda_i] - 1 == top_index);
+        //max interval should have at least one element
+        assert(interval_count[lambda_i] > 0);    
 #endif
 
 
-    //remove i from its interval
-    interval_count[lambda_i]--;
-    
+        //remove i from its interval
+        interval_count[lambda_i]--;
 
-    if(NodeSets[i] == F_NODE){
-      continue;
-    } else {
-      assert(NodeSets[i] == U_NODE);
 
-      NodeSets[i] = C_NODE;
+        if(NodeSets[i] == F_NODE){
+            continue;
+        } else {
+            assert(NodeSets[i] == U_NODE);
 
-      //For each j in S^T_i /\ U
-      for(int jj = Tp[i]; jj < Tp[i+1]; jj++){
-	int j = Tj[jj];
+            NodeSets[i] = C_NODE;
 
-	if(NodeSets[j] == U_NODE){
-	  NodeSets[j] = F_NODE;
-	  
-	  //For each k in S_j /\ U
-	  for(int kk = Sp[j]; kk < Sp[j+1]; kk++){
-	    int k = Sj[kk];
+            //For each j in S^T_i /\ U
+            for(I jj = Tp[i]; jj < Tp[i+1]; jj++){
+                I j = Tj[jj];
 
-	    if(NodeSets[k] == U_NODE){	      
-	      //move k to the end of its current interval
-	      assert(lambda[j] < n_nodes - 1);//this would cause problems!
+                if(NodeSets[j] == U_NODE){
+                    NodeSets[j] = F_NODE;
 
-	      int lambda_k = lambda[k];
-	      int old_pos  = node_to_index[k];
-	      int new_pos  = interval_ptr[lambda_k] + interval_count[lambda_k] - 1;
+                    //For each k in S_j /\ U
+                    for(I kk = Sp[j]; kk < Sp[j+1]; kk++){
+                        I k = Sj[kk];
 
-	      node_to_index[index_to_node[old_pos]] = new_pos;
-	      node_to_index[index_to_node[new_pos]] = old_pos;
-	      std::swap(index_to_node[old_pos],index_to_node[new_pos]);
-	      
-	      //update intervals
-	      interval_count[lambda_k]   -= 1;
-	      interval_count[lambda_k+1] += 1;
-	      interval_ptr[lambda_k+1]    = new_pos;
+                        if(NodeSets[k] == U_NODE){	      
+                            //move k to the end of its current interval
+                            assert(lambda[j] < n_nodes - 1);//this would cause problems!
 
-	      //increment lambda_k
-	      lambda[k]++;
+                            I lambda_k = lambda[k];
+                            I old_pos  = node_to_index[k];
+                            I new_pos  = interval_ptr[lambda_k] + interval_count[lambda_k] - 1;
 
+                            node_to_index[index_to_node[old_pos]] = new_pos;
+                            node_to_index[index_to_node[new_pos]] = old_pos;
+                            std::swap(index_to_node[old_pos],index_to_node[new_pos]);
+
+                            //update intervals
+                            interval_count[lambda_k]   -= 1;
+                            interval_count[lambda_k+1] += 1;
+                            interval_ptr[lambda_k+1]    = new_pos;
+
+                            //increment lambda_k
+                            lambda[k]++;
+
 #ifdef DEBUG
-	      assert(interval_count[lambda_k]   >= 0);
-	      assert(interval_count[lambda_k+1] >  0);
-	      assert(interval_ptr[lambda[k]] <= node_to_index[k]);
-	      assert(node_to_index[k] < interval_ptr[lambda[k]] + interval_count[lambda[k]]);
+                            assert(interval_count[lambda_k]   >= 0);
+                            assert(interval_count[lambda_k+1] >  0);
+                            assert(interval_ptr[lambda[k]] <= node_to_index[k]);
+                            assert(node_to_index[k] < interval_ptr[lambda[k]] + interval_count[lambda[k]]);
 #endif
-	    }
-	  }
-	}
-      }
+                        }
+                    }
+                }
+            }
 
-      //For each j in S_i /\ U
-      for(int jj = Sp[i]; jj < Sp[i+1]; jj++){
-	int j = Sj[jj];
-	if(NodeSets[j] == U_NODE){            //decrement lambda for node j
-	  assert(lambda[j] > 0);//this would cause problems!
+            //For each j in S_i /\ U
+            for(I jj = Sp[i]; jj < Sp[i+1]; jj++){
+                I j = Sj[jj];
+                if(NodeSets[j] == U_NODE){            //decrement lambda for node j
+                    assert(lambda[j] > 0);//this would cause problems!
 
-	  //move j to the beginning of its current interval
-	  int lambda_j = lambda[j];
-	  int old_pos  = node_to_index[j];
-	  int new_pos  = interval_ptr[lambda_j]; 
-	      
-	  node_to_index[index_to_node[old_pos]] = new_pos;
-	  node_to_index[index_to_node[new_pos]] = old_pos;
-	  std::swap(index_to_node[old_pos],index_to_node[new_pos]);
-	      
-	  //update intervals
-	  interval_count[lambda_j]   -= 1;
-	  interval_count[lambda_j-1] += 1;
-	  interval_ptr[lambda_j]     += 1;
-	  interval_ptr[lambda_j-1]    = interval_ptr[lambda_j] - interval_count[lambda_j-1];
+                    //move j to the beginning of its current interval
+                    I lambda_j = lambda[j];
+                    I old_pos  = node_to_index[j];
+                    I new_pos  = interval_ptr[lambda_j]; 
 
-	  //decrement lambda_j
-	  lambda[j]--;
+                    node_to_index[index_to_node[old_pos]] = new_pos;
+                    node_to_index[index_to_node[new_pos]] = old_pos;
+                    std::swap(index_to_node[old_pos],index_to_node[new_pos]);
 
+                    //update intervals
+                    interval_count[lambda_j]   -= 1;
+                    interval_count[lambda_j-1] += 1;
+                    interval_ptr[lambda_j]     += 1;
+                    interval_ptr[lambda_j-1]    = interval_ptr[lambda_j] - interval_count[lambda_j-1];
+
+                    //decrement lambda_j
+                    lambda[j]--;
+
 #ifdef DEBUG
-	  assert(interval_count[lambda_j]   >= 0);
-	  assert(interval_count[lambda_j-1] >  0);
-	  assert(interval_ptr[lambda[j]] <= node_to_index[j]);
-	  assert(node_to_index[j] < interval_ptr[lambda[j]] + interval_count[lambda[j]]);
+                    assert(interval_count[lambda_j]   >= 0);
+                    assert(interval_count[lambda_j-1] >  0);
+                    assert(interval_ptr[lambda[j]] <= node_to_index[j]);
+                    assert(node_to_index[j] < interval_ptr[lambda[j]] + interval_count[lambda[j]]);
 #endif
-	}
-      }
+                }
+            }
+        }
     }
-  }
 
 
 
 
 #ifdef DEBUG
-  //make sure each f-node has at least one strong c-node neighbor
-  for(int i = 0; i < n_nodes; i++){
-    if(NodeSets[i] == F_NODE){
-      int row_start = Sp[i];
-      int row_end   = Sp[i+1];
-      bool has_c_neighbor = false;
-      for(int jj = row_start; jj < row_end; jj++){
-	if(NodeSets[Sj[jj]] == C_NODE){
-	  has_c_neighbor = true;
-	  break;
-	}
-      }
-      assert(has_c_neighbor);
-    }   
-  }
+    //make sure each f-node has at least one strong c-node neighbor
+    for(I i = 0; i < n_nodes; i++){
+        if(NodeSets[i] == F_NODE){
+            I row_start = Sp[i];
+            I row_end   = Sp[i+1];
+            bool has_c_neighbor = false;
+            for(I jj = row_start; jj < row_end; jj++){
+                if(NodeSets[Sj[jj]] == C_NODE){
+                    has_c_neighbor = true;
+                    break;
+                }
+            }
+            assert(has_c_neighbor);
+        }   
+    }
 #endif
 
-  //Now construct interpolation operator
-  std::vector<T> d_k(n_nodes,0);
-  std::vector<bool> C_i(n_nodes,0);
-  Bp->push_back(0);
-  for(int i = 0; i < n_nodes; i++){
-    if(NodeSets[i] == C_NODE){
-      //interpolate directly
-      Bj->push_back(i);
-      Bx->push_back(1);      
-      Bp->push_back(Bj->size());
-    } else {
-      //F_NODE
-      
-      //Step 4
-      T d_i = 0; //denominator for this row
-      for(int jj = Ap[i]; jj < Ap[i+1]; jj++){ d_i += Ax[jj]; }
-      for(int jj = Sp[i]; jj < Sp[i+1]; jj++){ d_i -= Sx[jj]; }
-      
-      //Create C_i, initialize d_k
-      for(int jj = Sp[i]; jj < Sp[i+1]; jj++){ 
-	int j = Sj[jj];
-	if(NodeSets[j] == C_NODE){
-	  C_i[j] = true;
-	  d_k[j] = Sx[jj];
-	}
-      }
+    //Now construct interpolation operator
+    std::vector<T> d_k(n_nodes,0);
+    std::vector<bool> C_i(n_nodes,0);
+    Bp->push_back(0);
+    for(I i = 0; i < n_nodes; i++){
+        if(NodeSets[i] == C_NODE){
+            //interpolate directly
+            Bj->push_back(i);
+            Bx->push_back(1);      
+            Bp->push_back(Bj->size());
+        } else {
+            //F_NODE
 
-      bool Sj_intersects_Ci = true; //in the case that i has no F-neighbors
-      for(int jj = Sp[i]; jj < Sp[i+1]; jj++){ //for j in D^s_i
-	int    j = Sj[jj];
-	T   a_ij = Sx[jj];
-	T   a_jl = 0;
+            //Step 4
+            T d_i = 0; //denominator for this row
+            for(I jj = Ap[i]; jj < Ap[i+1]; jj++){ d_i += Ax[jj]; }
+            for(I jj = Sp[i]; jj < Sp[i+1]; jj++){ d_i -= Sx[jj]; }
 
-	if(NodeSets[j] != F_NODE){continue;}
+            //Create C_i, initialize d_k
+            for(I jj = Sp[i]; jj < Sp[i+1]; jj++){ 
+                I j = Sj[jj];
+                if(NodeSets[j] == C_NODE){
+                    C_i[j] = true;
+                    d_k[j] = Sx[jj];
+                }
+            }
 
-	//Step 5
-	Sj_intersects_Ci = false;
+            bool Sj_intersects_Ci = true; //in the case that i has no F-neighbors
+            for(I jj = Sp[i]; jj < Sp[i+1]; jj++){ //for j in D^s_i
+                I    j = Sj[jj];
+                T   a_ij = Sx[jj];
+                T   a_jl = 0;
 
-	//compute sum a_jl
-	for(int ll = Sp[j]; ll < Sp[j+1]; ll++){
-	  if(C_i[Sj[ll]]){
-	    Sj_intersects_Ci = true;
-	    a_jl += Sx[ll];
-	  }	    
-	}
+                if(NodeSets[j] != F_NODE){continue;}
 
-	if(!Sj_intersects_Ci){ break; }
+                //Step 5
+                Sj_intersects_Ci = false;
 
-	for(int kk = Sp[j]; kk < Sp[j+1]; kk++){
-	  int   k = Sj[kk];
-	  T  a_jk = Sx[kk];
-	  if(C_i[k]){
-	    d_k[k] += a_ij*a_jk / a_jl;
-	  }	    
-	}
-      }
+                //compute sum a_jl
+                for(I ll = Sp[j]; ll < Sp[j+1]; ll++){
+                    if(C_i[Sj[ll]]){
+                        Sj_intersects_Ci = true;
+                        a_jl += Sx[ll];
+                    }	    
+                }
 
-      //Step 6
-      if(Sj_intersects_Ci){
-	for(int jj = Sp[i]; jj < Sp[i+1]; jj++){ 
-	  int j = Sj[jj];
-	  if(NodeSets[j] == C_NODE){
-	    Bj->push_back(j);
-	    Bx->push_back(-d_k[j]/d_i);      
-	  }
-	}	
-	Bp->push_back(Bj->size());
-      } else { //make i a C_NODE
-	NodeSets[i] = C_NODE;
-	Bj->push_back(i);
-	Bx->push_back(1);      
-	Bp->push_back(Bj->size());
-      }
-      
+                if(!Sj_intersects_Ci){ break; }
 
-      //Clear C_i,d_k
-      for(int jj = Sp[i]; jj < Sp[i+1]; jj++){ 
-	int j = Sj[jj];
-	C_i[j] = false;
-	d_k[j] = 0;
-      }
+                for(I kk = Sp[j]; kk < Sp[j+1]; kk++){
+                    I   k = Sj[kk];
+                    T  a_jk = Sx[kk];
+                    if(C_i[k]){
+                        d_k[k] += a_ij*a_jk / a_jl;
+                    }	    
+                }
+            }
 
-    }
-      
-  }
+            //Step 6
+            if(Sj_intersects_Ci){
+                for(I jj = Sp[i]; jj < Sp[i+1]; jj++){ 
+                    I j = Sj[jj];
+                    if(NodeSets[j] == C_NODE){
+                        Bj->push_back(j);
+                        Bx->push_back(-d_k[j]/d_i);      
+                    }
+                }	
+                Bp->push_back(Bj->size());
+            } else { //make i a C_NODE
+                NodeSets[i] = C_NODE;
+                Bj->push_back(i);
+                Bx->push_back(1);      
+                Bp->push_back(Bj->size());
+            }
 
-  //for each c-node, determine its index in the coarser lvl
-  std::vector<int> cnode_index(n_nodes,-1);
-  int n_cnodes = 0;
-  for(int i = 0; i < n_nodes; i++){
-    if(NodeSets[i] == C_NODE)
-      cnode_index[i] = n_cnodes++;
-  }
-  //map old C indices to coarse indices
-  for(std::vector<int>::iterator i = Bj->begin(); i != Bj->end(); i++){
-    *i = cnode_index[*i];
-  }
 
-  
-  
+            //Clear C_i,d_k
+            for(I jj = Sp[i]; jj < Sp[i+1]; jj++){ 
+                I j = Sj[jj];
+                C_i[j] = false;
+                d_k[j] = 0;
+            }
+
+        }
+
+    }
+
+    //for each c-node, determine its index in the coarser lvl
+    std::vector<I> cnode_index(n_nodes,-1);
+    I n_cnodes = 0;
+    for(I i = 0; i < n_nodes; i++){
+        if(NodeSets[i] == C_NODE){
+            cnode_index[i] = n_cnodes++;
+        }
+    }
+    //map old C indices to coarse indices
+    for(typename std::vector<I>::iterator iter = Bj->begin(); iter != Bj->end(); iter++){
+        *iter = cnode_index[*iter];
+    }
 }
 
 #endif

Modified: trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h	2007-12-03 22:01:43 UTC (rev 3612)
+++ trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h	2007-12-03 23:30:01 UTC (rev 3613)
@@ -10,196 +10,196 @@
 //#define DEBUG
 
 
-template<class T>
-void sa_strong_connections(const int n_row,
-			   const T epsilon,
-			   const int Ap[], const int Aj[], const T Ax[],
-			   std::vector<int> * Sp, std::vector<int> * Sj, std::vector<T> * Sx){
-  //Sp,Sj form a CSR representation where the i-th row contains
-  //the indices of all the strong connections from node i
-  Sp->push_back(0);
+template<class I, class T>
+void sa_strong_connections(const I n_row,
+        const T epsilon,
+        const I Ap[], const I Aj[], const T Ax[],
+        std::vector<I> * Sp, std::vector<I> * Sj, std::vector<T> * Sx){
+    //Sp,Sj form a CSR representation where the i-th row contains
+    //the indices of all the strong connections from node i
+    Sp->push_back(0);
 
-  //compute diagonal values
-  std::vector<T> diags(n_row,T(0));
-  for(int i = 0; i < n_row; i++){
-    int row_start = Ap[i];
-    int row_end   = Ap[i+1];
-    for(int jj = row_start; jj < row_end; jj++){
-        if(Aj[jj] == i){
-        	diags[i] = Ax[jj];
-        	break;
-        }
-    }    
-  }
+    //compute diagonal values
+    std::vector<T> diags(n_row,T(0));
+    for(I i = 0; i < n_row; i++){
+        I row_start = Ap[i];
+        I row_end   = Ap[i+1];
+        for(I jj = row_start; jj < row_end; jj++){
+            if(Aj[jj] == i){
+                diags[i] = Ax[jj];
+                break;
+            }
+        }    
+    }
 
 #ifdef DEBUG
-  for(int i = 0; i < n_row; i++){ assert(diags[i] > 0); }
+    for(I i = 0; i < n_row; i++){ assert(diags[i] > 0); }
 #endif
 
-  for(int i = 0; i < n_row; i++){
-    int row_start = Ap[i];
-    int row_end   = Ap[i+1];
+    for(I i = 0; i < n_row; i++){
+        I row_start = Ap[i];
+        I row_end   = Ap[i+1];
 
-    T eps_Aii = epsilon*epsilon*diags[i];
+        T eps_Aii = epsilon*epsilon*diags[i];
 
-    T weak_sum = 0.0;
+        T weak_sum = 0.0;
 
-    for(int jj = row_start; jj < row_end; jj++){
-      const int   j = Aj[jj];
-      const T   Aij = Ax[jj];
+        for(I jj = row_start; jj < row_end; jj++){
+            const I   j = Aj[jj];
+            const T Aij = Ax[jj];
 
-      if(i == j){continue;} //skip diagonal until end of row
+            if(i == j){continue;} //skip diagonal until end of row
 
-      //  |A(i,j)| < epsilon * sqrt(|A(i,i)|*|A(j,j)|) 
-      if(Aij*Aij >= std::abs(eps_Aii * diags[j])){    
-          Sj->push_back(j);
-          Sx->push_back(Aij);
-      } else {
-          weak_sum += Aij;
-      }
+            //  |A(i,j)| < epsilon * sqrt(|A(i,i)|*|A(j,j)|) 
+            if(Aij*Aij >= std::abs(eps_Aii * diags[j])){    
+                Sj->push_back(j);
+                Sx->push_back(Aij);
+            } else {
+                weak_sum += Aij;
+            }
+        }
+        //Add modified diagonal entry
+        Sj->push_back(i);
+        Sx->push_back(diags[i] + weak_sum); //filtered matrix
+
+        Sp->push_back(Sj->size());
     }
-    //Add modified diagonal entry
-    Sj->push_back(i);
-    Sx->push_back(diags[i] + weak_sum); //filtered matrix
-
-    Sp->push_back(Sj->size());
-  }
 }
 
-
-void sa_get_aggregates(const int n_row,
-        		       const int Ap[], const int Aj[],
-		               std::vector<int> * Bj)
+    template <class I>
+void sa_get_aggregates(const I n_row,
+        const I Ap[], const I Aj[],
+        std::vector<I> * Bj)
 {
-  std::vector<int> aggregates(n_row,-1);
+    std::vector<I> aggregates(n_row,-1);
 
-  int num_aggregates = 0;
+    I num_aggregates = 0;
 
-  //Pass #1
-  for(int i = 0; i < n_row; i++){
-    if(aggregates[i] >= 0){ continue; } //already marked
+    //Pass #1
+    for(I i = 0; i < n_row; i++){
+        if(aggregates[i] >= 0){ continue; } //already marked
 
-    const int row_start = Ap[i];
-    const int row_end   = Ap[i+1];
-    
-    //Determine whether all neighbors of this node are free (not already aggregates)
-    bool free_neighborhood = true;
-    for(int jj = row_start; jj < row_end; jj++){
-      if(aggregates[Aj[jj]] >= 0){
-    	free_neighborhood = false;
-	    break;
-      }
-    }    
-    if(!free_neighborhood){ continue; } //bail out
+        const I row_start = Ap[i];
+        const I row_end   = Ap[i+1];
 
-    //Make an aggregate out of this node and its strong neigbors
-    aggregates[i] = num_aggregates;
-    for(int jj = row_start; jj < row_end; jj++){
-      aggregates[Aj[jj]] = num_aggregates;
+        //Determine whether all neighbors of this node are free (not already aggregates)
+        bool free_neighborhood = true;
+        for(I jj = row_start; jj < row_end; jj++){
+            if(aggregates[Aj[jj]] >= 0){
+                free_neighborhood = false;
+                break;
+            }
+        }    
+        if(!free_neighborhood){ continue; } //bail out
+
+        //Make an aggregate out of this node and its strong neigbors
+        aggregates[i] = num_aggregates;
+        for(I jj = row_start; jj < row_end; jj++){
+            aggregates[Aj[jj]] = num_aggregates;
+        }
+        num_aggregates++;
     }
-    num_aggregates++;
-  }
 
 
-  //Pass #2
-  std::vector<int> aggregates_copy(aggregates);
-  for(int i = 0; i < n_row; i++){
-    if(aggregates[i] >= 0){ continue; } //already marked
+    //Pass #2
+    std::vector<I> aggregates_copy(aggregates);
+    for(I i = 0; i < n_row; i++){
+        if(aggregates[i] >= 0){ continue; } //already marked
 
-    const int row_start = Ap[i];
-    const int row_end   = Ap[i+1];
-    
-    for(int jj = row_start; jj < row_end; jj++){
-      const int j = Aj[jj];
+        const I row_start = Ap[i];
+        const I row_end   = Ap[i+1];
 
-      if(aggregates_copy[j] >= 0){
-    	aggregates[i] = aggregates_copy[j];
-	    break;
-      }
-    }    
-  }
+        for(I jj = row_start; jj < row_end; jj++){
+            const I j = Aj[jj];
 
+            if(aggregates_copy[j] >= 0){
+                aggregates[i] = aggregates_copy[j];
+                break;
+            }
+        }    
+    }
 
-  //Pass #3
-  for(int i = 0; i < n_row; i++){
-    if(aggregates[i] >= 0){ continue; } //already marked
-  
-    const int row_start = Ap[i];
-    const int row_end   = Ap[i+1];
-    
-    aggregates[i] = num_aggregates;
 
-    for(int jj = row_start; jj < row_end; jj++){
-      const int j = Aj[jj];
+    //Pass #3
+    for(I i = 0; i < n_row; i++){
+        if(aggregates[i] >= 0){ continue; } //already marked
 
-      if(aggregates[j] < 0){ //unmarked neighbors
-    	aggregates[j] = num_aggregates;
-      }
-    }  
-    num_aggregates++;
-  }
+        const I row_start = Ap[i];
+        const I row_end   = Ap[i+1];
 
-#ifdef DEBUG
-  for(int i = 0; i < n_row; i++){ assert(aggregates[i] >= 0 && aggregates[i] < num_aggregates); }
-#endif
-  
-  aggregates.swap(*Bj);  
-}
+        aggregates[i] = num_aggregates;
 
+        for(I jj = row_start; jj < row_end; jj++){
+            const I j = Aj[jj];
 
+            if(aggregates[j] < 0){ //unmarked neighbors
+                aggregates[j] = num_aggregates;
+            }
+        }  
+        num_aggregates++;
+    }
 
+#ifdef DEBUG
+    for(I i = 0; i < n_row; i++){ assert(aggregates[i] >= 0 && aggregates[i] < num_aggregates); }
+#endif
 
+    aggregates.swap(*Bj);  
+}
 
 
-template<class T>
-void sa_smoother(const int n_row,
-		 const T   omega,
-		 const int Ap[], const int Aj[], const T Ax[],
-		 const int Sp[], const int Sj[], const T Sx[],
-		 std::vector<int> * Bp, std::vector<int> * Bj, std::vector<T> * Bx){
 
 
-  //compute filtered diagonal
-  std::vector<T> diags(n_row,0);
-  
-  for(int i = 0; i < n_row; i++){
-    int row_start = Ap[i];
-    int row_end   = Ap[i+1];
-    for(int jj = row_start; jj < row_end; jj++){
-      diags[i] += Ax[jj];
-    }    
-  }
-  for(int i = 0; i < n_row; i++){
-    int row_start = Sp[i];
-    int row_end   = Sp[i+1];
-    for(int jj = row_start; jj < row_end; jj++){
-      diags[i] -= Sx[jj];
-    }    
-  }
-  
-#ifdef DEBUG
-  for(int i = 0; i < n_row; i++){ assert(diags[i] > 0); }
-#endif
 
 
-  //compute omega Jacobi smoother
-  Bp->push_back(0);
-  for(int i = 0; i < n_row; i++){
-    int row_start = Sp[i];
-    int row_end   = Sp[i+1];
-    const T row_scale = -omega/diags[i];
+//template<class T>
+//void sa_smoother(const int n_row,
+//		 const T   omega,
+//		 const int Ap[], const int Aj[], const T Ax[],
+//		 const int Sp[], const int Sj[], const T Sx[],
+//		 std::vector<int> * Bp, std::vector<int> * Bj, std::vector<T> * Bx){
+//
+//
+//  //compute filtered diagonal
+//  std::vector<T> diags(n_row,0);
+//  
+//  for(int i = 0; i < n_row; i++){
+//    int row_start = Ap[i];
+//    int row_end   = Ap[i+1];
+//    for(int jj = row_start; jj < row_end; jj++){
+//      diags[i] += Ax[jj];
+//    }    
+//  }
+//  for(int i = 0; i < n_row; i++){
+//    int row_start = Sp[i];
+//    int row_end   = Sp[i+1];
+//    for(int jj = row_start; jj < row_end; jj++){
+//      diags[i] -= Sx[jj];
+//    }    
+//  }
+//  
+//#ifdef DEBUG
+//  for(int i = 0; i < n_row; i++){ assert(diags[i] > 0); }
+//#endif
+//
+//
+//  //compute omega Jacobi smoother
+//  Bp->push_back(0);
+//  for(int i = 0; i < n_row; i++){
+//    int row_start = Sp[i];
+//    int row_end   = Sp[i+1];
+//    const T row_scale = -omega/diags[i];
+//
+//    Bx->push_back(1.0);
+//    Bj->push_back( i );
+//    
+//    for(int jj = row_start; jj < row_end; jj++){
+//      Bx->push_back(row_scale*Sx[jj]);
+//      Bj->push_back(Sj[jj]);
+//    }    
+//    Bp->push_back(Bj->size());
+//  }
+//}
 
-    Bx->push_back(1.0);
-    Bj->push_back( i );
-    
-    for(int jj = row_start; jj < row_end; jj++){
-      Bx->push_back(row_scale*Sx[jj]);
-      Bj->push_back(Sj[jj]);
-    }    
-    Bp->push_back(Bj->size());
-  }
-}
 
 
-
 #endif




More information about the Scipy-svn mailing list