from __future__ import division, print_function

import contextlib
import gc
from itertools import product, cycle
import sys
import warnings
from numbers import Number, Integral

import numpy as np

from numba import unittest_support as unittest
from numba import jit, errors
from numba.numpy_support import version as numpy_version
from .support import TestCase, tag
from .matmul_usecase import matmul_usecase, needs_matmul, needs_blas

try:
    import scipy.linalg.cython_lapack
    has_lapack = True
except ImportError:
    has_lapack = False

needs_lapack = unittest.skipUnless(has_lapack,
                                   "LAPACK needs Scipy 0.16+")


def dot2(a, b):
    return np.dot(a, b)


def dot3(a, b, out):
    return np.dot(a, b, out=out)


def vdot(a, b):
    return np.vdot(a, b)


class TestProduct(TestCase):
    """
    Tests for dot products.
    """

    dtypes = (np.float64, np.float32, np.complex128, np.complex64)

    def setUp(self):
        # Collect leftovers from previous test cases before checking for leaks
        gc.collect()

    def sample_vector(self, n, dtype):
        # Be careful to generate only exactly representable float values,
        # to avoid rounding discrepancies between Numpy and Numba
        base = np.arange(n)
        if issubclass(dtype, np.complexfloating):
            return (base * (1 - 0.5j) + 2j).astype(dtype)
        else:
            return (base * 0.5 - 1).astype(dtype)

    def sample_matrix(self, m, n, dtype):
        return self.sample_vector(m * n, dtype).reshape((m, n))

    @contextlib.contextmanager
    def check_contiguity_warning(self, pyfunc):
        """
        Check performance warning(s) for non-contiguity.
        """
        with warnings.catch_warnings(record=True) as w:
            warnings.simplefilter('always', errors.PerformanceWarning)
            yield
        self.assertGreaterEqual(len(w), 1)
        self.assertIs(w[0].category, errors.PerformanceWarning)
        self.assertIn("faster on contiguous arrays", str(w[0].message))
        self.assertEqual(w[0].filename, pyfunc.__code__.co_filename)
        # This works because our functions are one-liners
        self.assertEqual(w[0].lineno, pyfunc.__code__.co_firstlineno + 1)

    def check_func(self, pyfunc, cfunc, args):
        with self.assertNoNRTLeak():
            expected = pyfunc(*args)
            got = cfunc(*args)
            self.assertPreciseEqual(got, expected, ignore_sign_on_zero=True)
            del got, expected

    def check_func_out(self, pyfunc, cfunc, args, out):
        with self.assertNoNRTLeak():
            expected = np.copy(out)
            got = np.copy(out)
            self.assertIs(pyfunc(*args, out=expected), expected)
            self.assertIs(cfunc(*args, out=got), got)
            self.assertPreciseEqual(got, expected, ignore_sign_on_zero=True)
            del got, expected

    def assert_mismatching_sizes(self, cfunc, args, is_out=False):
        with self.assertRaises(ValueError) as raises:
            cfunc(*args)
        msg = ("incompatible output array size" if is_out else
               "incompatible array sizes")
        self.assertIn(msg, str(raises.exception))

    def assert_mismatching_dtypes(self, cfunc, args, func_name="np.dot()"):
        with self.assertRaises(errors.TypingError) as raises:
            cfunc(*args)
        self.assertIn("%s arguments must all have the same dtype"
                      % (func_name,),
                      str(raises.exception))

    @needs_blas
    def check_dot_vv(self, pyfunc, func_name):
        n = 3
        cfunc = jit(nopython=True)(pyfunc)
        for dtype in self.dtypes:
            a = self.sample_vector(n, dtype)
            b = self.sample_vector(n, dtype)
            self.check_func(pyfunc, cfunc, (a, b))
            # Non-contiguous
            self.check_func(pyfunc, cfunc, (a[::-1], b[::-1]))

        # Mismatching sizes
        a = self.sample_vector(n - 1, np.float64)
        b = self.sample_vector(n, np.float64)
        self.assert_mismatching_sizes(cfunc, (a, b))
        # Mismatching dtypes
        a = self.sample_vector(n, np.float32)
        b = self.sample_vector(n, np.float64)
        self.assert_mismatching_dtypes(cfunc, (a, b), func_name=func_name)

    def test_dot_vv(self):
        """
        Test vector * vector np.dot()
        """
        self.check_dot_vv(dot2, "np.dot()")

    def test_vdot(self):
        """
        Test np.vdot()
        """
        self.check_dot_vv(vdot, "np.vdot()")

    @needs_blas
    def check_dot_vm(self, pyfunc2, pyfunc3, func_name):
        m, n = 2, 3

        def samples(m, n):
            for order in 'CF':
                a = self.sample_matrix(m, n, np.float64).copy(order=order)
                b = self.sample_vector(n, np.float64)
                yield a, b
            for dtype in self.dtypes:
                a = self.sample_matrix(m, n, dtype)
                b = self.sample_vector(n, dtype)
                yield a, b
            # Non-contiguous
            yield a[::-1], b[::-1]

        cfunc2 = jit(nopython=True)(pyfunc2)
        if pyfunc3 is not None:
            cfunc3 = jit(nopython=True)(pyfunc3)
        for a, b in samples(m, n):
            self.check_func(pyfunc2, cfunc2, (a, b))
            self.check_func(pyfunc2, cfunc2, (b, a.T))
        if pyfunc3 is not None:
            for a, b in samples(m, n):
                out = np.empty(m, dtype=a.dtype)
                self.check_func_out(pyfunc3, cfunc3, (a, b), out)
                self.check_func_out(pyfunc3, cfunc3, (b, a.T), out)

        # Mismatching sizes
        a = self.sample_matrix(m, n - 1, np.float64)
        b = self.sample_vector(n, np.float64)
        self.assert_mismatching_sizes(cfunc2, (a, b))
        self.assert_mismatching_sizes(cfunc2, (b, a.T))
        if pyfunc3 is not None:
            out = np.empty(m, np.float64)
            self.assert_mismatching_sizes(cfunc3, (a, b, out))
            self.assert_mismatching_sizes(cfunc3, (b, a.T, out))
            a = self.sample_matrix(m, m, np.float64)
            b = self.sample_vector(m, np.float64)
            out = np.empty(m - 1, np.float64)
            self.assert_mismatching_sizes(cfunc3, (a, b, out), is_out=True)
            self.assert_mismatching_sizes(cfunc3, (b, a.T, out), is_out=True)
        # Mismatching dtypes
        a = self.sample_matrix(m, n, np.float32)
        b = self.sample_vector(n, np.float64)
        self.assert_mismatching_dtypes(cfunc2, (a, b), func_name)
        if pyfunc3 is not None:
            a = self.sample_matrix(m, n, np.float64)
            b = self.sample_vector(n, np.float64)
            out = np.empty(m, np.float32)
            self.assert_mismatching_dtypes(cfunc3, (a, b, out), func_name)

    def test_dot_vm(self):
        """
        Test vector * matrix and matrix * vector np.dot()
        """
        self.check_dot_vm(dot2, dot3, "np.dot()")

    @needs_blas
    def check_dot_mm(self, pyfunc2, pyfunc3, func_name):

        def samples(m, n, k):
            for order_a, order_b in product('CF', 'CF'):
                a = self.sample_matrix(m, k, np.float64).copy(order=order_a)
                b = self.sample_matrix(k, n, np.float64).copy(order=order_b)
                yield a, b
            for dtype in self.dtypes:
                a = self.sample_matrix(m, k, dtype)
                b = self.sample_matrix(k, n, dtype)
                yield a, b
            # Non-contiguous
            yield a[::-1], b[::-1]

        cfunc2 = jit(nopython=True)(pyfunc2)
        if pyfunc3 is not None:
            cfunc3 = jit(nopython=True)(pyfunc3)

        # Test generic matrix * matrix as well as "degenerate" cases
        # where one of the outer dimensions is 1 (i.e. really represents
        # a vector, which may select a different implementation)
        for m, n, k in [(2, 3, 4),  # Generic matrix * matrix
                        (1, 3, 4),  # 2d vector * matrix
                        (1, 1, 4),  # 2d vector * 2d vector
                        ]:
            for a, b in samples(m, n, k):
                self.check_func(pyfunc2, cfunc2, (a, b))
                self.check_func(pyfunc2, cfunc2, (b.T, a.T))
            if pyfunc3 is not None:
                for a, b in samples(m, n, k):
                    out = np.empty((m, n), dtype=a.dtype)
                    self.check_func_out(pyfunc3, cfunc3, (a, b), out)
                    out = np.empty((n, m), dtype=a.dtype)
                    self.check_func_out(pyfunc3, cfunc3, (b.T, a.T), out)

        # Mismatching sizes
        m, n, k = 2, 3, 4
        a = self.sample_matrix(m, k - 1, np.float64)
        b = self.sample_matrix(k, n, np.float64)
        self.assert_mismatching_sizes(cfunc2, (a, b))
        if pyfunc3 is not None:
            out = np.empty((m, n), np.float64)
            self.assert_mismatching_sizes(cfunc3, (a, b, out))
            a = self.sample_matrix(m, k, np.float64)
            b = self.sample_matrix(k, n, np.float64)
            out = np.empty((m, n - 1), np.float64)
            self.assert_mismatching_sizes(cfunc3, (a, b, out), is_out=True)
        # Mismatching dtypes
        a = self.sample_matrix(m, k, np.float32)
        b = self.sample_matrix(k, n, np.float64)
        self.assert_mismatching_dtypes(cfunc2, (a, b), func_name)
        if pyfunc3 is not None:
            a = self.sample_matrix(m, k, np.float64)
            b = self.sample_matrix(k, n, np.float64)
            out = np.empty((m, n), np.float32)
            self.assert_mismatching_dtypes(cfunc3, (a, b, out), func_name)

    @tag('important')
    def test_dot_mm(self):
        """
        Test matrix * matrix np.dot()
        """
        self.check_dot_mm(dot2, dot3, "np.dot()")

    @needs_matmul
    def test_matmul_vv(self):
        """
        Test vector @ vector
        """
        self.check_dot_vv(matmul_usecase, "'@'")

    @needs_matmul
    def test_matmul_vm(self):
        """
        Test vector @ matrix and matrix @ vector
        """
        self.check_dot_vm(matmul_usecase, None, "'@'")

    @needs_matmul
    def test_matmul_mm(self):
        """
        Test matrix @ matrix
        """
        self.check_dot_mm(matmul_usecase, None, "'@'")

    @needs_blas
    def test_contiguity_warnings(self):
        m, k, n = 2, 3, 4
        dtype = np.float64
        a = self.sample_matrix(m, k, dtype)[::-1]
        b = self.sample_matrix(k, n, dtype)[::-1]
        out = np.empty((m, n), dtype)

        cfunc = jit(nopython=True)(dot2)
        with self.check_contiguity_warning(cfunc.py_func):
            cfunc(a, b)
        cfunc = jit(nopython=True)(dot3)
        with self.check_contiguity_warning(cfunc.py_func):
            cfunc(a, b, out)

        a = self.sample_vector(n, dtype)[::-1]
        b = self.sample_vector(n, dtype)[::-1]

        cfunc = jit(nopython=True)(vdot)
        with self.check_contiguity_warning(cfunc.py_func):
            cfunc(a, b)


