def test_spr_hpr(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES+COMPLEX_DTYPES):
            n = 3
            A = rand(n, n).astype(dtype)
            if ind > 1:
                A += rand(n, n)*1j
            A = A.astype(dtype)
            A = A + A.T if ind < 4 else A + A.conj().T
            c, r = tril_indices(n)
            Ap = A[r, c]
            x = rand(n).astype(dtype)
            alpha = (DTYPES+COMPLEX_DTYPES)[mod(ind, 4)](2.5)

            if ind > 3:
                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 ind > 3 else y1
            assert_array_almost_equal(y1f, y2)

    def test_spmv_hpmv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES+COMPLEX_DTYPES):
            n = 3
            A = rand(n, n).astype(dtype)
            if ind > 1:
                A += rand(n, n)*1j
            A = A.astype(dtype)
            A = A + A.T if ind < 4 else A + A.conj().T
            c, r = tril_indices(n)
            Ap = A[r, c]
            x = rand(n).astype(dtype)
            y = rand(n).astype(dtype)
            xlong = arange(2*n).astype(dtype)
            ylong = ones(2*n).astype(dtype)
            alpha, beta = dtype(1.25), dtype(2)

            if ind > 3:
                func, = get_blas_funcs(('hpmv',), dtype=dtype)
            else:
                func, = get_blas_funcs(('spmv',), 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])

def _have_blas_gemm():
    try:
        linalg.get_blas_funcs(["gemm"])
        return True
    except (AttributeError, ValueError):
        warnings.warn("Could not import BLAS, falling back to np.dot")
        return False

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(f1.module_name, 'cblas')
    assert_equal(f2.typecode, 'z')
    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')

    def test_sbmv_hbmv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            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)] = rand(n)
            for ind2 in range(1, k+1):
                temp = rand(n-ind2)
                A[arange(n-ind2), arange(ind2, n)] = temp
                Ab[-1-ind2, ind2:] = temp
            A = A.astype(dtype)
            A = A + A.T if ind < 2 else A + A.conj().T
            Ab[-1, :] = diag(A)
            x = rand(n).astype(dtype)
            y = rand(n).astype(dtype)
            alpha, beta = dtype(1.25), dtype(3)

            if ind > 1:
                func, = get_blas_funcs(('hbmv',), dtype=dtype)
            else:
                func, = get_blas_funcs(('sbmv',), dtype=dtype)
            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)

    def test_spr2_hpr2(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 3
            A = rand(n, n).astype(dtype)
            if ind > 1:
                A += rand(n, n)*1j
            A = A.astype(dtype)
            A = A + A.T if ind < 2 else A + A.conj().T
            c, r = tril_indices(n)
            Ap = A[r, c]
            x = rand(n).astype(dtype)
            y = rand(n).astype(dtype)
            alpha = dtype(2)

            if ind > 1:
                func, = get_blas_funcs(('hpr2',), dtype=dtype)
            else:
                func, = get_blas_funcs(('spr2',), dtype=dtype)

            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)

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 test_fast_dot():
    """Check fast dot blas wrapper function"""
    rng = np.random.RandomState(42)
    A = rng.random_sample([2, 10])
    B = rng.random_sample([2, 10])

    try:
        linalg.get_blas_funcs('gemm')
        has_blas = True
    except AttributeError, ValueError:
        has_blas = False

