first iteration with tests. Incomplete and some failing. Just a start…

…. Some bugs already fixed.
numpy · ovillellas · Oct 19, 2012 · Oct 19, 2012 · Oct 24, 2012 · Nov 7, 2012
commit fdeafbdeddb8cb2f161ffe13cac96440e2c3359b
diff --git a/numpy/core/src/umath/gufuncs_linalg.py b/numpy/core/src/umath/gufuncs_linalg.py
@@ -48,7 +48,7 @@ def solve(A,B,**kw_args):
     """
     Solves the systems of equations AX=B or Ax=b in the inner dimensions of A/B
     """
-    if (len(B.shape) == (len(A.shape) - 1)):
+    if len(B.shape) == (len(A.shape) - 1):
         gufunc = _impl.solve1
     else:
         gufunc = _impl.solve
@@ -57,6 +57,8 @@ def solve(A,B,**kw_args):
 
 
 def svd(a, full_matrices=1, compute_uv=1 ,**kw_args):
+    m = a.shape[-2]
+    n = a.shape[-1]
     if 1 == compute_uv:
         if 1 == full_matrices:
             if m < n:
@@ -77,24 +79,24 @@ def svd(a, full_matrices=1, compute_uv=1 ,**kw_args):
 
 
 def chosolve(A, B, UPLO='L', **kw_args):
-    if len(B.shape) == len(A.shape) + 1:
+    if len(B.shape) == (len(A.shape) - 1):
         if 'L' == UPLO:
-            gufunc = chosolve1_lo
+            gufunc = _impl.chosolve1_lo
         else:
-            gufunc = chosolve1_up
+            gufunc = _impl.chosolve1_up
     else:
         if 'L' == UPLO:
-            gufunc = chosolve1_lo
+            gufunc = _impl.chosolve_lo
         else:
-            gufunc = chosolve1_up
+            gufunc = _impl.chosolve_up
 
     return gufunc(A, B, **kw_args)
 
 
 def poinv(A, UPLO='L', **kw_args):
     if 'L' == UPLO:
-        gufunc = poinv_lo
+        gufunc = _impl.poinv_lo
     else:
-        gufunc = poinv_up
+        gufunc = _impl.poinv_up
 
     return gufunc(A, **kw_args);
diff --git a/numpy/core/src/umath/umath_linalg.c.src b/numpy/core/src/umath/umath_linalg.c.src
@@ -867,7 +867,6 @@ static void
      * use blas if the stride is a multiple of datatype size in the inputs
      * it should be the common case 
      */
-
     if ((0 == (steps[3] % sot)) &&
         (0 == (steps[4] % sot))) {
         /* use blas */
@@ -2151,6 +2150,7 @@ init_@lapack_func@(GEEV_PARAMS_t* params, char jobvl, char jobvr, fortran_int n)
     params->VL = vl;
     params->VR = vr;
     params->WORK = work;
+    params->W = w;
     params->N = n;
     params->LDA = n;
     params->LDVL = n;
@@ -2247,7 +2247,6 @@ static inline void
                                 steps[1], steps[0]);
             steps += 2;
         }