# Implementation definitions for the purpose of jitting.

def invert_matrix(a):
    return np.linalg.inv(a)


def cholesky_matrix(a):
    return np.linalg.cholesky(a)


def eig_matrix(a):
    return np.linalg.eig(a)


def eigvals_matrix(a):
    return np.linalg.eigvals(a)


def eigh_matrix(a):
    return np.linalg.eigh(a)


def eigvalsh_matrix(a):
    return np.linalg.eigvalsh(a)


def svd_matrix(a, full_matrices=1):
    return np.linalg.svd(a, full_matrices)


def qr_matrix(a):
    return np.linalg.qr(a)


def lstsq_system(A, B, rcond=-1):
    return np.linalg.lstsq(A, B, rcond)


def solve_system(A, B):
    return np.linalg.solve(A, B)


def pinv_matrix(A, rcond=1e-15):  # 1e-15 from numpy impl
    return np.linalg.pinv(A)


def slogdet_matrix(a):
    return np.linalg.slogdet(a)


def det_matrix(a):
    return np.linalg.det(a)


def norm_matrix(a, ord=None):
    return np.linalg.norm(a, ord)


def cond_matrix(a, p=None):
    return np.linalg.cond(a, p)


def matrix_rank_matrix(a, tol=None):
    return np.linalg.matrix_rank(a, tol)


def matrix_power_matrix(a, n):
    return np.linalg.matrix_power(a, n)


def trace_matrix(a, offset=0):
    return np.trace(a, offset)


def trace_matrix_no_offset(a):
    return np.trace(a)


if numpy_version >= (1, 9):
    def outer_matrix(a, b, out=None):
        return np.outer(a, b, out=out)
else:
    def outer_matrix(a, b):
        return np.outer(a, b)


def kron_matrix(a, b):
    return np.kron(a, b)


class TestLinalgBase(TestCase):
    """
    Provides setUp and common data/error modes for testing np.linalg functions.
    """

    # supported dtypes
    dtypes = (np.float64, np.float32, np.complex128, np.complex64)

    def setUp(self):
        # Collect leftovers from previous test cases before checking for leaks
        gc.collect()

    def sample_vector(self, n, dtype):
        # Be careful to generate only exactly representable float values,
        # to avoid rounding discrepancies between Numpy and Numba
        base = np.arange(n)
        if issubclass(dtype, np.complexfloating):
            return (base * (1 - 0.5j) + 2j).astype(dtype)
        else:
            return (base * 0.5 + 1).astype(dtype)

    def specific_sample_matrix(
            self, size, dtype, order, rank=None, condition=None):
        """
        Provides a sample matrix with an optionally specified rank or condition
        number.

        size: (rows, columns), the dimensions of the returned matrix.
        dtype: the dtype for the returned matrix.
        order: the memory layout for the returned matrix, 'F' or 'C'.
        rank: the rank of the matrix, an integer value, defaults to full rank.
        condition: the condition number of the matrix (defaults to 1.)

        NOTE: Only one of rank or condition may be set.
        """

        # default condition
        d_cond = 1.

        if len(size) != 2:
            raise ValueError("size must be a length 2 tuple.")

        if order not in ['F', 'C']:
            raise ValueError("order must be one of 'F' or 'C'.")

        if dtype not in [np.float32, np.float64, np.complex64, np.complex128]:
            raise ValueError("dtype must be a numpy floating point type.")

        if rank is not None and condition is not None:
            raise ValueError("Only one of rank or condition can be specified.")

        if condition is None:
            condition = d_cond

        if condition < 1:
            raise ValueError("Condition number must be >=1.")

        np.random.seed(0)  # repeatable seed
        m, n = size

        if m < 0 or n < 0:
            raise ValueError("Negative dimensions given for matrix shape.")

        minmn = min(m, n)
        if rank is None:
            rv = minmn
        else:
            if rank <= 0:
                raise ValueError("Rank must be greater than zero.")
            if not isinstance(rank, Integral):
                raise ValueError("Rank must an integer.")
            rv = rank
            if rank > minmn:
                raise ValueError("Rank given greater than full rank.")

        if m == 1 or n == 1:
            # vector, must be rank 1 (enforced above)
            # condition of vector is also 1
            if condition != d_cond:
                raise ValueError(
                    "Condition number was specified for a vector (always 1.).")
            maxmn = max(m, n)
            Q = self.sample_vector(maxmn, dtype).reshape(m, n)
        else:
            # Build a sample matrix via combining SVD like inputs.

            # Create matrices of left and right singular vectors.
            # This could use Modified Gram-Schmidt and perhaps be quicker,
            # at present it uses QR decompositions to obtain orthonormal
            # matrices.
            tmp = self.sample_vector(m * m, dtype).reshape(m, m)
            U, _ = np.linalg.qr(tmp)
            # flip the second array, else for m==n the identity matrix appears
            tmp = self.sample_vector(n * n, dtype)[::-1].reshape(n, n)
            V, _ = np.linalg.qr(tmp)
            # create singular values.
            sv = np.linspace(d_cond, condition, rv)
            S = np.zeros((m, n))
            idx = np.nonzero(np.eye(m, n))
            S[idx[0][:rv], idx[1][:rv]] = sv
            Q = np.dot(np.dot(U, S), V.T)  # construct
            Q = np.array(Q, dtype=dtype, order=order)  # sort out order/type

        return Q

    def assert_error(self, cfunc, args, msg, err=ValueError):
        with self.assertRaises(err) as raises:
            cfunc(*args)
        self.assertIn(msg, str(raises.exception))

    def assert_non_square(self, cfunc, args):
        msg = "Last 2 dimensions of the array must be square."
        self.assert_error(cfunc, args, msg, np.linalg.LinAlgError)

    def assert_wrong_dtype(self, name, cfunc, args):
        msg = "np.linalg.%s() only supported on float and complex arrays" % name
        self.assert_error(cfunc, args, msg, errors.TypingError)

    def assert_wrong_dimensions(self, name, cfunc, args, la_prefix=True):
        prefix = "np.linalg" if la_prefix else "np"
        msg = "%s.%s() only supported on 2-D arrays" % (prefix, name)
        self.assert_error(cfunc, args, msg, errors.TypingError)

    def assert_no_nan_or_inf(self, cfunc, args):
        msg = "Array must not contain infs or NaNs."
        self.assert_error(cfunc, args, msg, np.linalg.LinAlgError)

    def assert_contig_sanity(self, got, expected_contig):
        """
        This checks that in a computed result from numba (array, possibly tuple
        of arrays) all the arrays are contiguous in memory and that they are
        all at least one of "C_CONTIGUOUS" or "F_CONTIGUOUS". The computed
        result of the contiguousness is then compared against a hardcoded
        expected result.

        got: is the computed results from numba
        expected_contig: is "C" or "F" and is the expected type of
                        contiguousness across all input values
                        (and therefore tests).
        """

        if isinstance(got, tuple):
            # tuple present, check all results
            for a in got:
                self.assert_contig_sanity(a, expected_contig)
        else:
            if not isinstance(got, Number):
                # else a single array is present
                c_contig = got.flags.c_contiguous
                f_contig = got.flags.f_contiguous

                # check that the result (possible set of) is at least one of
                # C or F contiguous.
                msg = "Results are not at least one of all C or F contiguous."
                self.assertTrue(c_contig | f_contig, msg)

                msg = "Computed contiguousness does not match expected."
                if expected_contig == "C":
                    self.assertTrue(c_contig, msg)
                elif expected_contig == "F":
                    self.assertTrue(f_contig, msg)
                else:
                    raise ValueError("Unknown contig")

    def assert_raise_on_singular(self, cfunc, args):
        msg = "Matrix is singular to machine precision."
        self.assert_error(cfunc, args, msg, err=np.linalg.LinAlgError)

    def assert_is_identity_matrix(self, got, rtol=None, atol=None):
        """
        Checks if a matrix is equal to the identity matrix.
        """
        # check it is square
        self.assertEqual(got.shape[-1], got.shape[-2])
        # create identity matrix
        eye = np.eye(got.shape[-1], dtype=got.dtype)
        resolution = 5 * np.finfo(got.dtype).resolution
        if rtol is None:
            rtol = 10 * resolution
        if atol is None:
            atol = 100 * resolution  # zeros tend to be fuzzy
        # check it matches
        np.testing.assert_allclose(got, eye, rtol, atol)

    def assert_invalid_norm_kind(self, cfunc, args):
        """
        For use in norm() and cond() tests.
        """
        msg = "Invalid norm order for matrices."
        self.assert_error(cfunc, args, msg, ValueError)