def py_rect_maxvol(A, tol = 1., maxK = None, min_add_K = None, minK = None, start_maxvol_iters = 10, identity_submatrix = True):
    """Python implementation of rectangular 2-volume maximization. For information see :py:func:`rect_maxvol` function"""
    # tol2 - square of parameter tol
    tol2 = tol**2
    # N - number of rows, r - number of columns of matrix A
    N, r = A.shape
    if N <= r:
        return np.arange(N, dtype = np.int32), np.eye(N, dtype = A.dtype)
    if maxK is None or maxK > N:
        maxK = N
    if maxK < r:
        maxK = r
    if minK is None or minK < r:
        minK = r
    if minK > N:
        minK = N
    if min_add_K is not None:
        minK = max(minK, r + min_add_K) 
    if minK > maxK:
        minK = maxK
        #raise ValueError('minK value cannot be greater than maxK value')
    index = np.zeros(N, dtype = np.int32)
    chosen = np.ones(N)
    tmp_index, C = py_maxvol(A, tol = 1, max_iters = start_maxvol_iters)
    index[:r] = tmp_index
    chosen[tmp_index] = 0
    C = np.asfortranarray(C)
    # compute square 2-norms of each row in matrix C
    row_norm_sqr = np.array([chosen[i]*np.linalg.norm(C[i], 2)**2 for i in xrange(N)])
    # find maximum value in row_norm_sqr
    i = np.argmax(row_norm_sqr)
    K = r
    # set cgeru or zgeru for complex numbers and dger or sger for float numbers
    try:
        ger = get_blas_funcs('geru', [C])
    except:
        ger = get_blas_funcs('ger', [C])
    while (row_norm_sqr[i] > tol2 and K < maxK) or K < minK:
        # add i to index and recompute C and square norms of each row by SVM-formula
        index[K] = i
        chosen[i] = 0
        c = C[i].copy()
        v = C.dot(c.conj())
        l = 1.0/(1+v[i])
        ger(-l,v,c,a=C,overwrite_a=1)
        C = np.hstack([C, l*v.reshape(-1,1)])
        row_norm_sqr -= (l*v*v.conj()).real
        row_norm_sqr *= chosen
        # find maximum value in row_norm_sqr
        i = row_norm_sqr.argmax()
        K += 1
    if identity_submatrix:
        C[index[:K]] = np.eye(K, dtype = C.dtype)
    return index[:K].copy(), C

def test_inplace_swap_column():
    X = np.array([[0, 3, 0],
                  [2, 4, 0],
                  [0, 0, 0],
                  [9, 8, 7],
                  [4, 0, 5]], dtype=np.float64)
    X_csr = sp.csr_matrix(X)
    X_csc = sp.csc_matrix(X)

    swap = linalg.get_blas_funcs(('swap',), (X,))
    swap = swap[0]
    X[:, 0], X[:, -1] = swap(X[:, 0], X[:, -1])
    inplace_swap_column(X_csr, 0, -1)
    inplace_swap_column(X_csc, 0, -1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())

    X[:, 0], X[:, 1] = swap(X[:, 0], X[:, 1])
    inplace_swap_column(X_csr, 0, 1)
    inplace_swap_column(X_csc, 0, 1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())
    assert_raises(TypeError, inplace_swap_column, X_csr.tolil())

    X = np.array([[0, 3, 0],
                  [2, 4, 0],
                  [0, 0, 0],
                  [9, 8, 7],
                  [4, 0, 5]], dtype=np.float32)
    X_csr = sp.csr_matrix(X)
    X_csc = sp.csc_matrix(X)
    swap = linalg.get_blas_funcs(('swap',), (X,))
    swap = swap[0]
    X[:, 0], X[:, -1] = swap(X[:, 0], X[:, -1])
    inplace_swap_column(X_csr, 0, -1)
    inplace_swap_column(X_csc, 0, -1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())
    X[:, 0], X[:, 1] = swap(X[:, 0], X[:, 1])
    inplace_swap_column(X_csr, 0, 1)
    inplace_swap_column(X_csc, 0, 1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())
    assert_raises(TypeError, inplace_swap_column, X_csr.tolil())

def norm2(q):
    """
    Compute the euclidean norm of an array ``q`` by calling the BLAS routine
    """
    q = np.asarray(q)
    nrm2 = get_blas_funcs('nrm2', dtype=q.dtype)
    return nrm2(q)

    def test_trsv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 15
            A = (rand(n, n)+eye(n)).astype(dtype)
            x = rand(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)

    def test_tpsv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 10
            x = rand(n).astype(dtype)
            # Upper triangular array
            A = triu(rand(n, n)) if ind < 2 else triu(rand(n, n)+rand(n, n)*1j)
            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)

