[Scipy-svn] r3798 - in trunk/scipy/sandbox/multigrid: . multigridtools tests

scipy-svn@scip... scipy-svn@scip...
Tue Jan 8 00:19:24 CST 2008


Author: wnbell
Date: 2008-01-08 00:19:19 -0600 (Tue, 08 Jan 2008)
New Revision: 3798

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/relaxation.h
   trunk/scipy/sandbox/multigrid/relaxation.py
   trunk/scipy/sandbox/multigrid/tests/test_relaxation.py
Log:
added BSR gauss seidel method


Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i	2008-01-08 04:12:38 UTC (rev 3797)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i	2008-01-08 06:19:19 UTC (rev 3798)
@@ -170,6 +170,7 @@
 INSTANTIATE_BOTH(rs_strong_connections)
 INSTANTIATE_BOTH(rs_interpolation)
 INSTANTIATE_BOTH(sa_strong_connections)
+INSTANTIATE_BOTH(block_gauss_seidel)
 INSTANTIATE_BOTH(gauss_seidel)
 INSTANTIATE_BOTH(jacobi)
 

Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.py	2008-01-08 04:12:38 UTC (rev 3797)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.py	2008-01-08 06:19:19 UTC (rev 3798)
@@ -90,23 +90,32 @@
     """
   return _multigridtools.sa_strong_connections(*args)
 
+def block_gauss_seidel(*args):
+  """
+    block_gauss_seidel(int Ap, int Aj, float Ax, float x, float b, int row_start, 
+        int row_stop, int row_step, int blocksize)
+    block_gauss_seidel(int Ap, int Aj, double Ax, double x, double b, int row_start, 
+        int row_stop, int row_step, int blocksize)
+    """
+  return _multigridtools.block_gauss_seidel(*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 Ap, int Aj, float Ax, float x, float b, int row_start, 
+        int row_stop, int row_step)
+    gauss_seidel(int Ap, int Aj, double Ax, double x, double b, int row_start, 
+        int row_stop, int row_step)
     """
   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 Ap, int Aj, float Ax, float x, float b, float temp, 