class TestTestLinalgBase(TestCase):
    """
    The sample matrix code TestLinalgBase.specific_sample_matrix()
    is a bit involved, this class tests it works as intended.
    """

    def test_specific_sample_matrix(self):

        # add a default test to the ctor, it never runs so doesn't matter
        inst = TestLinalgBase('specific_sample_matrix')

        sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)]

        # test loop
        for size, dtype, order in product(sizes, inst.dtypes, 'FC'):

            m, n = size
            minmn = min(m, n)

            # test default full rank
            A = inst.specific_sample_matrix(size, dtype, order)
            self.assertEqual(A.shape, size)
            self.assertEqual(np.linalg.matrix_rank(A), minmn)

            # test reduced rank if a reduction is possible
            if minmn > 1:
                rank = minmn - 1
                A = inst.specific_sample_matrix(size, dtype, order, rank=rank)
                self.assertEqual(A.shape, size)
                self.assertEqual(np.linalg.matrix_rank(A), rank)

            resolution = 5 * np.finfo(dtype).resolution

            # test default condition
            A = inst.specific_sample_matrix(size, dtype, order)
            self.assertEqual(A.shape, size)
            np.testing.assert_allclose(np.linalg.cond(A),
                                       1.,
                                       rtol=resolution,
                                       atol=resolution)

            # test specified condition if matrix is > 1D
            if minmn > 1:
                condition = 10.
                A = inst.specific_sample_matrix(
                    size, dtype, order, condition=condition)
                self.assertEqual(A.shape, size)
                np.testing.assert_allclose(np.linalg.cond(A),
                                           10.,
                                           rtol=resolution,
                                           atol=resolution)

        # check errors are raised appropriately
        def check_error(args, msg, err=ValueError):
            with self.assertRaises(err) as raises:
                inst.specific_sample_matrix(*args)
            self.assertIn(msg, str(raises.exception))

        # check the checker runs ok
        with self.assertRaises(AssertionError) as raises:
            msg = "blank"
            check_error(((2, 3), np.float64, 'F'), msg, err=ValueError)

        # check invalid inputs...

        # bad size
        msg = "size must be a length 2 tuple."
        check_error(((1,), np.float64, 'F'), msg, err=ValueError)

        # bad order
        msg = "order must be one of 'F' or 'C'."
        check_error(((2, 3), np.float64, 'z'), msg, err=ValueError)

        # bad type
        msg = "dtype must be a numpy floating point type."
        check_error(((2, 3), np.int32, 'F'), msg, err=ValueError)

        # specifying both rank and condition
        msg = "Only one of rank or condition can be specified."
        check_error(((2, 3), np.float64, 'F', 1, 1), msg, err=ValueError)

        # specifying negative condition
        msg = "Condition number must be >=1."
        check_error(((2, 3), np.float64, 'F', None, -1), msg, err=ValueError)

        # specifying negative matrix dimension
        msg = "Negative dimensions given for matrix shape."
        check_error(((2, -3), np.float64, 'F'), msg, err=ValueError)

        # specifying negative rank
        msg = "Rank must be greater than zero."
        check_error(((2, 3), np.float64, 'F', -1), msg, err=ValueError)

        # specifying a rank greater than maximum rank
        msg = "Rank given greater than full rank."
        check_error(((2, 3), np.float64, 'F', 4), msg, err=ValueError)

        # specifying a condition number for a vector
        msg = "Condition number was specified for a vector (always 1.)."
        check_error(((1, 3), np.float64, 'F', None, 10), msg, err=ValueError)

        # specifying a non integer rank
        msg = "Rank must an integer."
        check_error(((2, 3), np.float64, 'F', 1.5), msg, err=ValueError)


class TestLinalgInv(TestLinalgBase):
    """
    Tests for np.linalg.inv.
    """

    @tag('important')
    @needs_lapack
    def test_linalg_inv(self):
        """
        Test np.linalg.inv
        """
        n = 10
        cfunc = jit(nopython=True)(invert_matrix)

        def check(a, **kwargs):
            expected = invert_matrix(a)
            got = cfunc(a)
            self.assert_contig_sanity(got, "F")

            use_reconstruction = False

            # try strict
            try:
                np.testing.assert_array_almost_equal_nulp(got, expected,
                                                          nulp=10)
            except AssertionError:
                # fall back to reconstruction
                use_reconstruction = True

            if use_reconstruction:
                rec = np.dot(got, a)
                self.assert_is_identity_matrix(rec)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a)

        for dtype, order in product(self.dtypes, 'CF'):
            a = self.specific_sample_matrix((n, n), dtype, order)
            check(a)

        # Non square matrix
        self.assert_non_square(cfunc, (np.ones((2, 3)),))

        # Wrong dtype
        self.assert_wrong_dtype("inv", cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions("inv", cfunc, (np.ones(10),))

        # Singular matrix
        self.assert_raise_on_singular(cfunc, (np.zeros((2, 2)),))


class TestLinalgCholesky(TestLinalgBase):
    """
    Tests for np.linalg.cholesky.
    """

    def sample_matrix(self, m, dtype, order):
        # pd. (positive definite) matrix has eigenvalues in Z+
        np.random.seed(0)  # repeatable seed
        A = np.random.rand(m, m)
        # orthonormal q needed to form up q^{-1}*D*q
        # no "orth()" in numpy
        q, _ = np.linalg.qr(A)
        L = np.arange(1, m + 1)  # some positive eigenvalues
        Q = np.dot(np.dot(q.T, np.diag(L)), q)  # construct
        Q = np.array(Q, dtype=dtype, order=order)  # sort out order/type
        return Q

    def assert_not_pd(self, cfunc, args):
        msg = "Matrix is not positive definite."
        self.assert_error(cfunc, args, msg, np.linalg.LinAlgError)

    @needs_lapack
    def test_linalg_cholesky(self):
        """
        Test np.linalg.cholesky
        """
        n = 10
        cfunc = jit(nopython=True)(cholesky_matrix)

        def check(a):
            expected = cholesky_matrix(a)
            got = cfunc(a)
            use_reconstruction = False
            # check that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "C")

            # try strict
            try:
                np.testing.assert_array_almost_equal_nulp(got, expected,
                                                          nulp=10)
            except AssertionError:
                # fall back to reconstruction
                use_reconstruction = True

            # try via reconstruction
            if use_reconstruction:
                rec = np.dot(got, np.conj(got.T))
                resolution = 5 * np.finfo(a.dtype).resolution
                np.testing.assert_allclose(
                    a,
                    rec,
                    rtol=resolution,
                    atol=resolution
                )

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a)

        for dtype, order in product(self.dtypes, 'FC'):
            a = self.sample_matrix(n, dtype, order)
            check(a)

        rn = "cholesky"
        # Non square matrices
        self.assert_non_square(cfunc, (np.ones((2, 3), dtype=np.float64),))

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions(rn, cfunc,
                                     (np.ones(10, dtype=np.float64),))

        # not pd
        self.assert_not_pd(cfunc,
                           (np.ones(4, dtype=np.float64).reshape(2, 2),))