def _update_dict_slow(X, A, B, G, Q, Q_idx, idx, fit_intercept,
                      components_range, norm, impute=True):
    Q_idx = Q[:, idx]

    if impute:
        old_sub_G = Q_idx.dot(Q_idx.T)

    ger, = linalg.get_blas_funcs(('ger',), (A, Q_idx))
    R = B[:, idx] - np.dot(Q_idx.T, A).T

    # norm = np.sqrt(np.sum(Q_idx ** 2, axis=1))
    norm = np.sqrt(np.sum(Q ** 2, axis=1))
    # Intercept on first column
    for j in components_range:
        ger(1.0, A[j], Q_idx[j], a=R, overwrite_a=True)
        Q_idx[j] = R[j] / A[j, j]
        # new_norm = np.sqrt(np.sum(Q_idx[j] ** 2))
        # if new_norm > norm[j]:
        #     Q_idx[j] /= new_norm / norm[j]
        Q[j, idx] = Q_idx[j]
        new_norm = np.sqrt(np.sum(Q[j] ** 2))
        if new_norm > 1:
            Q_idx[j] /= new_norm
            Q[j] /= new_norm

        ger(-1.0, A[j], Q_idx[j], a=R, overwrite_a=True)

    Q[:, idx] = Q_idx

    if impute:
        G += Q_idx.dot(Q_idx.T) - old_sub_G

