[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