class TestLinalgEigenSystems(TestLinalgBase):
    """
    Tests for np.linalg.eig/eigvals.
    """

    def sample_matrix(self, m, dtype, order):
        # This is a tridiag with the same but skewed values on the diagonals
        v = self.sample_vector(m, dtype)
        Q = np.diag(v)
        idx = np.nonzero(np.eye(Q.shape[0], Q.shape[1], 1))
        Q[idx] = v[1:]
        idx = np.nonzero(np.eye(Q.shape[0], Q.shape[1], -1))
        Q[idx] = v[:-1]
        Q = np.array(Q, dtype=dtype, order=order)
        return Q

    def assert_no_domain_change(self, name, cfunc, args):
        msg = name + "() argument must not cause a domain change."
        self.assert_error(cfunc, args, msg)

    def checker_for_linalg_eig(
            self, name, func, expected_res_len, check_for_domain_change=None):
        """
        Test np.linalg.eig
        """
        n = 10
        cfunc = jit(nopython=True)(func)

        def check(a):
            expected = func(a)
            got = cfunc(a)
            # check that the returned tuple is same length
            self.assertEqual(len(expected), len(got))
            # and that dimension is correct
            res_is_tuple = False
            if isinstance(got, tuple):
                res_is_tuple = True
                self.assertEqual(len(got), expected_res_len)
            else:  # its an array
                self.assertEqual(got.ndim, expected_res_len)

            # and that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "F")

            use_reconstruction = False
            # try plain match of each array to np first
            for k in range(len(expected)):
                try:
                    np.testing.assert_array_almost_equal_nulp(
                        got[k], expected[k], nulp=10)
                except AssertionError:
                    # plain match failed, test by reconstruction
                    use_reconstruction = True

            # If plain match fails then reconstruction is used.
            # this checks that A*V ~== V*diag(W)
            # i.e. eigensystem ties out
            # this is required as numpy uses only double precision lapack
            # routines and computation of eigenvectors is numerically
            # sensitive, numba uses the type specific routines therefore
            # sometimes comes out with a different (but entirely
            # valid) answer (eigenvectors are not unique etc.).
            # This is only applicable if eigenvectors are computed
            # along with eigenvalues i.e. result is a tuple.
            resolution = 5 * np.finfo(a.dtype).resolution
            if use_reconstruction:
                if res_is_tuple:
                    w, v = got
                    # modify 'a' if hermitian eigensystem functionality is
                    # being tested. 'L' for use lower part is default and
                    # the only thing used at present so we conjugate transpose
                    # the lower part into the upper for use in the
                    # reconstruction. By construction the sample matrix is
                    # tridiag so this is just a question of copying the lower
                    # diagonal into the upper and conjugating on the way.
                    if name[-1] == 'h':
                        idxl = np.nonzero(np.eye(a.shape[0], a.shape[1], -1))
                        idxu = np.nonzero(np.eye(a.shape[0], a.shape[1], 1))
                        cfunc(a)
                        # upper idx must match lower for default uplo="L"
                        # if complex, conjugate
                        a[idxu] = np.conj(a[idxl])
                        # also, only the real part of the diagonals is
                        # considered in the calculation so the imag is zeroed
                        # out for the purposes of use in reconstruction.
                        a[np.diag_indices(a.shape[0])] = np.real(np.diag(a))

                    lhs = np.dot(a, v)
                    rhs = np.dot(v, np.diag(w))

                    np.testing.assert_allclose(
                        lhs.real,
                        rhs.real,
                        rtol=resolution,
                        atol=resolution
                    )
                    if np.iscomplexobj(v):
                        np.testing.assert_allclose(
                            lhs.imag,
                            rhs.imag,
                            rtol=resolution,
                            atol=resolution
                        )
                else:
                    # This isn't technically reconstruction but is here to
                    # deal with that the order of the returned eigenvalues
                    # may differ in the case of routines just returning
                    # eigenvalues and there's no true reconstruction
                    # available with which to perform a check.
                    np.testing.assert_allclose(
                        np.sort(expected),
                        np.sort(got),
                        rtol=resolution,
                        atol=resolution
                    )

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a)

        # The main test loop
        for dtype, order in product(self.dtypes, 'FC'):
            a = self.sample_matrix(n, dtype, order)
            check(a)

        # Test both a real and complex type as the impls are different
        for ty in [np.float32, np.complex64]:
            # Non square matrices
            self.assert_non_square(cfunc, (np.ones((2, 3), dtype=ty),))

            # Wrong dtype
            self.assert_wrong_dtype(name, cfunc,
                                    (np.ones((2, 2), dtype=np.int32),))

            # Dimension issue
            self.assert_wrong_dimensions(name, cfunc, (np.ones(10, dtype=ty),))

            # no nans or infs
            self.assert_no_nan_or_inf(cfunc,
                                      (np.array([[1., 2., ], [np.inf, np.nan]],
                                                dtype=ty),))

        if check_for_domain_change:
            # By design numba does not support dynamic return types, numpy does
            # and uses this in the case of returning eigenvalues/vectors of
            # a real matrix. The return type of np.linalg.eig(), when
            # operating on a matrix in real space depends on the values present
            # in the matrix itself (recalling that eigenvalues are the roots of the
            # characteristic polynomial of the system matrix, which will by
            # construction depend on the values present in the system matrix).
            # This test asserts that if a domain change is required on the return
            # type, i.e. complex eigenvalues from a real input, an error is raised.
            # For complex types, regardless of the value of the imaginary part of
            # the returned eigenvalues, a complex type will be returned, this
            # follows numpy and fits in with numba.

            # First check that the computation is valid (i.e. in complex space)
            A = np.array([[1, -2], [2, 1]])
            check(A.astype(np.complex128))
            # and that the imaginary part is nonzero
            l, _ = func(A)
            self.assertTrue(np.any(l.imag))

            # Now check that the computation fails in real space
            for ty in [np.float32, np.float64]:
                self.assert_no_domain_change(name, cfunc, (A.astype(ty),))

    @needs_lapack
    def test_linalg_eig(self):
        self.checker_for_linalg_eig("eig", eig_matrix, 2, True)

    @needs_lapack
    def test_linalg_eigvals(self):
        self.checker_for_linalg_eig("eigvals", eigvals_matrix, 1, True)

    @needs_lapack
    def test_linalg_eigh(self):
        self.checker_for_linalg_eig("eigh", eigh_matrix, 2, False)

    @needs_lapack
    def test_linalg_eigvalsh(self):
        self.checker_for_linalg_eig("eigvalsh", eigvalsh_matrix, 1, False)


class TestLinalgSvd(TestLinalgBase):
    """
    Tests for np.linalg.svd.
    """

    @needs_lapack
    def test_linalg_svd(self):
        """
        Test np.linalg.svd
        """
        cfunc = jit(nopython=True)(svd_matrix)

        def check(a, **kwargs):
            expected = svd_matrix(a, **kwargs)
            got = cfunc(a, **kwargs)
            # check that the returned tuple is same length
            self.assertEqual(len(expected), len(got))
            # and that length is 3
            self.assertEqual(len(got), 3)
            # and that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "F")

            use_reconstruction = False
            # try plain match of each array to np first
            for k in range(len(expected)):

                try:
                    np.testing.assert_array_almost_equal_nulp(
                        got[k], expected[k], nulp=10)
                except AssertionError:
                    # plain match failed, test by reconstruction
                    use_reconstruction = True

            # if plain match fails then reconstruction is used.
            # this checks that A ~= U*S*V**H
            # i.e. SV decomposition ties out
            # this is required as numpy uses only double precision lapack
            # routines and computation of svd is numerically
            # sensitive, numba using the type specific routines therefore
            # sometimes comes out with a different answer (orthonormal bases
            # are not unique etc.).
            if use_reconstruction:
                u, sv, vt = got

                # check they are dimensionally correct
                for k in range(len(expected)):
                    self.assertEqual(got[k].shape, expected[k].shape)

                # regardless of full_matrices cols in u and rows in vt
                # dictates the working size of s
                s = np.zeros((u.shape[1], vt.shape[0]))
                np.fill_diagonal(s, sv)

                rec = np.dot(np.dot(u, s), vt)
                resolution = np.finfo(a.dtype).resolution
                np.testing.assert_allclose(
                    a,
                    rec,
                    rtol=10 * resolution,
                    atol=100 * resolution  # zeros tend to be fuzzy
                )

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, **kwargs)

        # test: column vector, tall, wide, square, row vector
        # prime sizes
        sizes = [(7, 1), (7, 5), (5, 7), (3, 3), (1, 7)]

        # flip on reduced or full matrices
        full_matrices = (True, False)

        # test loop
        for size, dtype, fmat, order in \
                product(sizes, self.dtypes, full_matrices, 'FC'):

            a = self.specific_sample_matrix(size, dtype, order)
            check(a, full_matrices=fmat)

        rn = "svd"

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions(rn, cfunc,
                                     (np.ones(10, dtype=np.float64),))

        # no nans or infs
        self.assert_no_nan_or_inf(cfunc,
                                  (np.array([[1., 2., ], [np.inf, np.nan]],
                                            dtype=np.float64),))