def cool_syrk(fact, X): 
    syrk = get_blas_funcs("syrk", [X])
    R = syrk(fact, X)
    d = np.diag(R).copy()
    size = mat_to_upper_F(R)
    R.resize([size,])
    return R,d

    def test_gbmv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 7
            m = 5
            kl = 1
            ku = 2
            # fake a banded matrix via toeplitz
            A = toeplitz(append(rand(kl+1), zeros(m-kl-1)),
                         append(rand(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 = rand(n).astype(dtype)
            y = rand(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)

    def _update_dict_slow(self, subset, D_range):
        """Update dictionary from statistic
        Parameters
        ----------
        subset: ndarray (len_subset),
            Mask used on X

        """
        D_subset = self.D_[:, subset]

        if self.projection == "full":
            norm = enet_norm(self.D_, self.l1_ratio)
        else:
            norm = enet_norm(D_subset, self.l1_ratio)
        R = self.B_[:, subset] - np.dot(D_subset.T, self.A_).T

        ger, = linalg.get_blas_funcs(("ger",), (self.A_, D_subset))
        for k in D_range:
            ger(1.0, self.A_[k], D_subset[k], a=R, overwrite_a=True)
            # R += np.dot(stat.A[:, j].reshape(n_components, 1),
            D_subset[k] = R[k] / (self.A_[k, k])
            if self.projection == "full":
                self.D_[k][subset] = D_subset[k]
                self.D_[k] = enet_projection(self.D_[k], norm[k], self.l1_ratio)
                D_subset[k] = self.D_[k][subset]
            else:
                D_subset[k] = enet_projection(D_subset[k], norm[k], self.l1_ratio)
            ger(-1.0, self.A_[k], D_subset[k], a=R, overwrite_a=True)
            # R -= np.dot(stat.A[:, j].reshape(n_components, 1),
        if self.projection == "partial":
            self.D_[:, subset] = D_subset

    def _solve(v, alpha, cs, ds):
        """Evaluate w = M^-1 v"""
        if len(cs) == 0:
            return v/alpha

        # (B + C D^H)^-1 = B^-1 - B^-1 C (I + D^H B^-1 C)^-1 D^H B^-1

        axpy, dotc = get_blas_funcs(['axpy', 'dotc'], cs[:1] + [v])

        c0 = cs[0]
        A = alpha * np.identity(len(cs), dtype=c0.dtype)
        for i, d in enumerate(ds):
            for j, c in enumerate(cs):
                A[i,j] += dotc(d, c)

        q = np.zeros(len(cs), dtype=c0.dtype)
        for j, d in enumerate(ds):
            q[j] = dotc(d, v)
        q /= alpha
        q = solve(A, q)

        w = v/alpha
        for c, qc in zip(cs, q):
            w = axpy(c, w, w.size, -qc)

        return w

def dot(A, B, out=None):
    """
    A drop in replacement for numpy.dot.
    Computes A.B optimized using fblas calls.
   """
    import scipy.linalg as sp
    gemm = sp.get_blas_funcs('gemm', arrays=(A,B))

    if out is None:
        lda, x, y, ldb = A.shape + B.shape
        if x != y:
            raise ValueError("matrices are not aligned")
        dtype = np.max([x.dtype for x in (A, B)])
        out = np.empty((lda, ldb), dtype, order='F')

    if A.flags.c_contiguous and B.flags.c_contiguous:
        gemm(alpha=1., a=A.T, b=B.T, trans_a=True, trans_b=True, c=out, overwrite_c=True)
    if A.flags.c_contiguous and B.flags.f_contiguous:
        gemm(alpha=1., a=A.T, b=B, trans_a=True, c=out, overwrite_c=True)
    if A.flags.f_contiguous and B.flags.c_contiguous:
        gemm(alpha=1., a=A, b=B.T, trans_b=True, c=out, overwrite_c=True)
    if A.flags.f_contiguous and B.flags.f_contiguous:
        gemm(alpha=1., a=A, b=B, c=out, overwrite_c=True)

    return out

def _mixed_norm_solver_bcd(M, G, alpha, lipschitz_constant, maxit=200,
                           tol=1e-8, verbose=None, init=None, n_orient=1,
                           dgap_freq=10):
    """Solve L21 inverse problem with block coordinate descent."""
    n_sensors, n_times = M.shape
    n_sensors, n_sources = G.shape
    n_positions = n_sources // n_orient

    if init is None:
        X = np.zeros((n_sources, n_times))
        R = M.copy()
    else:
        X = init
        R = M - np.dot(G, X)

    E = []  # track primal objective function
    highest_d_obj = - np.inf
    active_set = np.zeros(n_sources, dtype=np.bool)  # start with full AS

    alpha_lc = alpha / lipschitz_constant

    # First make G fortran for faster access to blocks of columns
    G = np.asfortranarray(G)
    # It is better to call gemm here
    # so it is called only once
    gemm = linalg.get_blas_funcs("gemm", [R.T, G[:, 0:n_orient]])
    one_ovr_lc = 1. / lipschitz_constant

    # assert that all the multiplied matrices are fortran contiguous
    assert X.T.flags.f_contiguous
    assert R.T.flags.f_contiguous
    assert G.flags.f_contiguous
    # storing list of contiguous arrays
    list_G_j_c = []
    for j in range(n_positions):
        idx = slice(j * n_orient, (j + 1) * n_orient)
        list_G_j_c.append(np.ascontiguousarray(G[:, idx]))

    for i in range(maxit):
        _bcd(G, X, R, active_set, one_ovr_lc, n_orient, n_positions,
             alpha_lc, gemm, list_G_j_c)

        if (i + 1) % dgap_freq == 0:
            _, p_obj, d_obj, _ = dgap_l21(M, G, X[active_set], active_set,
                                          alpha, n_orient)
            highest_d_obj = max(d_obj, highest_d_obj)
            gap = p_obj - highest_d_obj
            E.append(p_obj)
            logger.debug("Iteration %d :: p_obj %f :: dgap %f :: n_active %d" %
                         (i + 1, p_obj, gap, np.sum(active_set) / n_orient))

            if gap < tol:
                logger.debug('Convergence reached ! (gap: %s < %s)'
                             % (gap, tol))
                break

    X = X[active_set]

    return X, active_set, E