+        int row_start, int row_stop, int row_step, 
+        float omega)
+    jacobi(int Ap, int Aj, double Ax, double x, double b, double temp, 
+        int row_start, int row_stop, int row_step, 
+        double omega)
     """
   return _multigridtools.jacobi(*args)
 

Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx	2008-01-08 04:12:38 UTC (rev 3797)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx	2008-01-08 06:19:19 UTC (rev 3798)
@@ -4457,28 +4457,28 @@
 }
 
 
-SWIGINTERN PyObject *_wrap_gauss_seidel__SWIG_1(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_block_gauss_seidel__SWIG_1(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  int arg1 ;
+  int *arg1 ;
   int *arg2 ;
-  int *arg3 ;
+  float *arg3 ;
   float *arg4 ;
   float *arg5 ;
-  float *arg6 ;
+  int arg6 ;
   int arg7 ;
   int arg8 ;
   int arg9 ;
-  int val1 ;
-  int ecode1 = 0 ;
+  PyArrayObject *array1 = NULL ;
+  int is_new_object1 ;
   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 *temp4 = NULL ;
+  PyArrayObject *array5 = NULL ;
+  int is_new_object5 ;
+  int val6 ;
+  int ecode6 = 0 ;
   int val7 ;
   int ecode7 = 0 ;
   int val8 ;
@@ -4495,16 +4495,21 @@
   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_int(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "gauss_seidel" "', argument " "1"" of type '" "int""'");
-  } 
-  arg1 = static_cast< int >(val1);
+  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOO:block_gauss_seidel",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8)) SWIG_fail;
   {
     npy_intp size[1] = {
       -1
     };
+    array1 = obj_to_array_contiguous_allow_conversion(obj0, PyArray_INT, &is_new_object1);
+    if (!array1 || !require_dimensions(array1,1) || !require_size(array1,size,1)
+      || !require_contiguous(array1)   || !require_native(array1)) SWIG_fail;
+    
+    arg1 = (int*) array1->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
     array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_INT, &is_new_object2);
     if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
       || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
@@ -4515,106 +4520,101 @@
     npy_intp size[1] = {
       -1
     };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_INT, &is_new_object3);
+    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_FLOAT, &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;
+    arg3 = (float*) 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;
+    temp4 = obj_to_array_no_conversion(obj3,PyArray_FLOAT);
+    if (!temp4  || !require_contiguous(temp4) || !require_native(temp4)) SWIG_fail;
+    arg4 = (float*) array_data(temp4);
   }
   {
-    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;
+    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;
     
-    arg6 = (float*) array6->data;
+    arg5 = (float*) array5->data;
   }
+  ecode6 = SWIG_AsVal_int(obj5, &val6);
+  if (!SWIG_IsOK(ecode6)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode6), "in method '" "block_gauss_seidel" "', argument " "6"" of type '" "int""'");
+  } 
+  arg6 = static_cast< int >(val6);
   ecode7 = SWIG_AsVal_int(obj6, &val7);
   if (!SWIG_IsOK(ecode7)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "gauss_seidel" "', argument " "7"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "block_gauss_seidel" "', argument " "7"" of type '" "int""'");
   } 
   arg7 = static_cast< int >(val7);
   ecode8 = SWIG_AsVal_int(obj7, &val8);
   if (!SWIG_IsOK(ecode8)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "gauss_seidel" "', argument " "8"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "block_gauss_seidel" "', argument " "8"" of type '" "int""'");
   } 
   arg8 = static_cast< int >(val8);
   ecode9 = SWIG_AsVal_int(obj8, &val9);
   if (!SWIG_IsOK(ecode9)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "gauss_seidel" "', argument " "9"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "block_gauss_seidel" "', argument " "9"" of type '" "int""'");
   } 
   arg9 = static_cast< int >(val9);
-  gauss_seidel<int,float >(arg1,(int const (*))arg2,(int const (*))arg3,(float const (*))arg4,arg5,(float const (*))arg6,arg7,arg8,arg9);
+  block_gauss_seidel<int,float >((int const (*))arg1,(int const (*))arg2,(float const (*))arg3,arg4,(float const (*))arg5,arg6,arg7,arg8,arg9);
   resultobj = SWIG_Py_Void();
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return resultobj;
 fail:
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_gauss_seidel__SWIG_2(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_block_gauss_seidel__SWIG_2(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  int arg1 ;
+  int *arg1 ;
   int *arg2 ;
-  int *arg3 ;
+  double *arg3 ;
   double *arg4 ;
   double *arg5 ;
-  double *arg6 ;
+  int arg6 ;
   int arg7 ;
   int arg8 ;
   int arg9 ;
-  int val1 ;
-  int ecode1 = 0 ;
+  PyArrayObject *array1 = NULL ;
+  int is_new_object1 ;
   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 *temp4 = NULL ;
+  PyArrayObject *array5 = NULL ;
+  int is_new_object5 ;
+  int val6 ;
+  int ecode6 = 0 ;
   int val7 ;
   int ecode7 = 0 ;
   int val8 ;
@@ -4631,16 +4631,21 @@
   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_int(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "gauss_seidel" "', argument " "1"" of type '" "int""'");
-  } 
-  arg1 = static_cast< int >(val1);
+  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOO:block_gauss_seidel",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8)) SWIG_fail;
   {
     npy_intp size[1] = {
       -1
     };
+    array1 = obj_to_array_contiguous_allow_conversion(obj0, PyArray_INT, &is_new_object1);
+    if (!array1 || !require_dimensions(array1,1) || !require_size(array1,size,1)
+      || !require_contiguous(array1)   || !require_native(array1)) SWIG_fail;
+    
+    arg1 = (int*) array1->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
     array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_INT, &is_new_object2);
     if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
       || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
@@ -4651,85 +4656,80 @@
     npy_intp size[1] = {
       -1
     };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_INT, &is_new_object3);
+    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_DOUBLE, &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;
+    arg3 = (double*) 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;
+    temp4 = obj_to_array_no_conversion(obj3,PyArray_DOUBLE);
+    if (!temp4  || !require_contiguous(temp4) || !require_native(temp4)) SWIG_fail;
+    arg4 = (double*) array_data(temp4);
   }
   {
-    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;
+    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;
     
-    arg6 = (double*) array6->data;
+    arg5 = (double*) array5->data;
   }
+  ecode6 = SWIG_AsVal_int(obj5, &val6);
+  if (!SWIG_IsOK(ecode6)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode6), "in method '" "block_gauss_seidel" "', argument " "6"" of type '" "int""'");
+  } 
+  arg6 = static_cast< int >(val6);
   ecode7 = SWIG_AsVal_int(obj6, &val7);
   if (!SWIG_IsOK(ecode7)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "gauss_seidel" "', argument " "7"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "block_gauss_seidel" "', argument " "7"" of type '" "int""'");
   } 
   arg7 = static_cast< int >(val7);
   ecode8 = SWIG_AsVal_int(obj7, &val8);
   if (!SWIG_IsOK(ecode8)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "gauss_seidel" "', argument " "8"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "block_gauss_seidel" "', argument " "8"" of type '" "int""'");
   } 
   arg8 = static_cast< int >(val8);
   ecode9 = SWIG_AsVal_int(obj8, &val9);
   if (!SWIG_IsOK(ecode9)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "gauss_seidel" "', argument " "9"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "block_gauss_seidel" "', argument " "9"" of type '" "int""'");
   } 
   arg9 = static_cast< int >(val9);
-  gauss_seidel<int,double >(arg1,(int const (*))arg2,(int const (*))arg3,(double const (*))arg4,arg5,(double const (*))arg6,arg7,arg8,arg9);
+  block_gauss_seidel<int,double >((int const (*))arg1,(int const (*))arg2,(double const (*))arg3,arg4,(double const (*))arg5,arg6,arg7,arg8,arg9);
   resultobj = SWIG_Py_Void();
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return resultobj;
 fail:
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_gauss_seidel(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_block_gauss_seidel(PyObject *self, PyObject *args) {
   int argc;
   PyObject *argv[10];
   int ii;
@@ -4742,8 +4742,7 @@
   if (argc == 9) {
     int _v;
     {
-      int res = SWIG_AsVal_int(argv[0], NULL);
-      _v = SWIG_CheckState(res);
+      _v = (is_array(argv[0]) && PyArray_CanCastSafely(PyArray_TYPE(argv[0]),PyArray_INT)) ? 1 : 0;
     }
     if (_v) {
       {
@@ -4751,7 +4750,7 @@
       }
       if (_v) {
         {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_INT)) ? 1 : 0;
+          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_FLOAT)) ? 1 : 0;
         }
         if (_v) {
           {
@@ -4763,7 +4762,8 @@
             }
             if (_v) {
               {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_FLOAT)) ? 1 : 0;
+                int res = SWIG_AsVal_int(argv[5], NULL);
+                _v = SWIG_CheckState(res);
               }
               if (_v) {
                 {
@@ -4781,7 +4781,7 @@
                       _v = SWIG_CheckState(res);
                     }
                     if (_v) {
-                      return _wrap_gauss_seidel__SWIG_1(self, args);
+                      return _wrap_block_gauss_seidel__SWIG_1(self, args);
                     }
                   }
                 }
@@ -4795,8 +4795,7 @@
   if (argc == 9) {
     int _v;
     {
-      int res = SWIG_AsVal_int(argv[0], NULL);
-      _v = SWIG_CheckState(res);
+      _v = (is_array(argv[0]) && PyArray_CanCastSafely(PyArray_TYPE(argv[0]),PyArray_INT)) ? 1 : 0;
     }
     if (_v) {
       {
@@ -4804,7 +4803,7 @@
       }
       if (_v) {
         {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_INT)) ? 1 : 0;
+          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_DOUBLE)) ? 1 : 0;
         }
         if (_v) {
           {
@@ -4816,7 +4815,8 @@
             }
             if (_v) {
               {
-                _v = (is_array(argv[5]) && PyArray_CanCastSafely(PyArray_TYPE(argv[5]),PyArray_DOUBLE)) ? 1 : 0;
+                int res = SWIG_AsVal_int(argv[5], NULL);
+                _v = SWIG_CheckState(res);
               }
               if (_v) {
                 {
@@ -4834,7 +4834,7 @@
                       _v = SWIG_CheckState(res);
                     }
                     if (_v) {
-                      return _wrap_gauss_seidel__SWIG_2(self, args);
+                      return _wrap_block_gauss_seidel__SWIG_2(self, args);
                     }
                   }
                 }
@@ -4847,44 +4847,406 @@
   }
   
 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");
+  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number of arguments for overloaded function 'block_gauss_seidel'.\n  Possible C/C++ prototypes are:\n""    block_gauss_seidel<(int,float)>(int const [],int const [],float const [],float [],float const [],int const,int const,int const,int const)\n""    block_gauss_seidel<(int,double)>(int const [],int const [],double const [],double [],double const [],int const,int const,int const,int const)\n");
   return NULL;
 }
 
 
+SWIGINTERN PyObject *_wrap_gauss_seidel__SWIG_1(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  int *arg1 ;
+  int *arg2 ;
+  float *arg3 ;
+  float *arg4 ;
+  float *arg5 ;
+  int arg6 ;
+  int arg7 ;
+  int arg8 ;
+  PyArrayObject *array1 = NULL ;
+  int is_new_object1 ;
+  PyArrayObject *array2 = NULL ;
+  int is_new_object2 ;
+  PyArrayObject *array3 = NULL ;
+  int is_new_object3 ;
+  PyArrayObject *temp4 = NULL ;
+  PyArrayObject *array5 = NULL ;
+  int is_new_object5 ;
+  int val6 ;
+  int ecode6 = 0 ;
+  int val7 ;
+  int ecode7 = 0 ;
+  int val8 ;
+  int ecode8 = 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 ;
+  
+  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOO:gauss_seidel",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7)) SWIG_fail;
+  {
+    npy_intp size[1] = {
+      -1
+    };
+    array1 = obj_to_array_contiguous_allow_conversion(obj0, PyArray_INT, &is_new_object1);
+    if (!array1 || !require_dimensions(array1,1) || !require_size(array1,size,1)
+      || !require_contiguous(array1)   || !require_native(array1)) SWIG_fail;
+    
+    arg1 = (int*) array1->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
+    array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_INT, &is_new_object2);
+    if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
+      || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
+    
+    arg2 = (int*) array2->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
+    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_FLOAT, &is_new_object3);
+    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
+      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
+    
+    arg3 = (float*) array3->data;
+  }
+  {
+    temp4 = obj_to_array_no_conversion(obj3,PyArray_FLOAT);
+    if (!temp4  || !require_contiguous(temp4) || !require_native(temp4)) SWIG_fail;
+    arg4 = (float*) array_data(temp4);
+  }
+  {
+    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;
+  }
+  ecode6 = SWIG_AsVal_int(obj5, &val6);
+  if (!SWIG_IsOK(ecode6)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode6), "in method '" "gauss_seidel" "', argument " "6"" of type '" "int""'");
+  } 
+  arg6 = static_cast< int >(val6);
+  ecode7 = SWIG_AsVal_int(obj6, &val7);
+  if (!SWIG_IsOK(ecode7)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "gauss_seidel" "', argument " "7"" of type '" "int""'");
+  } 
+  arg7 = static_cast< int >(val7);
+  ecode8 = SWIG_AsVal_int(obj7, &val8);
+  if (!SWIG_IsOK(ecode8)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "gauss_seidel" "', argument " "8"" of type '" "int""'");
+  } 
+  arg8 = static_cast< int >(val8);
+  gauss_seidel<int,float >((int const (*))arg1,(int const (*))arg2,(float const (*))arg3,arg4,(float const (*))arg5,arg6,arg7,arg8);
+  resultobj = SWIG_Py_Void();
+  {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
+    if (is_new_object2 && array2) Py_DECREF(array2);
+  }
+  {
+    if (is_new_object3 && array3) Py_DECREF(array3);
+  }
+  {
+    if (is_new_object5 && array5) Py_DECREF(array5);
+  }
+  return resultobj;
+fail:
+  {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
+    if (is_new_object2 && array2) Py_DECREF(array2);
+  }
+  {
+    if (is_new_object3 && array3) Py_DECREF(array3);
+  }
+  {
+    if (is_new_object5 && array5) Py_DECREF(array5);
+  }
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_gauss_seidel__SWIG_2(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  int *arg1 ;
+  int *arg2 ;
+  double *arg3 ;
+  double *arg4 ;
+  double *arg5 ;
+  int arg6 ;
+  int arg7 ;
+  int arg8 ;
+  PyArrayObject *array1 = NULL ;
+  int is_new_object1 ;
+  PyArrayObject *array2 = NULL ;
+  int is_new_object2 ;
+  PyArrayObject *array3 = NULL ;
+  int is_new_object3 ;
+  PyArrayObject *temp4 = NULL ;
+  PyArrayObject *array5 = NULL ;
+  int is_new_object5 ;
+  int val6 ;
+  int ecode6 = 0 ;
+  int val7 ;
+  int ecode7 = 0 ;
+  int val8 ;
+  int ecode8 = 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 ;
+  
+  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOO:gauss_seidel",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7)) SWIG_fail;
+  {
+    npy_intp size[1] = {
+      -1
+    };
+    array1 = obj_to_array_contiguous_allow_conversion(obj0, PyArray_INT, &is_new_object1);
+    if (!array1 || !require_dimensions(array1,1) || !require_size(array1,size,1)
+      || !require_contiguous(array1)   || !require_native(array1)) SWIG_fail;
+    
+    arg1 = (int*) array1->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
+    array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_INT, &is_new_object2);
+    if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
+      || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
+    
+    arg2 = (int*) array2->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
+    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_DOUBLE, &is_new_object3);
+    if (!array3 || !require_dimensions(array3,1) || !require_size(array3,size,1)
+      || !require_contiguous(array3)   || !require_native(array3)) SWIG_fail;
+    
+    arg3 = (double*) array3->data;
+  }
+  {
+    temp4 = obj_to_array_no_conversion(obj3,PyArray_DOUBLE);
+    if (!temp4  || !require_contiguous(temp4) || !require_native(temp4)) SWIG_fail;
+    arg4 = (double*) array_data(temp4);
+  }
+  {
+    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;
+  }
+  ecode6 = SWIG_AsVal_int(obj5, &val6);
+  if (!SWIG_IsOK(ecode6)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode6), "in method '" "gauss_seidel" "', argument " "6"" of type '" "int""'");
+  } 
+  arg6 = static_cast< int >(val6);
+  ecode7 = SWIG_AsVal_int(obj6, &val7);
+  if (!SWIG_IsOK(ecode7)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "gauss_seidel" "', argument " "7"" of type '" "int""'");
+  } 
+  arg7 = static_cast< int >(val7);
+  ecode8 = SWIG_AsVal_int(obj7, &val8);
+  if (!SWIG_IsOK(ecode8)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "gauss_seidel" "', argument " "8"" of type '" "int""'");
+  } 
+  arg8 = static_cast< int >(val8);
+  gauss_seidel<int,double >((int const (*))arg1,(int const (*))arg2,(double const (*))arg3,arg4,(double const (*))arg5,arg6,arg7,arg8);
+  resultobj = SWIG_Py_Void();
+  {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
+    if (is_new_object2 && array2) Py_DECREF(array2);
+  }
+  {
+    if (is_new_object3 && array3) Py_DECREF(array3);
+  }
+  {
+    if (is_new_object5 && array5) Py_DECREF(array5);
+  }
+  return resultobj;
+fail:
+  {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
+    if (is_new_object2 && array2) Py_DECREF(array2);
+  }
+  {
+    if (is_new_object3 && array3) Py_DECREF(array3);
+  }
+  {
+    if (is_new_object5 && array5) Py_DECREF(array5);
+  }
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_gauss_seidel(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;
+    {
+      _v = (is_array(argv[0]) && PyArray_CanCastSafely(PyArray_TYPE(argv[0]),PyArray_INT)) ? 1 : 0;
+    }
+    if (_v) {
+      {
+        _v = (is_array(argv[1]) && PyArray_CanCastSafely(PyArray_TYPE(argv[1]),PyArray_INT)) ? 1 : 0;
+      }
+      if (_v) {
+        {
+          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_FLOAT)) ? 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) {
+              {
+                int res = SWIG_AsVal_int(argv[5], NULL);
+                _v = SWIG_CheckState(res);
+              }
+              if (_v) {
+                {
+                  int res = SWIG_AsVal_int(argv[6], NULL);
+                  _v = SWIG_CheckState(res);
+                }
+                if (_v) {
+                  {
+                    int res = SWIG_AsVal_int(argv[7], NULL);
+                    _v = SWIG_CheckState(res);
+                  }
+                  if (_v) {
+                    return _wrap_gauss_seidel__SWIG_1(self, args);
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+  if (argc == 8) {
+    int _v;
+    {
+      _v = (is_array(argv[0]) && PyArray_CanCastSafely(PyArray_TYPE(argv[0]),PyArray_INT)) ? 1 : 0;
+    }
+    if (_v) {
+      {
+        _v = (is_array(argv[1]) && PyArray_CanCastSafely(PyArray_TYPE(argv[1]),PyArray_INT)) ? 1 : 0;
+      }
+      if (_v) {
+        {
+          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_DOUBLE)) ? 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) {
+              {
+                int res = SWIG_AsVal_int(argv[5], NULL);
+                _v = SWIG_CheckState(res);
+              }
+              if (_v) {
+                {
+                  int res = SWIG_AsVal_int(argv[6], NULL);
+                  _v = SWIG_CheckState(res);
+                }
+                if (_v) {
+                  {
+                    int res = SWIG_AsVal_int(argv[7], NULL);
+                    _v = SWIG_CheckState(res);
+                  }
+                  if (_v) {
+                    return _wrap_gauss_seidel__SWIG_2(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 [],float const [],float [],float const [],int const,int const,int const)\n""    gauss_seidel<(int,double)>(int const [],int const [],double const [],double [],double const [],int const,int const,int const)\n");
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_jacobi__SWIG_1(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  int arg1 ;
+  int *arg1 ;
   int *arg2 ;
-  int *arg3 ;
+  float *arg3 ;
   float *arg4 ;
   float *arg5 ;
   float *arg6 ;
-  float *arg7 ;
+  int arg7 ;
   int arg8 ;
   int arg9 ;
-  int arg10 ;
-  float arg11 ;
-  int val1 ;
-  int ecode1 = 0 ;
+  float arg10 ;
+  PyArrayObject *array1 = NULL ;
+  int is_new_object1 ;
   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 ;
+  PyArrayObject *temp4 = NULL ;
+  PyArrayObject *array5 = NULL ;
+  int is_new_object5 ;
+  PyArrayObject *temp6 = NULL ;
+  int val7 ;
+  int ecode7 = 0 ;
   int val8 ;
   int ecode8 = 0 ;
   int val9 ;
   int ecode9 = 0 ;
-  int val10 ;
+  float val10 ;
   int ecode10 = 0 ;
-  float val11 ;
-  int ecode11 = 0 ;
   PyObject * obj0 = 0 ;
   PyObject * obj1 = 0 ;
   PyObject * obj2 = 0 ;
@@ -4895,18 +5257,22 @@
   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_int(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "jacobi" "', argument " "1"" of type '" "int""'");
-  } 
-  arg1 = static_cast< int >(val1);
+  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOOO:jacobi",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8,&obj9)) SWIG_fail;
   {
     npy_intp size[1] = {
       -1
     };
+    array1 = obj_to_array_contiguous_allow_conversion(obj0, PyArray_INT, &is_new_object1);
+    if (!array1 || !require_dimensions(array1,1) || !require_size(array1,size,1)
+      || !require_contiguous(array1)   || !require_native(array1)) SWIG_fail;
+    
+    arg1 = (int*) array1->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
     array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_INT, &is_new_object2);
     if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
       || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
@@ -4917,42 +5283,37 @@
     npy_intp size[1] = {
       -1
     };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_INT, &is_new_object3);
+    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_FLOAT, &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;
+    arg3 = (float*) 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;
+    temp4 = obj_to_array_no_conversion(obj3,PyArray_FLOAT);
+    if (!temp4  || !require_contiguous(temp4) || !require_native(temp4)) SWIG_fail;
+    arg4 = (float*) array_data(temp4);
   }
   {
-    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;
+    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;
     
-    arg6 = (float*) array6->data;
+    arg5 = (float*) array5->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);
+    temp6 = obj_to_array_no_conversion(obj5,PyArray_FLOAT);
+    if (!temp6  || !require_contiguous(temp6) || !require_native(temp6)) SWIG_fail;
+    arg6 = (float*) array_data(temp6);
   }
+  ecode7 = SWIG_AsVal_int(obj6, &val7);
+  if (!SWIG_IsOK(ecode7)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "jacobi" "', argument " "7"" of type '" "int""'");
+  } 
+  arg7 = static_cast< int >(val7);
   ecode8 = SWIG_AsVal_int(obj7, &val8);
   if (!SWIG_IsOK(ecode8)) {
     SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "jacobi" "', argument " "8"" of type '" "int""'");
@@ -4963,81 +5324,73 @@
     SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "jacobi" "', argument " "9"" of type '" "int""'");
   } 
   arg9 = static_cast< int >(val9);
-  ecode10 = SWIG_AsVal_int(obj9, &val10);
+  ecode10 = SWIG_AsVal_float(obj9, &val10);
   if (!SWIG_IsOK(ecode10)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode10), "in method '" "jacobi" "', argument " "10"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode10), "in method '" "jacobi" "', argument " "10"" of type '" "float""'");
   } 
-  arg10 = static_cast< int >(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<int,float >(arg1,(int const (*))arg2,(int const (*))arg3,(float const (*))arg4,arg5,(float const (*))arg6,arg7,arg8,arg9,arg10,arg11);
+  arg10 = static_cast< float >(val10);
+  jacobi<int,float >((int const (*))arg1,(int const (*))arg2,(float const (*))arg3,arg4,(float const (*))arg5,arg6,arg7,arg8,arg9,arg10);
   resultobj = SWIG_Py_Void();
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return resultobj;
 fail:
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return NULL;
 }
 
 
 SWIGINTERN PyObject *_wrap_jacobi__SWIG_2(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  int arg1 ;
+  int *arg1 ;
   int *arg2 ;
-  int *arg3 ;
+  double *arg3 ;
   double *arg4 ;
   double *arg5 ;
   double *arg6 ;
-  double *arg7 ;
+  int arg7 ;
   int arg8 ;
   int arg9 ;
-  int arg10 ;
-  double arg11 ;
-  int val1 ;
-  int ecode1 = 0 ;
+  double arg10 ;
+  PyArrayObject *array1 = NULL ;
+  int is_new_object1 ;
   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 ;
+  PyArrayObject *temp4 = NULL ;
+  PyArrayObject *array5 = NULL ;
+  int is_new_object5 ;
+  PyArrayObject *temp6 = NULL ;
+  int val7 ;
+  int ecode7 = 0 ;
   int val8 ;
   int ecode8 = 0 ;
   int val9 ;
   int ecode9 = 0 ;
-  int val10 ;
+  double val10 ;
   int ecode10 = 0 ;
-  double val11 ;
-  int ecode11 = 0 ;
   PyObject * obj0 = 0 ;
   PyObject * obj1 = 0 ;
   PyObject * obj2 = 0 ;
@@ -5048,18 +5401,22 @@
   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_int(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "jacobi" "', argument " "1"" of type '" "int""'");
-  } 
-  arg1 = static_cast< int >(val1);
+  if (!PyArg_ParseTuple(args,(char *)"OOOOOOOOOO:jacobi",&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6,&obj7,&obj8,&obj9)) SWIG_fail;
   {
     npy_intp size[1] = {
       -1
     };
+    array1 = obj_to_array_contiguous_allow_conversion(obj0, PyArray_INT, &is_new_object1);
+    if (!array1 || !require_dimensions(array1,1) || !require_size(array1,size,1)
+      || !require_contiguous(array1)   || !require_native(array1)) SWIG_fail;
+    
+    arg1 = (int*) array1->data;
+  }
+  {
+    npy_intp size[1] = {
+      -1
+    };
     array2 = obj_to_array_contiguous_allow_conversion(obj1, PyArray_INT, &is_new_object2);
     if (!array2 || !require_dimensions(array2,1) || !require_size(array2,size,1)
       || !require_contiguous(array2)   || !require_native(array2)) SWIG_fail;
@@ -5070,42 +5427,37 @@
     npy_intp size[1] = {
       -1
     };
-    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_INT, &is_new_object3);
+    array3 = obj_to_array_contiguous_allow_conversion(obj2, PyArray_DOUBLE, &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;
+    arg3 = (double*) 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;
+    temp4 = obj_to_array_no_conversion(obj3,PyArray_DOUBLE);
+    if (!temp4  || !require_contiguous(temp4) || !require_native(temp4)) SWIG_fail;
+    arg4 = (double*) array_data(temp4);
   }
   {
-    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;
+    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;
     
-    arg6 = (double*) array6->data;
+    arg5 = (double*) array5->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);
+    temp6 = obj_to_array_no_conversion(obj5,PyArray_DOUBLE);
+    if (!temp6  || !require_contiguous(temp6) || !require_native(temp6)) SWIG_fail;
+    arg6 = (double*) array_data(temp6);
   }
+  ecode7 = SWIG_AsVal_int(obj6, &val7);
+  if (!SWIG_IsOK(ecode7)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "jacobi" "', argument " "7"" of type '" "int""'");
+  } 
+  arg7 = static_cast< int >(val7);
   ecode8 = SWIG_AsVal_int(obj7, &val8);
   if (!SWIG_IsOK(ecode8)) {
     SWIG_exception_fail(SWIG_ArgError(ecode8), "in method '" "jacobi" "', argument " "8"" of type '" "int""'");
@@ -5116,63 +5468,57 @@
     SWIG_exception_fail(SWIG_ArgError(ecode9), "in method '" "jacobi" "', argument " "9"" of type '" "int""'");
   } 
   arg9 = static_cast< int >(val9);
-  ecode10 = SWIG_AsVal_int(obj9, &val10);
+  ecode10 = SWIG_AsVal_double(obj9, &val10);
   if (!SWIG_IsOK(ecode10)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode10), "in method '" "jacobi" "', argument " "10"" of type '" "int""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode10), "in method '" "jacobi" "', argument " "10"" of type '" "double""'");
   } 
-  arg10 = static_cast< int >(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<int,double >(arg1,(int const (*))arg2,(int const (*))arg3,(double const (*))arg4,arg5,(double const (*))arg6,arg7,arg8,arg9,arg10,arg11);
+  arg10 = static_cast< double >(val10);
+  jacobi<int,double >((int const (*))arg1,(int const (*))arg2,(double const (*))arg3,arg4,(double const (*))arg5,arg6,arg7,arg8,arg9,arg10);
   resultobj = SWIG_Py_Void();
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return resultobj;
 fail:
   {
+    if (is_new_object1 && array1) Py_DECREF(array1);
+  }
+  {
     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_object5 && array5) Py_DECREF(array5);
   }
-  {
-    if (is_new_object6 && array6) Py_DECREF(array6);
-  }
   return NULL;
 }
 
 
 SWIGINTERN PyObject *_wrap_jacobi(PyObject *self, PyObject *args) {
   int argc;
-  PyObject *argv[12];
+  PyObject *argv[11];
   int ii;
   
   if (!PyTuple_Check(args)) SWIG_fail;
   argc = (int)PyObject_Length(args);
-  for (ii = 0; (ii < argc) && (ii < 11); ii++) {
+  for (ii = 0; (ii < argc) && (ii < 10); ii++) {
     argv[ii] = PyTuple_GET_ITEM(args,ii);
   }
-  if (argc == 11) {
+  if (argc == 10) {
     int _v;
     {
-      int res = SWIG_AsVal_int(argv[0], NULL);
-      _v = SWIG_CheckState(res);
+      _v = (is_array(argv[0]) && PyArray_CanCastSafely(PyArray_TYPE(argv[0]),PyArray_INT)) ? 1 : 0;
     }
     if (_v) {
       {
@@ -5180,7 +5526,7 @@
       }
       if (_v) {
         {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_INT)) ? 1 : 0;
+          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_FLOAT)) ? 1 : 0;
         }
         if (_v) {
           {
@@ -5196,7 +5542,8 @@
               }
               if (_v) {
                 {
-                  _v = (is_array(argv[6]) && PyArray_CanCastSafely(PyArray_TYPE(argv[6]),PyArray_FLOAT)) ? 1 : 0;
+                  int res = SWIG_AsVal_int(argv[6], NULL);
+                  _v = SWIG_CheckState(res);
                 }
                 if (_v) {
                   {
@@ -5210,17 +5557,11 @@
                     }
                     if (_v) {
                       {
-                        int res = SWIG_AsVal_int(argv[9], NULL);
+                        int res = SWIG_AsVal_float(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_1(self, args);
-                        }
+                        return _wrap_jacobi__SWIG_1(self, args);
                       }
                     }
                   }
@@ -5232,11 +5573,10 @@
       }
     }
   }
-  if (argc == 11) {
+  if (argc == 10) {
     int _v;
     {
-      int res = SWIG_AsVal_int(argv[0], NULL);
-      _v = SWIG_CheckState(res);
+      _v = (is_array(argv[0]) && PyArray_CanCastSafely(PyArray_TYPE(argv[0]),PyArray_INT)) ? 1 : 0;
     }
     if (_v) {
       {
@@ -5244,7 +5584,7 @@
       }
       if (_v) {
         {
-          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_INT)) ? 1 : 0;
+          _v = (is_array(argv[2]) && PyArray_CanCastSafely(PyArray_TYPE(argv[2]),PyArray_DOUBLE)) ? 1 : 0;
         }
         if (_v) {
           {
@@ -5260,7 +5600,8 @@
               }
               if (_v) {
                 {
-                  _v = (is_array(argv[6]) && PyArray_CanCastSafely(PyArray_TYPE(argv[6]),PyArray_DOUBLE)) ? 1 : 0;
+                  int res = SWIG_AsVal_int(argv[6], NULL);
+                  _v = SWIG_CheckState(res);
                 }
                 if (_v) {
                   {
@@ -5274,17 +5615,11 @@
                     }
                     if (_v) {
                       {
-                        int res = SWIG_AsVal_int(argv[9], NULL);
+                        int res = SWIG_AsVal_double(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_2(self, args);
-                        }
+                        return _wrap_jacobi__SWIG_2(self, args);
                       }
                     }
                   }
@@ -5298,7 +5633,7 @@
   }
   
 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");
+  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 [],float const [],float [],float const [],float [],int const,int const,int const,float const)\n""    jacobi<(int,double)>(int const [],int const [],double const [],double [],double const [],double [],int const,int const,int const,double const)\n");
   return NULL;
 }
 
@@ -5329,19 +5664,25 @@
 		"    std::vector<(int)> Sp, std::vector<(int)> Sj, \n"
 		"    std::vector<(double)> Sx)\n"
 		""},
+	 { (char *)"block_gauss_seidel", _wrap_block_gauss_seidel, METH_VARARGS, (char *)"\n"
+		"block_gauss_seidel(int Ap, int Aj, float Ax, float x, float b, int row_start, \n"
+		"    int row_stop, int row_step, int blocksize)\n"
+		"block_gauss_seidel(int Ap, int Aj, double Ax, double x, double b, int row_start, \n"
+		"    int row_stop, int row_step, int blocksize)\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(int Ap, int Aj, float Ax, float x, float b, int row_start, \n"
+		"    int row_stop, int row_step)\n"
+		"gauss_seidel(int Ap, int Aj, double Ax, double x, double b, int row_start, \n"
+		"    int row_stop, int 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"
-		"    float temp, int row_start, int row_stop, \n"
-		"    int row_step, float omega)\n"
-		"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(int Ap, int Aj, float Ax, float x, float b, float temp, \n"
+		"    int row_start, int row_stop, int row_step, \n"
+		"    float omega)\n"
+		"jacobi(int Ap, int Aj, double Ax, double x, double b, double temp, \n"
+		"    int row_start, int row_stop, int row_step, \n"
+		"    double omega)\n"
 		""},
 	 { NULL, NULL, 0, NULL }
 };

Modified: trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h	2008-01-08 04:12:38 UTC (rev 3797)
+++ trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h	2008-01-08 06:19:19 UTC (rev 3798)
@@ -5,8 +5,7 @@
 #include <iostream>
 
 template<class I, class T>
-void gauss_seidel(const I n_row,
-			      const I Ap[], 
+void gauss_seidel(const I Ap[], 
                   const I Aj[], 
                   const T Ax[],
                         T  x[],
@@ -36,9 +35,59 @@
     }
 }
 
+
 template<class I, class T>
-void jacobi(const I n_row,
-			const I Ap[], 
+void block_gauss_seidel(const I Ap[], 
+                        const I Aj[], 
+                        const T Ax[],
+                              T  x[],
+                        const T  b[],
+                        const I row_start,
+                        const I row_stop,
+                        const I row_step,
+                        const I blocksize)
+{
+    const I B2 = blocksize * blocksize;
+    
+    for(I i = row_start; i != row_stop; i += row_step) {
+        I start = Ap[i];
+        I end   = Ap[i+1];
+
+        for(I bi = 0; bi < blocksize; bi++){
+            T rsum = 0;
+            T diag = 0;
+
+            for(I jj = start; jj < end; jj++){
+                I j = Aj[jj];
+                const T * block_row = Ax + B2*jj + blocksize*bi;
+                const T * block_x   = x + blocksize * j;
+
+                if (i == j){
+                    //diagonal block
+                    diag = block_row[bi];
+                    for(I bj = 0; bj < bi; bj++){
+                        rsum += block_row[bj] * block_x[bj];
+                    }
+                    for(I bj = bi+1; bj < blocksize; bj++){
+                        rsum += block_row[bj] * block_x[bj];
+                    }
+                } else {
+                    for(I bj = 0; bj < blocksize; bj++){
+                        rsum += block_row[bj] * block_x[bj];
+                    }
+                }
+            }
+
+            //TODO raise error? inform user?
+            if (diag != 0){
+                x[blocksize*i + bi] = (b[blocksize*i + bi] - rsum)/diag;
+            }
+        }
+    }
+}
+
+template<class I, class T>
+void jacobi(const I Ap[], 
             const I Aj[], 
             const T Ax[],
                   T  x[],
@@ -49,7 +98,9 @@
             const I row_step,
             const T omega)
 {
-    std::copy(x,x+n_row,temp);
+    for(I i = row_start; i != row_stop; i += row_step) {
+        temp[i] = x[i];
+    }
     
     for(I i = row_start; i != row_stop; i += row_step) {
         I start = Ap[i];

Modified: trunk/scipy/sandbox/multigrid/relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/relaxation.py	2008-01-08 04:12:38 UTC (rev 3797)
+++ trunk/scipy/sandbox/multigrid/relaxation.py	2008-01-08 06:19:19 UTC (rev 3798)
@@ -1,7 +1,30 @@
+from warnings import warn
+
+from numpy import empty_like, asarray, arange, ravel
+
 import multigridtools
-from numpy import empty_like, asarray
+from scipy.sparse import isspmatrix_csr, isspmatrix_csc, isspmatrix_bsr, \
+        csr_matrix, coo_matrix, bsr_matrix, SparseEfficiencyWarning
 
 
+#def split_ldu(A):
+#    """
+#    Return the lower triangle, diagonal, and upper triangular portions of a matrix.
+#    """
+#    if isspmatrix_csr(A) or isspmatrix_csc(A):
+#        coo = A.tocoo(copy=False)
+#    elif isspmatrix_bsr(A):
+#        M,N = A.shape
+#        R,C = A.blocksize 
+#        data = arange(len(A.indices),dtype='intc')
+#        proxy = csr_matrix((data,A.indices,A.indptr),shape=(M/R,N/C))
+#        L,D,U = split_ldu(proxy)
+#        L = bsr_matrix((A.data[L.data],L.indices,L.indptr),shape=A.shape)
+#        U = bsr_matrix((A.data[U.data],U.indices,U.indptr),shape=A.shape)
+#        D = A.data[D]
+
+
+
 def sor(A,x,b,omega,iterations=1,sweep='forward'):
     """
     Perform SOR iteration on the linear system Ax=b
@@ -30,37 +53,58 @@
          sweep      - direction of sweep:
                         'forward' (default), 'backward', or 'symmetric'
     """