class TestLinalgQr(TestLinalgBase):
    """
    Tests for np.linalg.qr.
    """

    @needs_lapack
    def test_linalg_qr(self):
        """
        Test np.linalg.qr
        """
        cfunc = jit(nopython=True)(qr_matrix)

        def check(a, **kwargs):
            expected = qr_matrix(a, **kwargs)
            got = cfunc(a, **kwargs)

            # check that the returned tuple is same length
            self.assertEqual(len(expected), len(got))
            # and that length is 2
            self.assertEqual(len(got), 2)
            # and that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "F")

            use_reconstruction = False
            # try plain match of each array to np first
            for k in range(len(expected)):
                try:
                    np.testing.assert_array_almost_equal_nulp(
                        got[k], expected[k], nulp=10)
                except AssertionError:
                    # plain match failed, test by reconstruction
                    use_reconstruction = True

            # if plain match fails then reconstruction is used.
            # this checks that A ~= Q*R and that (Q^H)*Q = I
            # i.e. QR decomposition ties out
            # this is required as numpy uses only double precision lapack
            # routines and computation of qr is numerically
            # sensitive, numba using the type specific routines therefore
            # sometimes comes out with a different answer (orthonormal bases
            # are not unique etc.).
            if use_reconstruction:
                q, r = got

                # check they are dimensionally correct
                for k in range(len(expected)):
                    self.assertEqual(got[k].shape, expected[k].shape)

                # check A=q*r
                rec = np.dot(q, r)
                resolution = np.finfo(a.dtype).resolution
                np.testing.assert_allclose(
                    a,
                    rec,
                    rtol=10 * resolution,
                    atol=100 * resolution  # zeros tend to be fuzzy
                )

                # check q is orthonormal
                self.assert_is_identity_matrix(np.dot(np.conjugate(q.T), q))

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, **kwargs)

        # test: column vector, tall, wide, square, row vector
        # prime sizes
        sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)]

        # test loop
        for size, dtype, order in \
                product(sizes, self.dtypes, 'FC'):
            a = self.specific_sample_matrix(size, dtype, order)
            check(a)

        rn = "qr"

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions(rn, cfunc,
                                     (np.ones(10, dtype=np.float64),))

        # no nans or infs
        self.assert_no_nan_or_inf(cfunc,
                                  (np.array([[1., 2., ], [np.inf, np.nan]],
                                            dtype=np.float64),))


class TestLinalgSystems(TestLinalgBase):
    """
    Base class for testing "system" solvers from np.linalg.
    Namely np.linalg.solve() and np.linalg.lstsq().
    """

    # check for RHS with dimension > 2 raises
    def assert_wrong_dimensions_1D(self, name, cfunc, args, la_prefix=True):
        prefix = "np.linalg" if la_prefix else "np"
        msg = "%s.%s() only supported on 1 and 2-D arrays" % (prefix, name)
        self.assert_error(cfunc, args, msg, errors.TypingError)

    # check that a dimensionally invalid system raises
    def assert_dimensionally_invalid(self, cfunc, args):
        msg = "Incompatible array sizes, system is not dimensionally valid."
        self.assert_error(cfunc, args, msg, np.linalg.LinAlgError)

    # check that args with differing dtypes raise
    def assert_homogeneous_dtypes(self, name, cfunc, args):
        msg = "np.linalg.%s() only supports inputs that have homogeneous dtypes." % name
        self.assert_error(cfunc, args, msg, errors.TypingError)


class TestLinalgLstsq(TestLinalgSystems):
    """
    Tests for np.linalg.lstsq.
    """

    # NOTE: The testing of this routine is hard as it has to handle numpy
    # using double precision routines on single precision input, this has
    # a knock on effect especially in rank deficient cases and cases where
    # conditioning is generally poor. As a result computed ranks can differ
    # and consequently the calculated residual can differ.
    # The tests try and deal with this as best as they can through the use
    # of reconstruction and measures like residual norms.
    # Suggestions for improvements are welcomed!

    @needs_lapack
    def test_linalg_lstsq(self):
        """
        Test np.linalg.lstsq
        """
        cfunc = jit(nopython=True)(lstsq_system)

        def check(A, B, **kwargs):
            expected = lstsq_system(A, B, **kwargs)
            got = cfunc(A, B, **kwargs)

            # check that the returned tuple is same length
            self.assertEqual(len(expected), len(got))
            # and that length is 4
            self.assertEqual(len(got), 4)
            # and that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "C")

            use_reconstruction = False

            # check the ranks are the same and continue to a standard
            # match if that is the case (if ranks differ, then output
            # in e.g. residual array is of different size!).
            try:
                self.assertEqual(got[2], expected[2])
                # try plain match of each array to np first
                for k in range(len(expected)):
                    try:
                        np.testing.assert_array_almost_equal_nulp(
                            got[k], expected[k], nulp=10)
                    except AssertionError:
                        # plain match failed, test by reconstruction
                        use_reconstruction = True
            except AssertionError:
                use_reconstruction = True

            if use_reconstruction:
                x, res, rank, s = got

                # indicies in the output which are ndarrays
                out_array_idx = [0, 1, 3]

                try:
                    # check the ranks are the same
                    self.assertEqual(rank, expected[2])
                    # check they are dimensionally correct, skip [2] = rank.
                    for k in out_array_idx:
                        if isinstance(expected[k], np.ndarray):
                            self.assertEqual(got[k].shape, expected[k].shape)
                except AssertionError:
                    # check the rank differs by 1. (numerical fuzz)
                    self.assertTrue(abs(rank - expected[2]) < 2)

                # check if A*X = B
                resolution = np.finfo(A.dtype).resolution
                try:
                    # this will work so long as the conditioning is
                    # ok and the rank is full
                    rec = np.dot(A, x)
                    np.testing.assert_allclose(
                        B,
                        rec,
                        rtol=10 * resolution,
                        atol=10 * resolution
                    )
                except AssertionError:
                    # system is probably under/over determined and/or
                    # poorly conditioned. Check slackened equality
                    # and that the residual norm is the same.
                    for k in out_array_idx:
                        try:
                            np.testing.assert_allclose(
                                expected[k],
                                got[k],
                                rtol=100 * resolution,
                                atol=100 * resolution
                            )
                        except AssertionError:
                            # check the fail is likely due to bad conditioning
                            c = np.linalg.cond(A)
                            self.assertGreater(10 * c, (1. / resolution))

                        # make sure the residual 2-norm is ok
                        # if this fails its probably due to numpy using double
                        # precision LAPACK routines for singles.
                        res_expected = np.linalg.norm(
                            B - np.dot(A, expected[0]))
                        res_got = np.linalg.norm(B - np.dot(A, x))
                        # rtol = 10. as all the systems are products of orthonormals
                        # and on the small side (rows, cols) < 100.
                        np.testing.assert_allclose(
                            res_expected, res_got, rtol=10.)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(A, B, **kwargs)

        # test: column vector, tall, wide, square, row vector
        # prime sizes, the A's
        sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)]
        # compatible B's for Ax=B must have same number of rows and 1 or more
        # columns

        # This test takes ages! So combinations are trimmed via cycling

        # gets a dtype
        cycle_dt = cycle(self.dtypes)

        orders = ['F', 'C']
        # gets a memory order flag
        cycle_order = cycle(orders)

        # a specific condition number to use in the following tests
        # there is nothing special about it other than it is not magic
        specific_cond = 10.

        # inner test loop, extracted as there's additional logic etc required
        # that'd end up with this being repeated a lot
        def inner_test_loop_fn(A, dt, **kwargs):
            # test solve Ax=B for (column, matrix) B, same dtype as A
            b_sizes = (1, 13)

            for b_size in b_sizes:

                # check 2D B
                b_order = next(cycle_order)
                B = self.specific_sample_matrix(
                    (A.shape[0], b_size), dt, b_order)
                check(A, B, **kwargs)

                # check 1D B
                b_order = next(cycle_order)
                tmp = B[:, 0].copy(order=b_order)
                check(A, tmp, **kwargs)

        # test loop
        for a_size in sizes:

            dt = next(cycle_dt)
            a_order = next(cycle_order)

            # A full rank, well conditioned system
            A = self.specific_sample_matrix(a_size, dt, a_order)

            # run the test loop
            inner_test_loop_fn(A, dt)

            m, n = a_size
            minmn = min(m, n)

            # operations that only make sense with a 2D matrix system
            if m != 1 and n != 1:

                # Test a rank deficient system
                r = minmn - 1
                A = self.specific_sample_matrix(
                    a_size, dt, a_order, rank=r)
                # run the test loop
                inner_test_loop_fn(A, dt)

                # Test a system with a given condition number for use in
                # testing the rcond parameter.
                # This works because the singular values in the
                # specific_sample_matrix code are linspace (1, cond, [0... if
                # rank deficient])
                A = self.specific_sample_matrix(
                    a_size, dt, a_order, condition=specific_cond)
                # run the test loop
                rcond = 1. / specific_cond
                approx_half_rank_rcond = minmn * rcond
                inner_test_loop_fn(A, dt,
                                   rcond=approx_half_rank_rcond)

        # Test input validation
        ok = np.array([[1., 2.], [3., 4.]], dtype=np.float64)

        # check ok input is ok
        cfunc, (ok, ok)

        # check bad inputs
        rn = "lstsq"

        # Wrong dtype
        bad = np.array([[1, 2], [3, 4]], dtype=np.int32)
        self.assert_wrong_dtype(rn, cfunc, (ok, bad))
        self.assert_wrong_dtype(rn, cfunc, (bad, ok))

        # different dtypes
        bad = np.array([[1, 2], [3, 4]], dtype=np.float32)
        self.assert_homogeneous_dtypes(rn, cfunc, (ok, bad))
        self.assert_homogeneous_dtypes(rn, cfunc, (bad, ok))

        # Dimension issue
        bad = np.array([1, 2], dtype=np.float64)
        self.assert_wrong_dimensions(rn, cfunc, (bad, ok))

        # no nans or infs
        bad = np.array([[1., 2., ], [np.inf, np.nan]], dtype=np.float64)
        self.assert_no_nan_or_inf(cfunc, (ok, bad))
        self.assert_no_nan_or_inf(cfunc, (bad, ok))

        # check 1D is accepted for B (2D is done previously)
        # and then that anything of higher dimension raises
        oneD = np.array([1., 2.], dtype=np.float64)
        cfunc, (ok, oneD)
        bad = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]], dtype=np.float64)
        self.assert_wrong_dimensions_1D(rn, cfunc, (ok, bad))

        # check a dimensionally invalid system raises (1D and 2D cases
        # checked)
        bad1D = np.array([1.], dtype=np.float64)
        bad2D = np.array([[1.], [2.], [3.]], dtype=np.float64)
        self.assert_dimensionally_invalid(cfunc, (ok, bad1D))
        self.assert_dimensionally_invalid(cfunc, (ok, bad2D))


