# # Created by: Pearu Peterson, April 2002 # import math import pytest import numpy as np from numpy.testing import (assert_equal, assert_almost_equal, assert_array_almost_equal, assert_allclose) from pytest import raises as assert_raises from numpy import (arange, triu, tril, zeros, tril_indices, ones, diag, append, eye, nonzero) import scipy from scipy.linalg import _fblas as fblas, get_blas_funcs, toeplitz, solve try: from scipy.linalg import _cblas as cblas except ImportError: cblas = None REAL_DTYPES = [np.float32, np.float64] COMPLEX_DTYPES = [np.complex64, np.complex128] DTYPES = REAL_DTYPES + COMPLEX_DTYPES def test_get_blas_funcs(): # check that it returns Fortran code for arrays that are # fortran-ordered f1, f2, f3 = get_blas_funcs( ('axpy', 'axpy', 'axpy'), (np.empty((2, 2), dtype=np.complex64, order='F'), np.empty((2, 2), dtype=np.complex128, order='C')) ) # get_blas_funcs will choose libraries depending on most generic # array assert_equal(f1.typecode, 'z') assert_equal(f2.typecode, 'z') if cblas is not None: assert_equal(f1.module_name, 'cblas') assert_equal(f2.module_name, 'cblas') # check defaults. f1 = get_blas_funcs('rotg') assert_equal(f1.typecode, 'd') # check also dtype interface f1 = get_blas_funcs('gemm', dtype=np.complex64) assert_equal(f1.typecode, 'c') f1 = get_blas_funcs('gemm', dtype='F') assert_equal(f1.typecode, 'c') # extended precision complex f1 = get_blas_funcs('gemm', dtype=np.clongdouble) assert_equal(f1.typecode, 'z') # check safe complex upcasting f1 = get_blas_funcs('axpy', (np.empty((2, 2), dtype=np.float64), np.empty((2, 2), dtype=np.complex64)) ) assert_equal(f1.typecode, 'z') def test_get_blas_funcs_alias(): # check alias for get_blas_funcs f, g = get_blas_funcs(('nrm2', 'dot'), dtype=np.complex64) assert f.typecode == 'c' assert g.typecode == 'c' f, g, h = get_blas_funcs(('dot', 'dotc', 'dotu'), dtype=np.float64) assert f is g assert f is h def parametrize_blas(mod, func_name, prefixes): if mod is None: return pytest.mark.skip(reason="cblas not available") params = [] for prefix in prefixes: if 'z' in prefix: dtype = np.complex128 elif 'c' in prefix: dtype = np.complex64 elif 'd' in prefix: dtype = np.float64 else: assert 's' in prefix dtype = np.float32 f = getattr(mod, prefix + func_name) params.append(pytest.param(f, dtype, id=prefix + func_name)) return pytest.mark.parametrize("f,dtype", params) class TestCBLAS1Simple: @parametrize_blas(cblas, "axpy", "sdcz") def test_axpy(self, f, dtype): assert_array_almost_equal(f([1, 2, 3], [2, -1, 3], a=5), [7, 9, 18]) if dtype in COMPLEX_DTYPES: assert_array_almost_equal(f([1, 2j, 3], [2, -1, 3], a=5), [7, 10j-1, 18]) class TestFBLAS1Simple: @parametrize_blas(fblas, "axpy", "sdcz") def test_axpy(self, f, dtype): assert_array_almost_equal(f([1, 2, 3], [2, -1, 3], a=5), [7, 9, 18]) if dtype in COMPLEX_DTYPES: assert_array_almost_equal(f([1, 2j, 3], [2, -1, 3], a=5), [7, 10j-1, 18]) @parametrize_blas(fblas, "copy", "sdcz") def test_copy(self, f, dtype): assert_array_almost_equal(f([3, 4, 5], [8]*3), [3, 4, 5]) if dtype in COMPLEX_DTYPES: assert_array_almost_equal(f([3, 4j, 5+3j], [8]*3), [3, 4j, 5+3j]) @parametrize_blas(fblas, "asum", ["s", "d", "sc", "dz"]) def test_asum(self, f, dtype): assert_almost_equal(f([3, -4, 5]), 12) if dtype in COMPLEX_DTYPES: assert_almost_equal(f([3j, -4, 3-4j]), 14) @parametrize_blas(fblas, "dot", "sd") def test_dot(self, f, dtype): assert_almost_equal(f([3, -4, 5], [2, 5, 1]), -9) @parametrize_blas(fblas, "dotu", "cz") def test_dotu(self, f, dtype): assert_almost_equal(f([3j, -4, 3-4j], [2, 3, 1]), -9+2j) @parametrize_blas(fblas, "dotc", "cz") def test_dotc(self, f, dtype): assert_almost_equal(f([3j, -4, 3-4j], [2, 3j, 1]), 3-14j) @parametrize_blas(fblas, "nrm2", ["s", "d", "sc", "dz"]) def test_nrm2(self, f, dtype): assert_almost_equal(f([3, -4, 5]), math.sqrt(50)) if dtype in COMPLEX_DTYPES: assert_almost_equal(f([3j, -4, 3-4j]), math.sqrt(50)) @parametrize_blas(fblas, "scal", ["s", "d", "cs", "zd"]) def test_scal(self, f, dtype): assert_array_almost_equal(f(2, [3, -4, 5]), [6, -8, 10]) if dtype in COMPLEX_DTYPES: assert_array_almost_equal(f(3, [3j, -4, 3-4j]), [9j, -12, 9-12j]) @parametrize_blas(fblas, "swap", "sdcz") def test_swap(self, f, dtype): x, y = [2, 3, 1], [-2, 3, 7] x1, y1 = f(x, y) assert_array_almost_equal(x1, y) assert_array_almost_equal(y1, x) if dtype in COMPLEX_DTYPES: x, y = [2, 3j, 1], [-2, 3, 7-3j] x1, y1 = f(x, y) assert_array_almost_equal(x1, y) assert_array_almost_equal(y1, x) @parametrize_blas(fblas, "amax", ["is", "id", "ic", "iz"]) def test_amax(self, f, dtype): assert_equal(f([-2, 4, 3]), 1) if dtype in COMPLEX_DTYPES: assert_equal(f([-5, 4+3j, 6]), 1) # XXX: need tests for rot,rotm,rotg,rotmg class TestFBLAS2Simple: @parametrize_blas(fblas, "gemv", "sdcz") def test_gemv(self, f, dtype): assert_array_almost_equal(f(3, [[3]], [-4]), [-36]) assert_array_almost_equal(f(3, [[3]], [-4], 3, [5]), [-21]) if dtype in COMPLEX_DTYPES: assert_array_almost_equal(f(3j, [[3-4j]], [-4]), [-48-36j]) assert_array_almost_equal(f(3j, [[3-4j]], [-4], 3, [5j]), [-48-21j]) @parametrize_blas(fblas, "ger", "sd") def test_ger(self, f, dtype): assert_array_almost_equal(f(1, [1, 2], [3, 4]), [[3, 4], [6, 8]]) assert_array_almost_equal(f(2, [1, 2, 3], [3, 4]), [[6, 8], [12, 16], [18, 24]]) assert_array_almost_equal(f(1, [1, 2], [3, 4], a=[[1, 2], [3, 4]]), [[4, 6], [9, 12]]) if dtype in COMPLEX_DTYPES: assert_array_almost_equal(f(1, [1j, 2], [3, 4]), [[3j, 4j], [6, 8]]) assert_array_almost_equal(f(2, [1j, 2j, 3j], [3j, 4j]), [[6, 8], [12, 16], [18, 24]]) @parametrize_blas(fblas, "geru", "cz") def test_geru(self, f, dtype): assert_array_almost_equal(f(1, [1j, 2], [3, 4]), [[3j, 4j], [6, 8]]) assert_array_almost_equal(f(-2, [1j, 2j, 3j], [3j, 4j]), [[6, 8], [12, 16], [18, 24]]) @parametrize_blas(fblas, "gerc", "cz") def test_gerc(self, f, dtype): assert_array_almost_equal(f(1, [1j, 2], [3, 4]), [[3j, 4j], [6, 8]]) assert_array_almost_equal(f(2, [1j, 2j, 3j], [3j, 4j]), [[6, 8], [12, 16], [18, 24]]) @parametrize_blas(fblas, "syr", "sdcz") def test_syr(self, f, dtype): x = np.arange(1, 5, dtype='d') resx = np.triu(x[:, np.newaxis] * x) resx_reverse = np.triu(x[::-1, np.newaxis] * x[::-1]) y = np.linspace(0, 8.5, 17, endpoint=False) z = np.arange(1, 9, dtype='d').view('D') resz = np.triu(z[:, np.newaxis] * z) resz_reverse = np.triu(z[::-1, np.newaxis] * z[::-1]) w = np.c_[np.zeros(4), z, np.zeros(4)].ravel() rtol = np.finfo(dtype).eps assert_allclose(f(1.0, x), resx, rtol=rtol) assert_allclose(f(1.0, x, lower=True), resx.T, rtol=rtol) assert_allclose(f(1.0, y, incx=2, offx=2, n=4), resx, rtol=rtol) # negative increments imply reversed vectors in blas assert_allclose(f(1.0, y, incx=-2, offx=2, n=4), resx_reverse, rtol=rtol) if dtype in COMPLEX_DTYPES: assert_allclose(f(1.0, z), resz, rtol=rtol) assert_allclose(f(1.0, z, lower=True), resz.T, rtol=rtol) assert_allclose(f(1.0, w, incx=3, offx=1, n=4), resz, rtol=rtol) # negative increments imply reversed vectors in blas assert_allclose(f(1.0, w, incx=-3, offx=1, n=4), resz_reverse, rtol=rtol) a = np.zeros((4, 4), dtype, 'F') b = f(1.0, z, a=a, overwrite_a=True) assert_allclose(a, resz, rtol=rtol) b = f(2.0, z, a=a) assert a is not b assert_allclose(b, 3*resz, rtol=rtol) else: a = np.zeros((4, 4), dtype, 'F') b = f(1.0, x, a=a, overwrite_a=True) assert_allclose(a, resx, rtol=rtol) b = f(2.0, x, a=a) assert a is not b assert_allclose(b, 3*resx, rtol=rtol) assert_raises(Exception, f, 1.0, x, incx=0) assert_raises(Exception, f, 1.0, x, offx=5) assert_raises(Exception, f, 1.0, x, offx=-2) assert_raises(Exception, f, 1.0, x, n=-2) assert_raises(Exception, f, 1.0, x, n=5) assert_raises(Exception, f, 1.0, x, lower=2) assert_raises(Exception, f, 1.0, x, a=np.zeros((2, 2), 'd', 'F')) @parametrize_blas(fblas, "her", "cz") def test_her(self, f, dtype): x = np.arange(1, 5, dtype='d') z = np.arange(1, 9, dtype='d').view('D') rehz = np.triu(z[:, np.newaxis] * z.conj()) rehz_reverse = np.triu(z[::-1, np.newaxis] * z[::-1].conj()) w = np.c_[np.zeros(4), z, np.zeros(4)].ravel() rtol = np.finfo(dtype).eps assert_allclose(f(1.0, z), rehz, rtol=rtol) assert_allclose(f(1.0, z, lower=True), rehz.T.conj(), rtol=rtol) assert_allclose(f(1.0, w, incx=3, offx=1, n=4), rehz, rtol=rtol) # negative increments imply reversed vectors in blas assert_allclose(f(1.0, w, incx=-3, offx=1, n=4), rehz_reverse, rtol=rtol) a = np.zeros((4, 4), dtype, 'F') b = f(1.0, z, a=a, overwrite_a=True) assert_allclose(a, rehz, rtol=rtol) b = f(2.0, z, a=a) assert a is not b assert_allclose(b, 3*rehz, rtol=rtol) assert_raises(Exception, f, 1.0, x, incx=0) assert_raises(Exception, f, 1.0, x, offx=5) assert_raises(Exception, f, 1.0, x, offx=-2) assert_raises(Exception, f, 1.0, x, n=-2) assert_raises(Exception, f, 1.0, x, n=5) assert_raises(Exception, f, 1.0, x, lower=2) assert_raises(Exception, f, 1.0, x, a=np.zeros((2, 2), 'd', 'F')) @parametrize_blas(fblas, "syr2", "sd") def test_syr2(self, f, dtype): x = np.arange(1, 5, dtype='d') y = np.arange(5, 9, dtype='d') resxy = np.triu(x[:, np.newaxis] * y + y[:, np.newaxis] * x) resxy_reverse = np.triu(x[::-1, np.newaxis] * y[::-1] + y[::-1, np.newaxis] * x[::-1]) q = np.linspace(0, 8.5, 17, endpoint=False) rtol = np.finfo(dtype).eps assert_allclose(f(1.0, x, y), resxy, rtol=rtol) assert_allclose(f(1.0, x, y, n=3), resxy[:3, :3], rtol=rtol) assert_allclose(f(1.0, x, y, lower=True), resxy.T, rtol=rtol) assert_allclose(f(1.0, q, q, incx=2, offx=2, incy=2, offy=10), resxy, rtol=rtol) assert_allclose(f(1.0, q, q, incx=2, offx=2, incy=2, offy=10, n=3), resxy[:3, :3], rtol=rtol) # negative increments imply reversed vectors in blas assert_allclose(f(1.0, q, q, incx=-2, offx=2, incy=-2, offy=10), resxy_reverse, rtol=rtol) a = np.zeros((4, 4), dtype, 'F') b = f(1.0, x, y, a=a, overwrite_a=True) assert_allclose(a, resxy, rtol=rtol) b = f(2.0, x, y, a=a) assert a is not b assert_allclose(b, 3*resxy, rtol=rtol) assert_raises(Exception, f, 1.0, x, y, incx=0) assert_raises(Exception, f, 1.0, x, y, offx=5) assert_raises(Exception, f, 1.0, x, y, offx=-2) assert_raises(Exception, f, 1.0, x, y, incy=0) assert_raises(Exception, f, 1.0, x, y, offy=5) assert_raises(Exception, f, 1.0, x, y, offy=-2) assert_raises(Exception, f, 1.0, x, y, n=-2) assert_raises(Exception, f, 1.0, x, y, n=5) assert_raises(Exception, f, 1.0, x, y, lower=2) assert_raises(Exception, f, 1.0, x, y, a=np.zeros((2, 2), 'd', 'F')) @parametrize_blas(fblas, "her2", "cz") def test_her2(self, f, dtype): x = np.arange(1, 9, dtype='d').view('D') y = np.arange(9, 17, dtype='d').view('D') resxy = x[:, np.newaxis] * y.conj() + y[:, np.newaxis] * x.conj() resxy = np.triu(resxy) resxy_reverse = x[::-1, np.newaxis] * y[::-1].conj() resxy_reverse += y[::-1, np.newaxis] * x[::-1].conj() resxy_reverse = np.triu(resxy_reverse) u = np.c_[np.zeros(4), x, np.zeros(4)].ravel() v = np.c_[np.zeros(4), y, np.zeros(4)].ravel() rtol = np.finfo(dtype).eps assert_allclose(f(1.0, x, y), resxy, rtol=rtol) assert_allclose(f(1.0, x, y, n=3), resxy[:3, :3], rtol=rtol) assert_allclose(f(1.0, x, y, lower=True), resxy.T.conj(), rtol=rtol) assert_allclose(f(1.0, u, v, incx=3, offx=1, incy=3, offy=1), resxy, rtol=rtol) assert_allclose(f(1.0, u, v, incx=3, offx=1, incy=3, offy=1, n=3), resxy[:3, :3], rtol=rtol) # negative increments imply reversed vectors in blas assert_allclose(f(1.0, u, v, incx=-3, offx=1, incy=-3, offy=1), resxy_reverse, rtol=rtol) a = np.zeros((4, 4), dtype, 'F') b = f(1.0, x, y, a=a, overwrite_a=True) assert_allclose(a, resxy, rtol=rtol) b = f(2.0, x, y, a=a) assert a is not b assert_allclose(b, 3*resxy, rtol=rtol) assert_raises(Exception, f, 1.0, x, y, incx=0) assert_raises(Exception, f, 1.0, x, y, offx=5) assert_raises(Exception, f, 1.0, x, y, offx=-2) assert_raises(Exception, f, 1.0, x, y, incy=0) assert_raises(Exception, f, 1.0, x, y, offy=5) assert_raises(Exception, f, 1.0, x, y, offy=-2) assert_raises(Exception, f, 1.0, x, y, n=-2) assert_raises(Exception, f, 1.0, x, y, n=5) assert_raises(Exception, f, 1.0, x, y, lower=2) assert_raises(Exception, f, 1.0, x, y, a=np.zeros((2, 2), 'd', 'F')) @pytest.mark.parametrize("dtype", DTYPES) def test_gbmv(self, dtype): rng = np.random.default_rng(1234) n = 7 m = 5 kl = 1 ku = 2 # fake a banded matrix via toeplitz A = toeplitz(append(rng.random(kl+1), zeros(m-kl-1)), append(rng.random(ku+1), zeros(n-ku-1))) A = A.astype(dtype) Ab = zeros((kl+ku+1, n), dtype=dtype) # Form the banded storage Ab[2, :5] = A[0, 0] # diag Ab[1, 1:6] = A[0, 1] # sup1 Ab[0, 2:7] = A[0, 2] # sup2 Ab[3, :4] = A[1, 0] # sub1 x = rng.random(n).astype(dtype) y = rng.random(m).astype(dtype) alpha, beta = dtype(3), dtype(-5) func, = get_blas_funcs(('gbmv',), dtype=dtype) y1 = func(m=m, n=n, ku=ku, kl=kl, alpha=alpha, a=Ab, x=x, y=y, beta=beta) y2 = alpha * A.dot(x) + beta * y assert_array_almost_equal(y1, y2) y1 = func(m=m, n=n, ku=ku, kl=kl, alpha=alpha, a=Ab, x=y, y=x, beta=beta, trans=1) y2 = alpha * A.T.dot(y) + beta * x assert_array_almost_equal(y1, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_sbmv_hbmv(self, dtype): rng = np.random.default_rng(1234) n = 6 k = 2 A = zeros((n, n), dtype=dtype) Ab = zeros((k+1, n), dtype=dtype) # Form the array and its packed banded storage A[arange(n), arange(n)] = rng.random(n) for ind2 in range(1, k+1): temp = rng.random(n-ind2) A[arange(n-ind2), arange(ind2, n)] = temp Ab[-1-ind2, ind2:] = temp A = A.astype(dtype) if dtype in COMPLEX_DTYPES: A += A.conj().T func, = get_blas_funcs(('hbmv',), dtype=dtype) else: A += A.T func, = get_blas_funcs(('sbmv',), dtype=dtype) Ab[-1, :] = diag(A) x = rng.random(n).astype(dtype) y = rng.random(n).astype(dtype) alpha, beta = dtype(1.25), dtype(3) y1 = func(k=k, alpha=alpha, a=Ab, x=x, y=y, beta=beta) y2 = alpha * A.dot(x) + beta * y assert_array_almost_equal(y1, y2) @pytest.mark.parametrize("fname,dtype", [ *[('spmv', dtype) for dtype in REAL_DTYPES + COMPLEX_DTYPES], *[('hpmv', dtype) for dtype in COMPLEX_DTYPES], ]) def test_spmv_hpmv(self, fname, dtype): rng = np.random.default_rng(1234) n = 3 A = rng.random((n, n)).astype(dtype) if dtype in COMPLEX_DTYPES: A += rng.random((n, n))*1j A += A.T if fname == 'spmv' else A.conj().T c, r = tril_indices(n) Ap = A[r, c] x = rng.random(n).astype(dtype) y = rng.random(n).astype(dtype) xlong = arange(2*n).astype(dtype) ylong = ones(2*n).astype(dtype) alpha, beta = dtype(1.25), dtype(2) func, = get_blas_funcs((fname,), dtype=dtype) y1 = func(n=n, alpha=alpha, ap=Ap, x=x, y=y, beta=beta) y2 = alpha * A.dot(x) + beta * y assert_array_almost_equal(y1, y2) # Test inc and offsets y1 = func(n=n-1, alpha=alpha, beta=beta, x=xlong, y=ylong, ap=Ap, incx=2, incy=2, offx=n, offy=n) y2 = (alpha * A[:-1, :-1]).dot(xlong[3::2]) + beta * ylong[3::2] assert_array_almost_equal(y1[3::2], y2) assert_almost_equal(y1[4], ylong[4]) @pytest.mark.parametrize("fname,dtype", [ *[('spr', dtype) for dtype in REAL_DTYPES + COMPLEX_DTYPES], *[('hpr', dtype) for dtype in COMPLEX_DTYPES], ]) def test_spr_hpr(self, fname, dtype): rng = np.random.default_rng(1234) n = 3 A = rng.random((n, n)).astype(dtype) if dtype in COMPLEX_DTYPES: A += rng.random((n, n))*1j A += A.T if fname == 'spr' else A.conj().T c, r = tril_indices(n) Ap = A[r, c] x = rng.random(n).astype(dtype) alpha = np.finfo(dtype).dtype.type(2.5) if fname == 'hpr': func, = get_blas_funcs(('hpr',), dtype=dtype) y2 = alpha * x[:, None].dot(x[None, :].conj()) + A else: func, = get_blas_funcs(('spr',), dtype=dtype) y2 = alpha * x[:, None].dot(x[None, :]) + A y1 = func(n=n, alpha=alpha, ap=Ap, x=x) y1f = zeros((3, 3), dtype=dtype) y1f[r, c] = y1 y1f[c, r] = y1.conj() if fname == 'hpr' else y1 assert_array_almost_equal(y1f, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_spr2_hpr2(self, dtype): rng = np.random.default_rng(1234) n = 3 A = rng.random((n, n)).astype(dtype) if dtype in COMPLEX_DTYPES: A += rng.random((n, n))*1j A += A.conj().T func, = get_blas_funcs(('hpr2',), dtype=dtype) else: A += A.T func, = get_blas_funcs(('spr2',), dtype=dtype) c, r = tril_indices(n) Ap = A[r, c] x = rng.random(n).astype(dtype) y = rng.random(n).astype(dtype) alpha = dtype(2) u = alpha.conj() * x[:, None].dot(y[None, :].conj()) y2 = A + u + u.conj().T y1 = func(n=n, alpha=alpha, x=x, y=y, ap=Ap) y1f = zeros((3, 3), dtype=dtype) y1f[r, c] = y1 y1f[[1, 2, 2], [0, 0, 1]] = y1[[1, 3, 4]].conj() assert_array_almost_equal(y1f, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_tbmv(self, dtype): rng = np.random.default_rng(1234) n = 10 k = 3 x = rng.random(n).astype(dtype) A = zeros((n, n), dtype=dtype) # Banded upper triangular array for sup in range(k+1): A[arange(n-sup), arange(sup, n)] = rng.random(n-sup) # Add complex parts for c,z if dtype in COMPLEX_DTYPES: A[nonzero(A)] += 1j * rng.random((k+1)*n-(k*(k+1)//2)).astype(dtype) # Form the banded storage Ab = zeros((k+1, n), dtype=dtype) for row in range(k+1): Ab[-row-1, row:] = diag(A, k=row) func, = get_blas_funcs(('tbmv',), dtype=dtype) y1 = func(k=k, a=Ab, x=x) y2 = A.dot(x) assert_array_almost_equal(y1, y2) y1 = func(k=k, a=Ab, x=x, diag=1) A[arange(n), arange(n)] = dtype(1) y2 = A.dot(x) assert_array_almost_equal(y1, y2) y1 = func(k=k, a=Ab, x=x, diag=1, trans=1) y2 = A.T.dot(x) assert_array_almost_equal(y1, y2) y1 = func(k=k, a=Ab, x=x, diag=1, trans=2) y2 = A.conj().T.dot(x) assert_array_almost_equal(y1, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_tbsv(self, dtype): rng = np.random.default_rng(12345) n = 6 k = 3 x = rng.random(n).astype(dtype) A = zeros((n, n), dtype=dtype) # Banded upper triangular array for sup in range(k+1): A[arange(n-sup), arange(sup, n)] = rng.random(n-sup) # Add complex parts for c,z if dtype in COMPLEX_DTYPES: A[nonzero(A)] += 1j * rng.random((k+1)*n-(k*(k+1)//2)).astype(dtype) # Form the banded storage Ab = zeros((k+1, n), dtype=dtype) for row in range(k+1): Ab[-row-1, row:] = diag(A, k=row) func, = get_blas_funcs(('tbsv',), dtype=dtype) y1 = func(k=k, a=Ab, x=x) y2 = solve(A, x) assert_array_almost_equal(y1, y2) y1 = func(k=k, a=Ab, x=x, diag=1) A[arange(n), arange(n)] = dtype(1) y2 = solve(A, x) assert_array_almost_equal(y1, y2) y1 = func(k=k, a=Ab, x=x, diag=1, trans=1) y2 = solve(A.T, x) assert_array_almost_equal(y1, y2) y1 = func(k=k, a=Ab, x=x, diag=1, trans=2) y2 = solve(A.conj().T, x) assert_array_almost_equal(y1, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_tpmv(self, dtype): rng = np.random.default_rng(1234) n = 10 x = rng.random(n).astype(dtype) # Upper triangular array if dtype in COMPLEX_DTYPES: A = triu(rng.random((n, n)) + rng.random((n, n))*1j) else: A = triu(rng.random((n, n))) # Form the packed storage c, r = tril_indices(n) Ap = A[r, c] func, = get_blas_funcs(('tpmv',), dtype=dtype) y1 = func(n=n, ap=Ap, x=x) y2 = A.dot(x) assert_array_almost_equal(y1, y2) y1 = func(n=n, ap=Ap, x=x, diag=1) A[arange(n), arange(n)] = dtype(1) y2 = A.dot(x) assert_array_almost_equal(y1, y2) y1 = func(n=n, ap=Ap, x=x, diag=1, trans=1) y2 = A.T.dot(x) assert_array_almost_equal(y1, y2) y1 = func(n=n, ap=Ap, x=x, diag=1, trans=2) y2 = A.conj().T.dot(x) assert_array_almost_equal(y1, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_tpsv(self, dtype): rng = np.random.default_rng(1234) n = 10 x = rng.random(n).astype(dtype) # Upper triangular array if dtype in COMPLEX_DTYPES: A = triu(rng.random((n, n)) + rng.random((n, n))*1j) else: A = triu(rng.random((n, n))) A += eye(n) # Form the packed storage c, r = tril_indices(n) Ap = A[r, c] func, = get_blas_funcs(('tpsv',), dtype=dtype) y1 = func(n=n, ap=Ap, x=x) y2 = solve(A, x) assert_array_almost_equal(y1, y2) y1 = func(n=n, ap=Ap, x=x, diag=1) A[arange(n), arange(n)] = dtype(1) y2 = solve(A, x) assert_array_almost_equal(y1, y2) y1 = func(n=n, ap=Ap, x=x, diag=1, trans=1) y2 = solve(A.T, x) assert_array_almost_equal(y1, y2) y1 = func(n=n, ap=Ap, x=x, diag=1, trans=2) y2 = solve(A.conj().T, x) assert_array_almost_equal(y1, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_trmv(self, dtype): rng = np.random.default_rng(1234) n = 3 A = (rng.random((n, n))+eye(n)).astype(dtype) x = rng.random(3).astype(dtype) func, = get_blas_funcs(('trmv',), dtype=dtype) y1 = func(a=A, x=x) y2 = triu(A).dot(x) assert_array_almost_equal(y1, y2) y1 = func(a=A, x=x, diag=1) A[arange(n), arange(n)] = dtype(1) y2 = triu(A).dot(x) assert_array_almost_equal(y1, y2) y1 = func(a=A, x=x, diag=1, trans=1) y2 = triu(A).T.dot(x) assert_array_almost_equal(y1, y2) y1 = func(a=A, x=x, diag=1, trans=2) y2 = triu(A).conj().T.dot(x) assert_array_almost_equal(y1, y2) @pytest.mark.parametrize("dtype", DTYPES) def test_trsv(self, dtype): rng = np.random.default_rng(1234) n = 15 A = (rng.random((n, n))+eye(n)).astype(dtype) x = rng.random(n).astype(dtype) func, = get_blas_funcs(('trsv',), dtype=dtype) y1 = func(a=A, x=x) y2 = solve(triu(A), x) assert_array_almost_equal(y1, y2) y1 = func(a=A, x=x, lower=1) y2 = solve(tril(A), x) assert_array_almost_equal(y1, y2) y1 = func(a=A, x=x, diag=1) A[arange(n), arange(n)] = dtype(1) y2 = solve(triu(A), x) assert_array_almost_equal(y1, y2) y1 = func(a=A, x=x, diag=1, trans=1) y2 = solve(triu(A).T, x) assert_array_almost_equal(y1, y2) y1 = func(a=A, x=x, diag=1, trans=2) y2 = solve(triu(A).conj().T, x) assert_array_almost_equal(y1, y2) class TestFBLAS3Simple: @parametrize_blas(fblas, "gemm", "sdcz") def test_gemm(self, f, dtype): assert_array_almost_equal(f(3, [3], [-4]), [[-36]]) assert_array_almost_equal(f(3, [3], [-4], 3, [5]), [-21]) if dtype in COMPLEX_DTYPES: assert_array_almost_equal(f(3j, [3-4j], [-4]), [[-48-36j]]) assert_array_almost_equal(f(3j, [3-4j], [-4], 3, [5j]), [-48-21j]) class TestBLAS3Symm: def setup_method(self): self.a = np.array([[1., 2.], [0., 1.]]) self.b = np.array([[1., 0., 3.], [0., -1., 2.]]) self.c = np.ones((2, 3)) self.t = np.array([[2., -1., 8.], [3., 0., 9.]]) @parametrize_blas(fblas, "symm", "sdcz") def test_symm(self, f, dtype): res = f(a=self.a, b=self.b, c=self.c, alpha=1., beta=1.) assert_array_almost_equal(res, self.t) res = f(a=self.a.T, b=self.b, lower=1, c=self.c, alpha=1., beta=1.) assert_array_almost_equal(res, self.t) res = f(a=self.a, b=self.b.T, side=1, c=self.c.T, alpha=1., beta=1.) assert_array_almost_equal(res, self.t.T) @parametrize_blas(fblas, "symm", "sdcz") def test_symm_wrong_side(self, f, dtype): """`side=1` means C <- B*A, hence shapes of A and B are to be compatible. Otherwise, f2py exception is raised. """ # FIXME narrow down to _fblas.error with pytest.raises(Exception): f(a=self.a, b=self.b, alpha=1, side=1) @parametrize_blas(fblas, "symm", "sdcz") def test_symm_wrong_uplo(self, f, dtype): """SYMM only considers the upper/lower part of A. Hence setting wrong value for `lower` (default is lower=0, meaning upper triangle) gives a wrong result. """ res = f(a=self.a, b=self.b, c=self.c, alpha=1., beta=1.) assert np.allclose(res, self.t) res = f(a=self.a, b=self.b, lower=1, c=self.c, alpha=1., beta=1.) assert not np.allclose(res, self.t) class TestBLAS3Syrk: def setup_method(self): self.a = np.array([[1., 0.], [0., -2.], [2., 3.]]) self.t = np.array([[1., 0., 2.], [0., 4., -6.], [2., -6., 13.]]) self.tt = np.array([[5., 6.], [6., 13.]]) @parametrize_blas(fblas, "syrk", "sdcz") def test_syrk(self, f, dtype): c = f(a=self.a, alpha=1.) assert_array_almost_equal(np.triu(c), np.triu(self.t)) c = f(a=self.a, alpha=1., lower=1) assert_array_almost_equal(np.tril(c), np.tril(self.t)) c0 = np.ones(self.t.shape) c = f(a=self.a, alpha=1., beta=1., c=c0) assert_array_almost_equal(np.triu(c), np.triu(self.t+c0)) c = f(a=self.a, alpha=1., trans=1) assert_array_almost_equal(np.triu(c), np.triu(self.tt)) # prints '0-th dimension must be fixed to 3 but got 5', # FIXME: suppress? @parametrize_blas(fblas, "syrk", "sdcz") def test_syrk_wrong_c(self, f, dtype): # FIXME narrow down to _fblas.error with pytest.raises(Exception): f(a=self.a, alpha=1., c=np.ones((5, 8))) # if C is supplied, it must have compatible dimensions class TestBLAS3Syr2k: def setup_method(self): self.a = np.array([[1., 0.], [0., -2.], [2., 3.]]) self.b = np.array([[0., 1.], [1., 0.], [0, 1.]]) self.t = np.array([[0., -1., 3.], [-1., 0., 0.], [3., 0., 6.]]) self.tt = np.array([[0., 1.], [1., 6]]) @parametrize_blas(fblas, "syr2k", "sdcz") def test_syr2k(self, f, dtype): c = f(a=self.a, b=self.b, alpha=1.) assert_array_almost_equal(np.triu(c), np.triu(self.t)) c = f(a=self.a, b=self.b, alpha=1., lower=1) assert_array_almost_equal(np.tril(c), np.tril(self.t)) c0 = np.ones(self.t.shape) c = f(a=self.a, b=self.b, alpha=1., beta=1., c=c0) assert_array_almost_equal(np.triu(c), np.triu(self.t+c0)) c = f(a=self.a, b=self.b, alpha=1., trans=1) assert_array_almost_equal(np.triu(c), np.triu(self.tt)) # prints '0-th dimension must be fixed to 3 but got 5', FIXME: suppress? @parametrize_blas(fblas, "syr2k", "sdcz") def test_syr2k_wrong_c(self, f, dtype): with pytest.raises(Exception): f(a=self.a, b=self.b, alpha=1., c=np.zeros((15, 8))) # if C is supplied, it must have compatible dimensions class TestSyHe: """Quick and simple tests for (zc)-symm, syrk, syr2k.""" def setup_method(self): self.sigma_y = np.array([[0., -1.j], [1.j, 0.]]) @parametrize_blas(fblas, "symm", "zc") def test_symm(self, f, dtype): # NB: a is symmetric w/upper diag of ONLY res = f(a=self.sigma_y, b=self.sigma_y, alpha=1.) assert_array_almost_equal(np.triu(res), np.diag([1, -1])) @parametrize_blas(fblas, "hemm", "zc") def test_hemm(self, f, dtype): # NB: a is hermitian w/upper diag of ONLY res = f(a=self.sigma_y, b=self.sigma_y, alpha=1.) assert_array_almost_equal(np.triu(res), np.diag([1, 1])) @parametrize_blas(fblas, "syrk", "zc") def test_syrk(self, f, dtype): res = f(a=self.sigma_y, alpha=1.) assert_array_almost_equal(np.triu(res), np.diag([-1, -1])) @parametrize_blas(fblas, "herk", "zc") def test_herk(self, f, dtype): res = f(a=self.sigma_y, alpha=1.) assert_array_almost_equal(np.triu(res), np.diag([1, 1])) @parametrize_blas(fblas, "syr2k", "zc") def test_syr2k_zr(self, f, dtype): res = f(a=self.sigma_y, b=self.sigma_y, alpha=1.) assert_array_almost_equal(np.triu(res), 2.*np.diag([-1, -1])) @parametrize_blas(fblas, "her2k", "zc") def test_her2k_zr(self, f, dtype): res = f(a=self.sigma_y, b=self.sigma_y, alpha=1.) assert_array_almost_equal(np.triu(res), 2.*np.diag([1, 1])) class TestTRMM: """Quick and simple tests for *trmm.""" def setup_method(self): self.a = np.array([[1., 2., ], [-2., 1.]]) self.b = np.array([[3., 4., -1.], [5., 6., -2.]]) self.a2 = np.array([[1, 1, 2, 3], [0, 1, 4, 5], [0, 0, 1, 6], [0, 0, 0, 1]], order="f") self.b2 = np.array([[1, 4], [2, 5], [3, 6], [7, 8], [9, 10]], order="f") @pytest.mark.parametrize("dtype", DTYPES) def test_side(self, dtype): trmm = get_blas_funcs("trmm", dtype=dtype) # Provide large A array that works for side=1 but not 0 (see gh-10841) assert_raises(Exception, trmm, 1.0, self.a2, self.b2) res = trmm(1.0, self.a2.astype(dtype), self.b2.astype(dtype), side=1) k = self.b2.shape[1] assert_allclose(res, self.b2 @ self.a2[:k, :k], rtol=0., atol=100*np.finfo(dtype).eps) @parametrize_blas(fblas, "trmm", "sdcz") def test_ab(self, f, dtype): result = f(1., self.a, self.b) # default a is upper triangular expected = np.array([[13., 16., -5.], [ 5., 6., -2.]]) assert_array_almost_equal(result, expected) @parametrize_blas(fblas, "trmm", "sdcz") def test_ab_lower(self, f, dtype): result = f(1., self.a, self.b, lower=True) expected = np.array([[ 3., 4., -1.], [-1., -2., 0.]]) # now a is lower triangular assert_array_almost_equal(result, expected) @parametrize_blas(fblas, "trmm", "sdcz") def test_b_overwrites(self, f, dtype): # BLAS *trmm modifies B argument in-place. # Here the default is to copy, but this can be overridden b = self.b.astype(dtype) for overwr in [True, False]: bcopy = b.copy() result = f(1., self.a, bcopy, overwrite_b=overwr) # C-contiguous arrays are copied assert not bcopy.flags.f_contiguous assert not np.may_share_memory(bcopy, result) assert_equal(bcopy, b) bcopy = np.asfortranarray(b.copy()) # or just transpose it result = f(1., self.a, bcopy, overwrite_b=True) assert bcopy.flags.f_contiguous assert np.may_share_memory(bcopy, result) assert_array_almost_equal(bcopy, result) @pytest.mark.parametrize("dtype", DTYPES) def test_trsm(dtype): rng = np.random.default_rng(1234) tol = np.finfo(dtype).eps*1000 func, = get_blas_funcs(('trsm',), dtype=dtype) # Test protection against size mismatches A = rng.random((4, 5)).astype(dtype) B = rng.random((4, 4)).astype(dtype) alpha = dtype(1) assert_raises(Exception, func, alpha, A, B) assert_raises(Exception, func, alpha, A.T, B) n = 8 m = 7 alpha = dtype(-2.5) if dtype in COMPLEX_DTYPES: A = (rng.random((m, m)) + rng.random((m, m))*1j) + eye(m) else: A = rng.random((m, m)) + eye(m) A = A.astype(dtype) Au = triu(A) Al = tril(A) B1 = rng.random((m, n)).astype(dtype) B2 = rng.random((n, m)).astype(dtype) x1 = func(alpha=alpha, a=A, b=B1) assert_equal(B1.shape, x1.shape) x2 = solve(Au, alpha*B1) assert_allclose(x1, x2, atol=tol) x1 = func(alpha=alpha, a=A, b=B1, trans_a=1) x2 = solve(Au.T, alpha*B1) assert_allclose(x1, x2, atol=tol) x1 = func(alpha=alpha, a=A, b=B1, trans_a=2) x2 = solve(Au.conj().T, alpha*B1) assert_allclose(x1, x2, atol=tol) x1 = func(alpha=alpha, a=A, b=B1, diag=1) Au[arange(m), arange(m)] = dtype(1) x2 = solve(Au, alpha*B1) assert_allclose(x1, x2, atol=tol) x1 = func(alpha=alpha, a=A, b=B2, diag=1, side=1) x2 = solve(Au.conj().T, alpha*B2.conj().T) assert_allclose(x1, x2.conj().T, atol=tol) x1 = func(alpha=alpha, a=A, b=B2, diag=1, side=1, lower=1) Al[arange(m), arange(m)] = dtype(1) x2 = solve(Al.conj().T, alpha*B2.conj().T) assert_allclose(x1, x2.conj().T, atol=tol) @pytest.mark.xfail(run=False, reason="gh-16930") def test_gh_169309(): x = np.repeat(10, 9) actual = scipy.linalg.blas.dnrm2(x, 5, 3, -1) expected = math.sqrt(500) assert_allclose(actual, expected) def test_dnrm2_neg_incx(): # check that dnrm2(..., incx < 0) raises # XXX: remove the test after the lowest supported BLAS implements # negative incx (new in LAPACK 3.10) x = np.repeat(10, 9) incx = -1 with assert_raises(fblas.__fblas_error): scipy.linalg.blas.dnrm2(x, 5, 3, incx)