-    x = asarray(x).reshape(-1)
-    b = asarray(b).reshape(-1)
+    #TODO replace pointwise BSR with block BSR
 
+    x = x.reshape(-1) #TODO warn if not inplace
+    b = ravel(b)
+
+    if isspmatrix_csr(A):
+        pass
+    elif isspmatrix_bsr(A):
+        R,C = A.blocksize
+        if R != C:
+            raise ValueError,'BSR blocks must be square'
+    else:
+        warn('implicit conversion to CSR',SparseEfficiencyWarning)
+        A = csr_matrix(A)
+
     if A.shape[0] != A.shape[1]:
-        raise ValueError,'expected symmetric matrix'
+        raise ValueError,'expected square matrix'
 
     if A.shape[1] != len(x) or len(x) != len(b):
         raise ValueError,'unexpected number of unknowns'
 
+
+    
     if sweep == 'forward':
         row_start,row_stop,row_step = 0,len(x),1
-        for iter in xrange(iterations):
-            multigridtools.gauss_seidel(A.shape[0],
-                                        A.indptr, A.indices, A.data,
-                                        x, b,
-                                        row_start, row_stop, row_step)
     elif sweep == 'backward':
-        row_start,row_stop,row_step = len(x)-1,-1,-1
-        for iter in xrange(iterations):
-            multigridtools.gauss_seidel(A.shape[0],
-                                        A.indptr, A.indices, A.data,
-                                        x, b,
-                                        row_start, row_stop, row_step)
+        row_start,row_stop,row_step = len(x)-1,-1,-1 
     elif sweep == 'symmetric':
         for iter in xrange(iterations):
             gauss_seidel(A,x,b,iterations=1,sweep='forward')
             gauss_seidel(A,x,b,iterations=1,sweep='backward')