class TestLinalgSolve(TestLinalgSystems):
    """
    Tests for np.linalg.solve.
    """

    @needs_lapack
    def test_linalg_solve(self):
        """
        Test np.linalg.solve
        """
        cfunc = jit(nopython=True)(solve_system)

        def check(a, b, **kwargs):
            expected = solve_system(a, b, **kwargs)
            got = cfunc(a, b, **kwargs)

            # check that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "F")

            use_reconstruction = False
            # try plain match of the result first
            try:
                np.testing.assert_array_almost_equal_nulp(
                    got, expected, nulp=10)
            except AssertionError:
                # plain match failed, test by reconstruction
                use_reconstruction = True

            # If plain match fails then reconstruction is used,
            # this checks that AX ~= B.
            # Plain match can fail due to numerical fuzziness associated
            # with system size and conditioning, or more simply from
            # numpy using double precision routines for computation that
            # could be done in single precision (which is what numba does).
            # Therefore minor differences in results can appear due to
            # e.g. numerical roundoff being different between two precisions.
            if use_reconstruction:
                # check they are dimensionally correct
                self.assertEqual(got.shape, expected.shape)

                # check AX=B
                rec = np.dot(a, got)
                resolution = np.finfo(a.dtype).resolution
                np.testing.assert_allclose(
                    b,
                    rec,
                    rtol=10 * resolution,
                    atol=100 * resolution  # zeros tend to be fuzzy
                )

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, b, **kwargs)

        # test: prime size squares
        sizes = [(1, 1), (3, 3), (7, 7)]

        # test loop
        for size, dtype, order in \
                product(sizes, self.dtypes, 'FC'):
            A = self.specific_sample_matrix(size, dtype, order)

            b_sizes = (1, 13)

            for b_size, b_order in product(b_sizes, 'FC'):
                # check 2D B
                B = self.specific_sample_matrix(
                    (A.shape[0], b_size), dtype, b_order)
                check(A, B)

                # check 1D B
                tmp = B[:, 0].copy(order=b_order)
                check(A, tmp)

        # Test input validation
        ok = np.array([[1., 0.], [0., 1.]], dtype=np.float64)

        # check ok input is ok
        cfunc(ok, ok)

        # check bad inputs
        rn = "solve"

        # Wrong dtype
        bad = np.array([[1, 0], [0, 1]], dtype=np.int32)
        self.assert_wrong_dtype(rn, cfunc, (ok, bad))
        self.assert_wrong_dtype(rn, cfunc, (bad, ok))

        # different dtypes
        bad = np.array([[1, 2], [3, 4]], dtype=np.float32)
        self.assert_homogeneous_dtypes(rn, cfunc, (ok, bad))
        self.assert_homogeneous_dtypes(rn, cfunc, (bad, ok))

        # Dimension issue
        bad = np.array([1, 0], dtype=np.float64)
        self.assert_wrong_dimensions(rn, cfunc, (bad, ok))

        # no nans or infs
        bad = np.array([[1., 0., ], [np.inf, np.nan]], dtype=np.float64)
        self.assert_no_nan_or_inf(cfunc, (ok, bad))
        self.assert_no_nan_or_inf(cfunc, (bad, ok))

        # check 1D is accepted for B (2D is done previously)
        # and then that anything of higher dimension raises
        ok_oneD = np.array([1., 2.], dtype=np.float64)
        cfunc(ok, ok_oneD)
        bad = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]], dtype=np.float64)
        self.assert_wrong_dimensions_1D(rn, cfunc, (ok, bad))

        # check an invalid system raises (1D and 2D cases checked)
        bad1D = np.array([1.], dtype=np.float64)
        bad2D = np.array([[1.], [2.], [3.]], dtype=np.float64)
        self.assert_dimensionally_invalid(cfunc, (ok, bad1D))
        self.assert_dimensionally_invalid(cfunc, (ok, bad2D))

        # check that a singular system raises
        bad2D = self.specific_sample_matrix((2, 2), np.float64, 'C', rank=1)
        self.assert_raise_on_singular(cfunc, (bad2D, ok))


class TestLinalgPinv(TestLinalgBase):
    """
    Tests for np.linalg.pinv.
    """

    @needs_lapack
    def test_linalg_pinv(self):
        """
        Test np.linalg.pinv
        """
        cfunc = jit(nopython=True)(pinv_matrix)

        def check(a, **kwargs):
            expected = pinv_matrix(a, **kwargs)
            got = cfunc(a, **kwargs)

            # check that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "F")

            use_reconstruction = False
            # try plain match of each array to np first

            try:
                np.testing.assert_array_almost_equal_nulp(
                    got, expected, nulp=10)
            except AssertionError:
                # plain match failed, test by reconstruction
                use_reconstruction = True

            # If plain match fails then reconstruction is used.
            # This can occur due to numpy using double precision
            # LAPACK when single can be used, this creates round off
            # problems. Also, if the matrix has machine precision level
            # zeros in its singular values then the singular vectors are
            # likely to vary depending on round off.
            if use_reconstruction:

                # check they are dimensionally correct
                self.assertEqual(got.shape, expected.shape)

                # check pinv(A)*A~=eye
                # if the problem is numerical fuzz then this will probably
                # work, if the problem is rank deficiency then it won't!
                rec = np.dot(got, a)
                try:
                    self.assert_is_identity_matrix(rec)
                except AssertionError:
                    # check A=pinv(pinv(A))
                    resolution = 5 * np.finfo(a.dtype).resolution
                    rec = cfunc(got)
                    np.testing.assert_allclose(
                        rec,
                        a,
                        rtol=10 * resolution,
                        atol=100 * resolution  # zeros tend to be fuzzy
                    )
                    if a.shape[0] >= a.shape[1]:
                        # if it is overdetermined or fully determined
                        # use numba lstsq function (which is type specific) to
                        # compute the inverse and check against that.
                        lstsq = jit(nopython=True)(lstsq_system)
                        lstsq_pinv = lstsq(
                            a, np.eye(
                                a.shape[0]).astype(
                                a.dtype), **kwargs)[0]
                        np.testing.assert_allclose(
                            got,
                            lstsq_pinv,
                            rtol=10 * resolution,
                            atol=100 * resolution  # zeros tend to be fuzzy
                        )
                    # check the 2 norm of the difference is small
                    self.assertLess(np.linalg.norm(got - expected), resolution)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, **kwargs)

        # test: column vector, tall, wide, square, row vector
        # prime sizes
        sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)]

        # When required, a specified condition number
        specific_cond = 10.

        # test loop
        for size, dtype, order in \
                product(sizes, self.dtypes, 'FC'):
            # check a full rank matrix
            a = self.specific_sample_matrix(size, dtype, order)
            check(a)

            m, n = size
            if m != 1 and n != 1:
                # check a rank deficient matrix
                minmn = min(m, n)
                a = self.specific_sample_matrix(size, dtype, order,
                                                condition=specific_cond)
                rcond = 1. / specific_cond
                approx_half_rank_rcond = minmn * rcond
                check(a, rcond=approx_half_rank_rcond)

        rn = "pinv"

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions(rn, cfunc,
                                     (np.ones(10, dtype=np.float64),))

        # no nans or infs
        self.assert_no_nan_or_inf(cfunc,
                                  (np.array([[1., 2., ], [np.inf, np.nan]],
                                            dtype=np.float64),))