-
         if ('V' == geev_params.JOBVR) {
             init_linearize_data(&vr_out,
                                 geev_params.N, geev_params.N,
@@ -2261,7 +2260,6 @@ static inline void
 
             call_@lapack_func@(&geev_params);
             process_@lapack_func@_results(&geev_params);
-
             delinearize_@COMPLEXTYPE@_matrix(*arg_iter++, geev_params.W, &w_out);
             if ('V' == geev_params.JOBVL)
                 delinearize_@COMPLEXTYPE@_matrix(*arg_iter++, geev_params.VL, &vl_out);
@@ -2400,7 +2398,7 @@ init_@lapack_func@(GESDD_PARAMS_t *params, char jobz, fortran_int m, fortran_int
     iwork = vt + vt_size;
 
     /* fix vt_column_count so that it is a valid lapack parameter (0 is not) */
-    vt_column_count = vt_column_count < 1? 1 : 0;
+    vt_column_count = vt_column_count < 1? 1 : vt_column_count;
     {
         /* compute optimal work size */
         @ftyp@ work_size_query;
@@ -2442,6 +2440,7 @@ init_@lapack_func@(GESDD_PARAMS_t *params, char jobz, fortran_int m, fortran_int
 
     return 1;
  error:
+    printf("%s failed init\n", __FUNCTION__);
     free(mem_buff2);
     free(mem_buff2);
     memset(params, 0, sizeof(*params));
@@ -2508,7 +2507,7 @@ init_@lapack_func@(GESDD_PARAMS_t *params, char jobz, fortran_int m, fortran_int
     iwork = rwork + rwork_size;                      
 
     /* fix vt_column_count so that it is a valid lapack parameter (0 is not) */
-    vt_column_count = vt_column_count < 1? 1 : 0;
+    vt_column_count = vt_column_count < 1? 1 : vt_column_count;
     {
         /* compute optimal work size */
         @ftyp@ work_size_query;
@@ -2549,6 +2548,7 @@ init_@lapack_func@(GESDD_PARAMS_t *params, char jobz, fortran_int m, fortran_int
 
     return 1;
  error:
+    printf("%s failed init\n", __FUNCTION__);
     free(mem_buff2);
     free(mem_buff);
     memset(params, 0, sizeof(*params));
@@ -2605,6 +2605,7 @@ static inline void
     if (init_@lapack_func@(&params, JOBZ, (fortran_int)dimensions[0], (fortran_int)dimensions[1])) {
         LINEARIZE_DATA_t a_in, u_out, s_out, v_out;
 
+	dump_gesdd_params("gesdd params", &params);
         init_linearize_data(&a_in, params.N, params.M, steps[1], steps[0]);
         if ('N' == params.JOBZ) {
             /* only the singular values are wanted */
@@ -3427,7 +3428,7 @@ GUFUNC_DESCRIPTOR_t gufunc_descriptors [] = {
         "    \"(i),(i)->()\" \n",
         2, 2, 1,
         inner1d_funcs,
-        equal_2_types
+        equal_3_types
     },
     {
         "innerwt",
@@ -3436,7 +3437,7 @@ GUFUNC_DESCRIPTOR_t gufunc_descriptors [] = {
         "    \"(i),(i),(i)->()\" \n",
         2, 3, 1,
         innerwt_funcs,
-        equal_3_types
+        equal_4_types
     },
     {
         "matrix_multiply",

diff --git a/numpy/core/tests/test_gufuncs_linalg.py b/numpy/core/tests/test_gufuncs_linalg.py
@@ -0,0 +1,230 @@
+"""
+Test functions for gufuncs_linalg module
+Heavily inspired (ripped in part) test_linalg
+"""
+
+# The following functions are implemented in the module "gufuncs_linalg"
+#
+# category "linalg"
+# - inv (TestInv)
+# - poinv (TestPoinv)
+# - det (TestDet)
+# - slogdet (TestDet)
+# - eig (TestEig)
+# - eigh (TestEigh)
+# - eigvals (TestEigvals)
+# - eigvalsh (TestEigvalsh)
+# - cholesky
+# - solve (TestSolve)
+# - chosolve (TestChosolve)
+# - svd (TestSVD)
+
+# ** unimplemented **
+# - qr
+# - matrix_power
+# - matrix_rank
+# - pinv
+# - lstsq
+# - tensorinv
+# - tensorsolve
+# - norm
+# - cond
+#
+# category "inspired by pdl"
+# - quadratic_form
+# - matrix_multiply3
+# - add3
+# - multiply3
+# - multiply3_add
+# - multiply_add
+# - multiply_add2
+# - multiply4
+# - multiply4_add
+#
+# category "others"
+# - convolve
+# - inner1d
+# - innerwt
+# - matrix_multiply
+
+
+import numpy as np
+
+from numpy.testing import (TestCase, assert_, assert_equal, assert_raises,
+                           assert_array_equal, assert_almost_equal,
+                           run_module_suite)
+
+from numpy import array, single, double, csingle, cdouble, dot, identity
+from numpy import multiply, inf
+import numpy.core.gufuncs_linalg as gula
+
+old_assert_almost_equal = assert_almost_equal
+
+def assert_almost_equal(a, b, **kw):
+    if a.dtype.type in (single, csingle):
+        decimal = 6
+    else:
+        decimal = 12
+    old_assert_almost_equal(a, b, decimal = decimal, **kw)
+
+
+class GeneralTestCase(object):
+    def test_single(self):
+        a = array([[[1.,2.], [3.,4.]]], dtype=single)
+        b = array([[2.,1.]], dtype=single)
+        self.do(a, b)
+
+    def test_double(self):
+        a = array([[[1.,2.], [3.,4.]]], dtype=double)
+        b = array([[2.,1.]], dtype=double)
+        self.do(a, b)
+
+    def test_double_2(self):
+        a = array([[[1.,2.], [3.,4.]]], dtype=double)
+        b = array([[[2.,1., 4.], [3.,4.,6.]]], dtype=double)
+        self.do(a, b)
+
+    def test_csingle(self):
+        a = array([[[1+2j,2+3j],[3+4j,4+5j]]], dtype=csingle)
+        b = array([[2+1j,1+2j]], dtype=csingle)
+        self.do(a, b)
+
+    def test_double(self):
+        a = array([[[1+2j, 2+3j], [3+4j, 4+5j]]], dtype=cdouble)
+        b = array([[2+1j,1+2j]], dtype=cdouble)
+        self.do(a,b)
+
+    def test_cdouble_2(self):
+        a = array([[[1+2j, 2+3j], [3+4j, 4+5j]]], dtype=cdouble)
+        b = array([[[2+1j, 1+2j, 1+3j], [1-2j, 1-3j, 1-6j]]], dtype=cdouble)
+        self.do(a,b)
+
+
+class HermitianTestCase(object):
+    def test_single(self):
+        a = array([[[1,2], [2,1]]], dtype=single)
+        b = array([[2.,1.]], dtype=single)
+        self.do(a,b)
+
+    def test_double(self):
+        a = array([[[1,2], [2,1]]], dtype=double)
+        b = array([[2.,1.]], dtype=double)
+        self.do(a,b)
+
+    def test_double_2(self):
+        a = array([[[1.,2.], [3.,4.]]], dtype=double)
+        b = array([[[2.,1., 4.], [3.,4.,6.]]], dtype=double)
+        self.do(a,b)
+
+    def test_csingle(self):
+        a = array([[[1,2+3j], [2-3j,1]]], dtype=csingle)
+        b = array([[2+1j,1+2j]], dtype=csingle)
+        self.do(a,b)
+
+    def test_cdouble(self):
+        a = array([[[1,2+3j], [2-3j,1]]], dtype=cdouble)
+        b = array([[2+1j,1+2j]], dtype=cdouble)
+        self.do(a,b)
+
+    def test_cdouble_2(self):
+        a = array([[[1,2+3j], [2-3j,1]]], dtype=cdouble)
+        b = array([[[2+1j, 1+2j, 1+3j], [1-2j, 1-3j, 1-6j]]], dtype=cdouble)
+        self.do(a,b)
+
+class TestInv(GeneralTestCase, TestCase):
+    def do(self, a, b):
+        a_inv = gula.inv(a)
+        assert_almost_equal(dot(a, a_inv), identity(a.shape[0]))
+
+
+class TestPoinv(HermitianTestCase, TestCase):
+    def do(self, a, b):
+        a_inv = gula.poinv(a)
+        assert_almost_equal(dot(a, a_inv), identity(a.shape[0]))
+
+
+class TestDet(GeneralTestCase, TestCase):
+    def do(self, a, b):
+        d = gula.det(a)
+        s, ld = gula.slogdet(a)
+        ev = gula.eigvals(a)
+        assert_almost_equal(d, multiply.reduce(ev, axis=(ev.ndim-1)))
+        assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=(ev.ndim-1)))
+        if s != 0:
+            assert_almost_equal(np.abs(s), 1)
+        else:
+            assert_equal(ld, -inf)
+
+    def test_zero(self):
+        assert_equal(gula.det([[0.0]]), 0.0)
+        assert_equal(gula.det([[0.0j]]), 0.0)
+        assert_equal(gula.slogdet([[0.0]]), (0.0, -inf))
+        assert_equal(gula.slogdet([[0-0j]]), (0.0j, -inf))
+
+
+class TestEig(GeneralTestCase, TestCase):
+    def do(self, a, b):
+        evalues, evectors = gula.eig(a)
+        assert_almost_equal(dot(a, evectors), multiply(evectors, evalues))
+
+
+class TestEigh(HermitianTestCase, TestCase):
+    def do(self, a, b):
+        evalues_lo, evectors_lo = gula.eigh(a,'L')
+        evalues_up, evectors_up = gula.eigh(a,'U')
+
+        assert_almost_equal(dot(a, evectors_lo), multiply(evectors_lo, evalues_lo))
+        assert_almost_equal(dot(a, evectors_up), multiply(evectors_up, evalues_up))
+        assert_almost_equal(evalues_lo, evalues_up)
+        assert_almost_equal(evectors_lo, evectors_up)
+
+
+class TestEigVals(GeneralTestCase, TestCase):
+    def do(self, a, b):
+        ev = gula.eigvals(a)
+        evalues, evectors = gula.eig(a)
+        assert_almost_equal(ev, evalues)
+
+
+class TestEigvalsh(HermitianTestCase, TestCase):
+    def do(self, a, b):
+        ev_lo = gula.eigvalsh(a, 'L')
+        ev_up = gula.eigvalsh(a, 'U')
+        evalues_lo, evectors_lo = gula.eigh(a, 'L')
+        evalues_up, evectors_up = gula.eigh(a, 'U')
+        assert_equal(ev_lo, evalues_lo)
+        assert_equal(ev_up, evalues_up)
+
+"""
+class TestSolve(GeneralTestCase,TestCase):
+    def do(self, a, b):
+        print a
+        print b
+        x = gula.solve(a,b)
+        assert_almost_equal(b,gula.matrix_multiply(a,x.T))
+"""
+
+class TestChosolve(HermitianTestCase, TestCase):
+    def do(self, a, b):
+        x_lo = gula.chosolve(a,b,'L')
+        x_up = gula.chosolve(a,b,'U')
+        assert_almost_equal(x_lo, x_up)
+        if a.dtype == single or a.dtype == double:
+            assert_almost_equal(b, gula.inner1d(a,x_lo))
+            assert_almost_equal(a, gula.inner1d(a,x_up))
+
+
+
+class TestSVD(GeneralTestCase, TestCase):
+    def do(self, a, b):
+        u, s, vt = gula.svd(a, 0)
+        assert_almost_equal(a, dot(multiply(u, s), vt))
+
+"""
+class TestCholesky(HermitianTestCase, TestCase):
+    def do(self, a, b):
+        pass
+"""
+
+if __name__ == "__main__":
+    run_module_suite()