+        return
     else:
         raise ValueError,'valid sweep directions are \'forward\', \'backward\', and \'symmetric\''
 
 
+    if isspmatrix_csr(A):
+        for iter in xrange(iterations):
+            multigridtools.gauss_seidel(A.indptr, A.indices, A.data,
+                                        x, b,
+                                        row_start, row_stop, row_step)
+    else:
+        blocksize = A.blocksize[0]
+        row_start = row_start/blocksize
+        row_stop  = row_stop/blocksize
+        for iter in xrange(iterations):
+            multigridtools.block_gauss_seidel(A.indptr, A.indices, ravel(A.data),
+                                              x, b,
+                                              row_start, row_stop, row_step,
+                                              blocksize)
+
+
 def jacobi(A,x,b,iterations=1,omega=1.0):
     """
     Perform Jacobi iteration on the linear system Ax=b
@@ -89,8 +133,7 @@
     temp = empty_like(x)
 
     for iter in xrange(iterations):
-        multigridtools.jacobi(A.shape[0],
-                              A.indptr, A.indices, A.data,
+        multigridtools.jacobi(A.indptr, A.indices, A.data,
                               x, b, temp,
                               row_start, row_stop, row_step,
                               omega)

Modified: trunk/scipy/sandbox/multigrid/tests/test_relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_relaxation.py	2008-01-08 04:12:38 UTC (rev 3797)
+++ trunk/scipy/sandbox/multigrid/tests/test_relaxation.py	2008-01-08 06:19:19 UTC (rev 3798)
@@ -2,7 +2,7 @@
 
 import numpy
 import scipy
-from scipy import arange,ones,zeros,array,allclose
+from scipy import arange,ones,zeros,array,allclose,zeros_like
 from scipy.sparse import spdiags
 
 
@@ -15,7 +15,7 @@
 class TestRelaxation(NumpyTestCase):
     def check_polynomial(self):
         N  = 3
-        A  = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A  = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x0 = arange(N).astype(numpy.float64)
         x  = x0.copy()
         b  = zeros(N)
@@ -39,87 +39,105 @@
 
     def check_jacobi(self):
         N = 1
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = zeros(N)
         jacobi(A,x,b)
         assert_almost_equal(x,array([0]))
 
         N = 3
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = zeros(N)
         b = arange(N).astype(numpy.float64)
         jacobi(A,x,b)
         assert_almost_equal(x,array([0.0,0.5,1.0]))
 
         N = 3
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = zeros(N)
         jacobi(A,x,b)
         assert_almost_equal(x,array([0.5,1.0,0.5]))
 
         N = 1
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = array([10])
         jacobi(A,x,b)
         assert_almost_equal(x,array([5]))
 
         N = 3
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = array([10,20,30])
         jacobi(A,x,b)
         assert_almost_equal(x,array([5.5,11.0,15.5]))
 
         N = 3
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         x_copy = x.copy()
         b = array([10,20,30])
         jacobi(A,x,b,omega=1.0/3.0)
         assert_almost_equal(x,2.0/3.0*x_copy + 1.0/3.0*array([5.5,11.0,15.5]))
 
+    def check_gauss_seidel_bsr(self):
+        cases = []
 
-    def check_gauss_seidel(self):
+        for N in [1,2,3,4,5,6,10]:
+            A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+            
+            divisors = [ n for n in range(1,N+1) if N % n == 0 ]
+
+            x_csr = arange(N).astype(numpy.float64)
+            b = x_csr**2
+            gauss_seidel(A,x_csr,b)
+
+            for D in divisors:
+                B = A.tobsr(blocksize=(D,D))
+                x_bsr = arange(N).astype(numpy.float64)
+                gauss_seidel(B,x_bsr,b)
+                assert_almost_equal(x_bsr,x_csr)
+               
+
+    def check_gauss_seidel_csr(self):
         N = 1
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = zeros(N)
         gauss_seidel(A,x,b)
         assert_almost_equal(x,array([0]))
 
         N = 3
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = zeros(N)
         gauss_seidel(A,x,b)
         assert_almost_equal(x,array([1.0/2.0,5.0/4.0,5.0/8.0]))
 
         N = 1
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = zeros(N)
         gauss_seidel(A,x,b,sweep='backward')
         assert_almost_equal(x,array([0]))
 
         N = 3
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = zeros(N)
         gauss_seidel(A,x,b,sweep='backward')
         assert_almost_equal(x,array([1.0/8.0,1.0/4.0,1.0/2.0]))
 
         N = 1
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = array([10])
         gauss_seidel(A,x,b)
         assert_almost_equal(x,array([5]))
 
         N = 3
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = arange(N).astype(numpy.float64)
         b = array([10,20,30])
         gauss_seidel(A,x,b)
@@ -128,7 +146,7 @@
 
         #forward and backward passes should give same result with x=ones(N),b=zeros(N)
         N = 100
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N,format='csr')
         x = ones(N)
         b = zeros(N)
         gauss_seidel(A,x,b,iterations=200,sweep='forward')



More information about the Scipy-svn mailing list