class TestLinalgDetAndSlogdet(TestLinalgBase):
    """
    Tests for np.linalg.det. and np.linalg.slogdet.
    Exactly the same inputs are used for both tests as
    det() is a trivial function of slogdet(), the tests
    are therefore combined.
    """

    def check_det(self, cfunc, a, **kwargs):
        expected = det_matrix(a, **kwargs)
        got = cfunc(a, **kwargs)

        resolution = 5 * np.finfo(a.dtype).resolution

        # check the determinants are the same
        np.testing.assert_allclose(got, expected, rtol=resolution)

        # Ensure proper resource management
        with self.assertNoNRTLeak():
            cfunc(a, **kwargs)

    def check_slogdet(self, cfunc, a, **kwargs):
        expected = slogdet_matrix(a, **kwargs)
        got = cfunc(a, **kwargs)

        # As numba returns python floats types and numpy returns
        # numpy float types, some more adjustment and different
        # types of comparison than those used with array based
        # results is required.

        # check that the returned tuple is same length
        self.assertEqual(len(expected), len(got))
        # and that length is 2
        self.assertEqual(len(got), 2)

        # check that the domain of the results match
        for k in range(2):
            self.assertEqual(
                np.iscomplexobj(got[k]),
                np.iscomplexobj(expected[k]))

        # turn got[0] into the same dtype as `a`
        # this is so checking with nulp will work
        got_conv = a.dtype.type(got[0])
        np.testing.assert_array_almost_equal_nulp(
            got_conv, expected[0], nulp=10)
        # compare log determinant magnitude with a more fuzzy value
        # as numpy values come from higher precision lapack routines
        resolution = 5 * np.finfo(a.dtype).resolution
        np.testing.assert_allclose(
            got[1], expected[1], rtol=resolution, atol=resolution)

        # Ensure proper resource management
        with self.assertNoNRTLeak():
            cfunc(a, **kwargs)

    def do_test(self, rn, check, cfunc):

        # test: 1x1 as it is unusual, 4x4 as it is even and 7x7 as is it odd!
        sizes = [(1, 1), (4, 4), (7, 7)]

        # test loop
        for size, dtype, order in \
                product(sizes, self.dtypes, 'FC'):
            # check a full rank matrix
            a = self.specific_sample_matrix(size, dtype, order)
            check(cfunc, a)

        # use a matrix of zeros to trip xgetrf U(i,i)=0 singular test
        for dtype, order in product(self.dtypes, 'FC'):
            a = np.zeros((3, 3), dtype=dtype)
            check(cfunc, a)

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions(rn, cfunc,
                                     (np.ones(10, dtype=np.float64),))

        # no nans or infs
        self.assert_no_nan_or_inf(cfunc,
                                  (np.array([[1., 2., ], [np.inf, np.nan]],
                                            dtype=np.float64),))

    @needs_lapack
    def test_linalg_det(self):
        cfunc = jit(nopython=True)(det_matrix)
        self.do_test("det", self.check_det, cfunc)

    @needs_lapack
    def test_linalg_slogdet(self):
        cfunc = jit(nopython=True)(slogdet_matrix)
        self.do_test("slogdet", self.check_slogdet, cfunc)

# Use TestLinalgSystems as a base to get access to additional
# testing for 1 and 2D inputs.


class TestLinalgNorm(TestLinalgSystems):
    """
    Tests for np.linalg.norm.
    """

    @needs_lapack
    def test_linalg_norm(self):
        """
        Test np.linalg.norm
        """
        cfunc = jit(nopython=True)(norm_matrix)

        def check(a, **kwargs):
            expected = norm_matrix(a, **kwargs)
            got = cfunc(a, **kwargs)

            # All results should be in the real domain
            self.assertTrue(not np.iscomplexobj(got))

            resolution = 5 * np.finfo(a.dtype).resolution

            # check the norms are the same to the arg `a` precision
            np.testing.assert_allclose(got, expected, rtol=resolution)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, **kwargs)

        # Check 1D inputs
        sizes = [1, 4, 7]
        nrm_types = [None, np.inf, -np.inf, 0, 1, -1, 2, -2, 5, 6.7, -4.3]

        # standard 1D input
        for size, dtype, nrm_type in \
                product(sizes, self.dtypes, nrm_types):
            a = self.sample_vector(size, dtype)
            check(a, ord=nrm_type)

        # sliced 1D input
        for dtype, nrm_type in \
                product(self.dtypes, nrm_types):
            a = self.sample_vector(10, dtype)[::3]
            check(a, ord=nrm_type)

        # check that numba returns zero for empty arrays. Numpy returns zero
        # for most norm types and raises ValueError for +/-np.inf.
        for dtype, nrm_type, order in \
                product(self.dtypes, nrm_types, 'FC'):
            a = np.array([], dtype=dtype, order=order)
            self.assertEqual(cfunc(a, nrm_type), 0.0)

        # Check 2D inputs:
        # test: column vector, tall, wide, square, row vector
        # prime sizes
        sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)]
        nrm_types = [None, np.inf, -np.inf, 1, -1, 2, -2]

        # standard 2D input
        for size, dtype, order, nrm_type in \
                product(sizes, self.dtypes, 'FC', nrm_types):
            # check a full rank matrix
            a = self.specific_sample_matrix(size, dtype, order)
            check(a, ord=nrm_type)

        # check 2D slices work for the case where xnrm2 is called from
        # BLAS (ord=None) to make sure it is working ok.
        nrm_types = [None]
        for dtype, nrm_type, order in \
                product(self.dtypes, nrm_types, 'FC'):
            a = self.specific_sample_matrix((17, 13), dtype, order)
            # contig for C order
            check(a[:3], ord=nrm_type)

            # contig for Fortran order
            check(a[:, 3:], ord=nrm_type)

            # contig for neither order
            check(a[1, 4::3], ord=nrm_type)

        # check that numba returns zero for empty arrays. Numpy returns zero
        # for some norm types and raises a variety of errors for others.
        for dtype, nrm_type, order in \
                product(self.dtypes, nrm_types, 'FC'):
            a = np.array([[]], dtype=dtype, order=order)
            self.assertEqual(cfunc(a, nrm_type), 0.0)

        rn = "norm"

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue, reuse the test from the TestLinalgSystems class
        self.assert_wrong_dimensions_1D(
            rn, cfunc, (np.ones(
                12, dtype=np.float64).reshape(
                2, 2, 3),))

        # no nans or infs for 2d case when SVD is used (e.g 2-norm)
        self.assert_no_nan_or_inf(cfunc,
                                  (np.array([[1., 2.], [np.inf, np.nan]],
                                            dtype=np.float64), 2))

        # assert 2D input raises for an invalid norm kind kwarg
        self.assert_invalid_norm_kind(cfunc, (np.array([[1., 2.], [3., 4.]],
                                                       dtype=np.float64), 6))


class TestLinalgCond(TestLinalgBase):
    """
    Tests for np.linalg.cond.
    """

    @needs_lapack
    def test_linalg_cond(self):
        """
        Test np.linalg.cond
        """

        cfunc = jit(nopython=True)(cond_matrix)

        def check(a, **kwargs):
            expected = cond_matrix(a, **kwargs)
            got = cfunc(a, **kwargs)

            # All results should be in the real domain
            self.assertTrue(not np.iscomplexobj(got))

            resolution = 5 * np.finfo(a.dtype).resolution

            # check the cond is the same to the arg `a` precision
            np.testing.assert_allclose(got, expected, rtol=resolution)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, **kwargs)

        # valid p values (used to indicate norm type)
        ps = [None, np.inf, -np.inf, 1, -1, 2, -2]
        sizes = [(3, 3), (7, 7)]

        for size, dtype, order, p in \
                product(sizes, self.dtypes, 'FC', ps):
            a = self.specific_sample_matrix(size, dtype, order)
            check(a, p=p)

        # When p=None non-square matrices are accepted.
        sizes = [(7, 1), (11, 5), (5, 11), (1, 7)]
        for size, dtype, order in \
                product(sizes, self.dtypes, 'FC'):
            a = self.specific_sample_matrix(size, dtype, order)
            check(a)

        # try an ill-conditioned system with 2-norm, make sure np raises an
        # overflow warning as the result is `+inf` and that the result from
        # numba matches.
        with warnings.catch_warnings():
            a = np.array([[1.e308, 0], [0, 0.1]], dtype=np.float64)
            warnings.simplefilter("error", RuntimeWarning)
            self.assertRaisesRegexp(RuntimeWarning,
                                    'overflow encountered in.*',
                                    check, a)
            warnings.simplefilter("ignore", RuntimeWarning)
            check(a)

        rn = "cond"

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions(rn, cfunc,
                                     (np.ones(10, dtype=np.float64),))

        # no nans or infs when p="None" (default for kwarg).
        self.assert_no_nan_or_inf(cfunc,
                                  (np.array([[1., 2., ], [np.inf, np.nan]],
                                            dtype=np.float64),))

        # assert raises for an invalid norm kind kwarg
        self.assert_invalid_norm_kind(cfunc, (np.array([[1., 2.], [3., 4.]],
                                                       dtype=np.float64), 6))


class TestLinalgMatrixRank(TestLinalgSystems):
    """
    Tests for np.linalg.matrix_rank.
    """

    @needs_lapack
    def test_linalg_matrix_rank(self):
        """
        Test np.linalg.matrix_rank
        """

        cfunc = jit(nopython=True)(matrix_rank_matrix)

        def check(a, **kwargs):
            expected = matrix_rank_matrix(a, **kwargs)
            got = cfunc(a, **kwargs)

            # Ranks are integral so comparison should be trivial.
            # check the rank is the same
            np.testing.assert_allclose(got, expected)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, **kwargs)

        sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)]

        for size, dtype, order in \
                product(sizes, self.dtypes, 'FC'):
            # check full rank system
            a = self.specific_sample_matrix(size, dtype, order)
            check(a)

            # If the system is a matrix, check rank deficiency is reported
            # correctly. Check all ranks from 0 to (full rank - 1).
            tol = 1e-13
            # first check 1 to (full rank - 1)
            for k in range(1, min(size) - 1):
                # check rank k
                a = self.specific_sample_matrix(size, dtype, order, rank=k)
                self.assertEqual(cfunc(a), k)
                check(a)
                # check provision of a tolerance works as expected
                # create a (m x n) diagonal matrix with a singular value
                # guaranteed below the tolerance 1e-13
                m, n = a.shape
                a[:, :] = 0.  # reuse `a`'s memory
                idx = np.nonzero(np.eye(m, n))
                if np.iscomplexobj(a):
                    b = 1. + np.random.rand(k) + 1.j +\
                        1.j * np.random.rand(k)
                    # min singular value is sqrt(2)*1e-14
                    b[0] = 1e-14 + 1e-14j
                else:
                    b = 1. + np.random.rand(k)
                    b[0] = 1e-14  # min singular value is 1e-14
                a[idx[0][:k], idx[1][:k]] = b.astype(dtype)
                # rank should be k-1 (as tol is present)
                self.assertEqual(cfunc(a, tol), k - 1)
                check(a, tol=tol)
            # then check zero rank
            a[:, :] = 0.
            self.assertEqual(cfunc(a), 0)
            check(a)
            # add in a singular value that is small
            if np.iscomplexobj(a):
                a[-1, -1] = 1e-14 + 1e-14j
            else:
                a[-1, -1] = 1e-14
            # check the system has zero rank to a given tolerance
            self.assertEqual(cfunc(a, tol), 0)
            check(a, tol=tol)

        # check the zero vector returns rank 0 and a nonzero vector
        # returns rank 1.
        for dt in self.dtypes:
            a = np.zeros((5), dtype=dt)
            self.assertEqual(cfunc(a), 0)
            check(a)
            # make it a nonzero vector
            a[0] = 1.
            self.assertEqual(cfunc(a), 1)
            check(a)

        rn = "matrix_rank"

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32),))

        # Dimension issue
        self.assert_wrong_dimensions_1D(
            rn, cfunc, (np.ones(
                12, dtype=np.float64).reshape(
                2, 2, 3),))

        # no nans or infs for 2D case
        self.assert_no_nan_or_inf(cfunc,
                                  (np.array([[1., 2., ], [np.inf, np.nan]],
                                            dtype=np.float64),))


class TestLinalgMatrixPower(TestLinalgBase):
    """
    Tests for np.linalg.matrix_power.
    """

    def assert_int_exponenent(self, cfunc, args):
        # validate first arg is ok
        cfunc(args[0], 1)
        # pass in both args and assert fail
        with self.assertRaises(errors.TypingError):
            cfunc(*args)

    @needs_lapack
    def test_linalg_matrix_power(self):
        cfunc = jit(nopython=True)(matrix_power_matrix)

        def check(a, pwr):
            expected = matrix_power_matrix(a, pwr)
            got = cfunc(a, pwr)

            # check that the computed results are contig and in the same way
            self.assert_contig_sanity(got, "C")

            res = 5 * np.finfo(a.dtype).resolution
            np.testing.assert_allclose(got, expected, rtol=res, atol=res)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, pwr)

        sizes = [(1, 1), (5, 5), (7, 7)]
        powers = [-33, -17] + list(range(-10, 10)) + [17, 33]

        for size, pwr, dtype, order in \
                product(sizes, powers, self.dtypes, 'FC'):
            a = self.specific_sample_matrix(size, dtype, order)
            check(a, pwr)

        rn = "matrix_power"

        # Wrong dtype
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32), 1))

        # not an int power
        self.assert_wrong_dtype(rn, cfunc,
                                (np.ones((2, 2), dtype=np.int32), 1))

        # Dimension issue
        self.assert_wrong_dimensions(rn, cfunc,
                                     (np.ones(10, dtype=np.float64), 1))

        # non-integer supplied as exponent
        self.assert_int_exponenent(cfunc, (np.ones((2, 2)), 1.2))

        # singular matrix is not invertible
        self.assert_raise_on_singular(cfunc, (np.array([[0., 0], [1, 1]]), -1))


class TestTrace(TestLinalgBase):
    """
    Tests for np.trace.
    """

    def setUp(self):
        super(TestTrace, self).setUp()
        # compile two versions, one with and one without the offset kwarg
        self.cfunc_w_offset = jit(nopython=True)(trace_matrix)
        self.cfunc_no_offset = jit(nopython=True)(trace_matrix_no_offset)

    def assert_int_offset(self, cfunc, a, **kwargs):
        # validate first arg is ok
        cfunc(a)
        # pass in kwarg and assert fail
        with self.assertRaises(errors.TypingError):
            cfunc(a, **kwargs)

    def test_trace(self):

        def check(a, **kwargs):
            if 'offset' in kwargs:
                expected = trace_matrix(a, **kwargs)
                cfunc = self.cfunc_w_offset
            else:
                expected = trace_matrix_no_offset(a, **kwargs)
                cfunc = self.cfunc_no_offset

            got = cfunc(a, **kwargs)

            res = 5 * np.finfo(a.dtype).resolution
            np.testing.assert_allclose(got, expected, rtol=res, atol=res)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, **kwargs)

        # test: column vector, tall, wide, square, row vector
        # prime sizes
        sizes = [(7, 1), (11, 5), (5, 11), (3, 3), (1, 7)]

        # offsets to cover the range of the matrix sizes above
        offsets = [-13, -12, -11] + list(range(-10, 10)) + [11, 12, 13]

        for size, offset, dtype, order in \
                product(sizes, offsets, self.dtypes, 'FC'):
            a = self.specific_sample_matrix(size, dtype, order)
            check(a, offset=offset)
            if offset == 0:
                check(a)

        rn = "trace"

        # Dimension issue
        self.assert_wrong_dimensions(rn, self.cfunc_w_offset,
                                     (np.ones(10, dtype=np.float64), 1), False)
        self.assert_wrong_dimensions(rn, self.cfunc_no_offset,
                                     (np.ones(10, dtype=np.float64),), False)

        # non-integer supplied as exponent
        self.assert_int_offset(
            self.cfunc_w_offset, np.ones(
                (2, 2)), offset=1.2)

    def test_trace_w_optional_input(self):
        "Issue 2314"
        @jit("(optional(float64[:,:]),)", nopython=True)
        def tested(a):
            return np.trace(a)

        a = np.ones((5, 5), dtype=np.float64)
        tested(a)

        with self.assertRaises(TypeError) as raises:
            tested(None)

        errmsg = str(raises.exception)
        self.assertEqual('expected array(float64, 2d, A), got None', errmsg)


class TestBasics(TestLinalgSystems):  # TestLinalgSystems for 1d test

    order1 = cycle(['F', 'C', 'C', 'F'])
    order2 = cycle(['C', 'F', 'C', 'F'])

    # test: column vector, matrix, row vector, 1d sizes
    # (7, 1, 3) and two scalars
    sizes = [(7, 1), (3, 3), (1, 7), (7,), (1,), (3,), 3., 5.]

    def _assert_wrong_dim(self, rn, cfunc):
        # Dimension issue
        self.assert_wrong_dimensions_1D(
            rn, cfunc, (np.array([[[1]]], dtype=np.float64), np.ones(1)), False)
        self.assert_wrong_dimensions_1D(
            rn, cfunc, (np.ones(1), np.array([[[1]]], dtype=np.float64)), False)

    def _gen_input(self, size, dtype, order):
        if not isinstance(size, tuple):
            return size
        else:
            if len(size) == 1:
                return self.sample_vector(size[0], dtype)
            else:
                return self.sample_vector(
                    size[0] * size[1],
                    dtype).reshape(
                    size, order=order)

    def _get_input(self, size1, size2, dtype):
        a = self._gen_input(size1, dtype, next(self.order1))
        b = self._gen_input(size2, dtype, next(self.order2))
        # force domain consistency as underlying ufuncs require it
        if np.iscomplexobj(a):
            b = b + 1j
        if np.iscomplexobj(b):
            a = a + 1j
        return (a, b)

    def test_outer(self):
        cfunc = jit(nopython=True)(outer_matrix)

        def check(a, b, **kwargs):

            # check without kwargs
            expected = outer_matrix(a, b)
            got = cfunc(a, b)

            res = 5 * np.finfo(np.asarray(a).dtype).resolution
            np.testing.assert_allclose(got, expected, rtol=res, atol=res)

            # if kwargs present check with them too
            if 'out' in kwargs:
                got = cfunc(a, b, **kwargs)
                np.testing.assert_allclose(got, expected, rtol=res,
                                           atol=res)
                np.testing.assert_allclose(kwargs['out'], expected,
                                           rtol=res, atol=res)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, b, **kwargs)

        for size1, size2, dtype in \
                product(self.sizes, self.sizes, self.dtypes):
            (a, b) = self._get_input(size1, size2, dtype)
            check(a, b)
            if numpy_version >= (1, 9):
                c = np.empty((np.asarray(a).size, np.asarray(b).size),
                             dtype=np.asarray(a).dtype)
                check(a, b, out=c)

        self._assert_wrong_dim("outer", cfunc)

    def test_kron(self):
        cfunc = jit(nopython=True)(kron_matrix)

        def check(a, b, **kwargs):

            expected = kron_matrix(a, b)
            got = cfunc(a, b)

            res = 5 * np.finfo(np.asarray(a).dtype).resolution
            np.testing.assert_allclose(got, expected, rtol=res, atol=res)

            # Ensure proper resource management
            with self.assertNoNRTLeak():
                cfunc(a, b)

        for size1, size2, dtype in \
                product(self.sizes, self.sizes, self.dtypes):
            (a, b) = self._get_input(size1, size2, dtype)
            check(a, b)

        self._assert_wrong_dim("kron", cfunc)


if __name__ == '__main__':
    unittest.main()