Python scipy.linalg.lu_factor() Examples

def build_krylov_space(frequency, r, side, a, b):

    if frequency == np.inf or frequency.real == np.inf:
        approx_type = 'partial_realisation'
        lu_a = a
    else:
        approx_type = 'Pade'
        lu_a = lu_factor(frequency, a)

    try:
        nu = b.shape[1]
    except IndexError:
        nu = 1

    if nu == 1:
        krylov_function = construct_krylov
    else:
        krylov_function = construct_mimo_krylov

    v = krylov_function(r, lu_a, b, approx_type, side)

    return v 

def lu_factor(sigma, A):
    """
    LU Factorisation wrapper of:

    .. math:: LU = (\sigma \mathbf{I} - \mathbf{A})

    In the case of ``A`` being a sparse matrix, the sparse methods in scipy are employed

    Args:
        sigma (float): Expansion frequency
        A (csc_matrix or np.ndarray): Dynamics matrix

    Returns:
        tuple or SuperLU: tuple (dense) or SuperLU (sparse) objects containing the LU factorisation
    """
    n = A.shape[0]
    if type(A) == libsp.csc_matrix:
        return scsp.linalg.splu(sigma * scsp.identity(n, dtype=complex, format='csc') - A)
    else:
        return sclalg.lu_factor(sigma * np.eye(n) - A) 

def __init__(self, inv_array, inv_lu_and_piv, inv_lu_transposed):
        """
        Args:
            inv_array (array): 2D array specifying inverse matrix entries.
            inv_lu_and_piv (Tuple[array, array]): Pivoted LU factorisation
                represented as a tuple with first element a 2D array containing
                the lower and upper triangular factors (with the unit diagonal
                of the lower triangular factor not stored) and the second
                element a 1D array containing the pivot indices. Corresponds
                to the output of `scipy.linalg.lu_factor` and input to
                `scipy.linalg.lu_solve`.
            inv_lu_transposed (bool): Whether LU factorisation is of inverse of
                array or transpose of inverse of array.
        """
        super().__init__(inv_array.shape)
        self._inv_array = inv_array
        self._inv_lu_and_piv = inv_lu_and_piv
        self._inv_lu_transposed = inv_lu_transposed 

def __init__(self, M):
        self.M_lu = lu_factor(M)
        self.shape = M.shape
        self.dtype = M.dtype 

def lu_decomp(A):
    LOG.debug(f"Compute LU decomposition of {A}.")
    return sl.lu_factor(A.full_matrix()) 

def lu_solve_AATI(A, rho, b, lu, piv, check_finite=True):
    r"""Solve the linear system :math:`(A A^T + \rho I)\mathbf{x} = \mathbf{b}`
    or :math:`(A A^T + \rho I)X = B` using :func:`scipy.linalg.lu_solve`.

    Parameters
    ----------
    A : array_like
      Matrix :math:`A`
    rho : float
      Scalar :math:`\rho`
    b : array_like
      Vector :math:`\mathbf{b}` or matrix :math:`B`
    lu : array_like
      Matrix containing U in its upper triangle, and L in its lower triangle,
      as returned by :func:`scipy.linalg.lu_factor`
    piv : array_like
      Pivot indices representing the permutation matrix P, as returned by
      :func:`scipy.linalg.lu_factor`
    check_finite : bool, optional (default False)
      Flag indicating whether the input array should be checked for Inf
      and NaN values

    Returns
    -------
    x : ndarray
      Solution to the linear system
    """

    N, M = A.shape
    if N >= M:
        x = (b - linalg.lu_solve((lu, piv), b.dot(A).T,
                                 check_finite=check_finite).T.dot(A.T)) / rho
    else:
        x = linalg.lu_solve((lu, piv), b.T, check_finite=check_finite).T
    return x 

def lu_solve_ATAI(A, rho, b, lu, piv, check_finite=True):
    r"""Solve the linear system :math:`(A^T A + \rho I)\mathbf{x} = \mathbf{b}`
    or :math:`(A^T A + \rho I)X = B` using :func:`scipy.linalg.lu_solve`.

    Parameters
    ----------
    A : array_like
      Matrix :math:`A`
    rho : float
      Scalar :math:`\rho`
    b : array_like
      Vector :math:`\mathbf{b}` or matrix :math:`B`
    lu : array_like
      Matrix containing U in its upper triangle, and L in its lower triangle,
      as returned by :func:`scipy.linalg.lu_factor`
    piv : array_like
      Pivot indices representing the permutation matrix P, as returned by
      :func:`scipy.linalg.lu_factor`
    check_finite : bool, optional (default False)
      Flag indicating whether the input array should be checked for Inf
      and NaN values

    Returns
    -------
    x : ndarray
      Solution to the linear system
    """

    N, M = A.shape
    if N >= M:
        x = linalg.lu_solve((lu, piv), b, check_finite=check_finite)
    else:
        x = (b - A.T.dot(linalg.lu_solve((lu, piv), A.dot(b), 1,
                                         check_finite=check_finite))) / rho
    return x 

def lu_factor(A, rho, check_finite=True):
    r"""Compute LU factorisation of either :math:`A^T A + \rho I` or
    :math:`A A^T + \rho I`, depending on which matrix is smaller.

    Parameters
    ----------
    A : array_like
      Array :math:`A`
    rho : float
      Scalar :math:`\rho`
    check_finite : bool, optional (default False)
      Flag indicating whether the input array should be checked for Inf
      and NaN values

    Returns
    -------
    lu : ndarray
      Matrix containing U in its upper triangle, and L in its lower
      triangle, as returned by :func:`scipy.linalg.lu_factor`
    piv : ndarray
      Pivot indices representing the permutation matrix P, as returned
      by :func:`scipy.linalg.lu_factor`
    """

    N, M = A.shape
    # If N < M it is cheaper to factorise A*A^T + rho*I and then use the
    # matrix inversion lemma to compute the inverse of A^T*A + rho*I
    if N >= M:
        lu, piv = linalg.lu_factor(A.T.dot(A) +
                                   rho * np.identity(M, dtype=A.dtype),
                                   check_finite=check_finite)
    else:
        lu, piv = linalg.lu_factor(A.dot(A.T) +
                                   rho * np.identity(N, dtype=A.dtype),
                                   check_finite=check_finite)
    return lu, piv 

def lu_solve_AATI(A, rho, b, lu, piv):
    r"""
    Solve the linear system :math:`(A A^T + \rho I)\mathbf{x} = \mathbf{b}`
    or :math:`(A A^T + \rho I)X = B` using :func:`scipy.linalg.lu_solve`.

    Parameters
    ----------
    A : array_like
      Matrix :math:`A`
    rho : float
      Scalar :math:`\rho`
    b : array_like
      Vector :math:`\mathbf{b}` or matrix :math:`B`
    lu : array_like
      Matrix containing U in its upper triangle, and L in its lower triangle,
      as returned by :func:`scipy.linalg.lu_factor`
    piv : array_like
      Pivot indices representing the permutation matrix P, as returned by
      :func:`scipy.linalg.lu_factor`

    Returns
    -------
    x : ndarray
      Solution to the linear system.
    """

    N, M = A.shape
    if N >= M:
        x = (b - linalg.lu_solve((lu, piv), b.dot(A).T).T.dot(A.T)) / rho
    else:
        x = linalg.lu_solve((lu, piv), b.T).T
    return x 

def lu_solve_ATAI(A, rho, b, lu, piv):
    r"""
    Solve the linear system :math:`(A^T A + \rho I)\mathbf{x} = \mathbf{b}`
    or :math:`(A^T A + \rho I)X = B` using :func:`scipy.linalg.lu_solve`.

    Parameters
    ----------
    A : array_like
      Matrix :math:`A`
    rho : float
      Scalar :math:`\rho`
    b : array_like
      Vector :math:`\mathbf{b}` or matrix :math:`B`
    lu : array_like
      Matrix containing U in its upper triangle, and L in its lower triangle,
      as returned by :func:`scipy.linalg.lu_factor`
    piv : array_like
      Pivot indices representing the permutation matrix P, as returned by
      :func:`scipy.linalg.lu_factor`

    Returns
    -------
    x : ndarray
      Solution to the linear system.
    """

    N, M = A.shape
    if N >= M:
        x = linalg.lu_solve((lu, piv), b)
    else:
        x = (b - A.T.dot(linalg.lu_solve((lu, piv), A.dot(b), 1))) / rho
    return x 

def lu_factor(A, rho):
    r"""
    Compute LU factorisation of either :math:`A^T A + \rho I` or
    :math:`A A^T + \rho I`, depending on which matrix is smaller.

    Parameters
    ----------
    A : array_like
      Array :math:`A`
    rho : float
      Scalar :math:`\rho`

    Returns
    -------
    lu : ndarray
      Matrix containing U in its upper triangle, and L in its lower triangle,
      as returned by :func:`scipy.linalg.lu_factor`
    piv : ndarray
      Pivot indices representing the permutation matrix P, as returned by
      :func:`scipy.linalg.lu_factor`
    """

    N, M = A.shape
    # If N < M it is cheaper to factorise A*A^T + rho*I and then use the
    # matrix inversion lemma to compute the inverse of A^T*A + rho*I
    if N >= M:
        lu, piv = linalg.lu_factor(A.T.dot(A) + rho*np.identity(M,
                                   dtype=A.dtype))
    else:
        lu, piv = linalg.lu_factor(A.dot(A.T) + rho*np.identity(N,
                                   dtype=A.dtype))
    return lu, piv 

def __init__(self):
        matrix_pairs = {}
        rng = np.random.RandomState(SEED)
        for sz in SIZES:
            for transposed in [True, False]:
                inverse_array = rng.standard_normal((sz, sz))
                inverse_lu_and_piv = sla.lu_factor(
                    inverse_array.T if transposed else inverse_array)
                array = nla.inv(inverse_array)
                matrix_pairs[(sz, transposed)] = (
                    matrices.InverseLUFactoredSquareMatrix(
                        inverse_array, inverse_lu_and_piv, transposed), array)
            super().__init__(matrix_pairs, rng) 

def lu_and_piv(self):
        """Pivoted LU factorisation of matrix."""
        if self._lu_and_piv is None:
            self._lu_and_piv = sla.lu_factor(self._array, check_finite=False)
            self._lu_transposed = False
        return self._lu_and_piv 

def __init__(self, array, lu_and_piv=None, lu_transposed=None):
        """
        Args:
            array (array): 2D array specifying matrix entries.
            lu_and_piv (Tuple[array, array]): Pivoted LU factorisation
                represented as a tuple with first element a 2D array containing
                the lower and upper triangular factors (with the unit diagonal
                of the lower triangular factor not stored) and the second
                element a 1D array containing the pivot indices. Corresponds
                to the output of `scipy.linalg.lu_factor` and input to
                `scipy.linalg.lu_solve`.
            lu_transposed (bool): Whether LU factorisation is of original array
                or its transpose.
        """
        super().__init__(array.shape, _array=array)
        self._lu_and_piv = lu_and_piv
        self._lu_transposed = lu_transposed 

def factored(self):
        "Caches the LU factorisation of the matrix"
        try:
            return self.lu_factored
        except AttributeError:
            self.lu_factored = la.lu_factor(self.val().simple_view())
            return self.lu_factored 

def factorise_interp_matrix(self):
        if not self.factorisation_current:
            A, self.left_scaling, self.right_scaling = self.interpolation_matrix()
            self.lu, self.piv = LA.lu_factor(A)
            self.factorisation_current = True
        return 

def statdist(generator):
    """Compute the stationary distribution of a Markov chain.

    Parameters
    ----------
    generator : array_like
        Infinitesimal generator matrix of the Markov chain.

    Returns
    -------
    dist : numpy.ndarray
        The unnormalized stationary distribution of the Markov chain.

    Raises
    ------
    ValueError
        If the Markov chain does not have a unique stationary distribution.
    """
    generator = np.asarray(generator)
    n = generator.shape[0]
    with warnings.catch_warnings():
        # The LU decomposition raises a warning when the generator matrix is
        # singular (which it, by construction, is!).
        warnings.filterwarnings("ignore")
        lu, piv = spl.lu_factor(generator.T, check_finite=False)
    # The last row contains 0's only.
    left = lu[:-1,:-1]
    right = -lu[:-1,-1]
    # Solves system `left * x = right`. Assumes that `left` is
    # upper-triangular (ignores lower triangle).
    try:
        res = spl.solve_triangular(left, right, check_finite=False)
    except:
        # Ideally we would like to catch `spl.LinAlgError` only, but there seems
        # to be a bug in scipy, in the code that raises the LinAlgError (!!).
        raise ValueError(
                "stationary distribution could not be computed. "
                "Perhaps the Markov chain has more than one absorbing class?")
    res = np.append(res, 1.0)
    return (n / res.sum()) * res 

def __init__(self, M):
        self.M_lu = lu_factor(M)
        self.shape = M.shape
        self.dtype = M.dtype 

def linearize(self, inputs, outputs, partials):
        system_size = self.system_size
        self.lu = lu_factor(inputs['mtx'])

        partials['circulations', 'circulations'] = inputs['mtx'].flatten()
        partials['circulations', 'mtx'] = \
            np.outer(np.ones(system_size), outputs['circulations']).flatten() 

def solve_nonlinear(self, inputs, outputs):
        self.lu = lu_factor(inputs['mtx'])

        outputs['circulations'] = lu_solve(self.lu, inputs['rhs']) 

def __init__(self, M):
        self.M_lu = lu_factor(M)
        LinearOperator.__init__(self, M.shape, self._matvec, dtype=M.dtype) 

def __init__(self, M):
        self.M_lu = lu_factor(M)
        self.shape = M.shape
        self.dtype = M.dtype 

def compute_b(G_lst, GtG_lst, beta_lst, Rc0, num_bands, a_ri, use_lu=False, GtG_inv_lst=None):
    """
    compute the uniform sinusoidal samples b from the updated annihilating
    filter coeffiients.

    Parameters
    ----------
    GtG_lst: 
        list of G^H G for different subbands
    beta_lst: 
        list of beta-s for different subbands
    Rc0: 
        right-dual matrix, here it is the convolution matrix associated with c
    num_bands: 
        number of bands
    a_ri: 
        a 2D numpy array. each column corresponds to the measurements within a subband
    """
    b_lst = []
    a_Gb_lst = []
    if use_lu:
        assert GtG_inv_lst is not None
        lu_lst = []
        for loop in range(num_bands):
            GtG_loop = GtG_lst[loop]
            beta_loop = beta_lst[loop]
            lu_loop = linalg.lu_factor(
                np.dot(Rc0, np.dot(GtG_inv_lst[loop], Rc0.T)),
                check_finite=False
            )
            b_loop = beta_loop - \
                     np.dot(GtG_inv_lst[loop],
                                     np.dot(Rc0.T,
                                            linalg.lu_solve(lu_loop, np.dot(Rc0, beta_loop),
                                                            check_finite=False)),
                                     )

            b_lst.append(b_loop)
            a_Gb_lst.append(a_ri[:, loop] - np.dot(G_lst[loop], b_loop))
            lu_lst.append(lu_loop)

        return np.column_stack(b_lst), linalg.norm(np.concatenate(a_Gb_lst)), lu_lst
    else:
        for loop in range(num_bands):
            GtG_loop = GtG_lst[loop]
            beta_loop = beta_lst[loop]
            b_loop = beta_loop - \
                     linalg.solve(GtG_loop,
                                  np.dot(Rc0.T,
                                         linalg.solve(np.dot(Rc0, linalg.solve(GtG_loop, Rc0.T)),
                                                      np.dot(Rc0, beta_loop)))
                                  )

            b_lst.append(b_loop)
            a_Gb_lst.append(a_ri[:, loop] - np.dot(G_lst[loop], b_loop))

        return np.column_stack(b_lst), linalg.norm(np.concatenate(a_Gb_lst)) 

def solve_constraint_lagrangian(x, jac_x, c_x, c_prime):
    """
    Function which solves the problem
    minimize || J.dot(x_mod - x) || 
    subject to C(x_mod) = 0 
    via the method of Lagrange multipliers.

    Parameters
    ----------
    x : 1D numpy array
        Parameter values at x
    jac_x : 2D numpy array.
        The (estimated, approximate or exact)
        value of the Jacobian J(x)
    c_x : 1D numpy array
        Values of the constraints at x
    c_prime : 2D array of floats
        The Jacobian of the constraints
        (A, where A.x + b = 0)

    Returns
    -------
    x_mod : 1D numpy array
        The parameter values which minimizes the L2-norm
        of any function which has the Jacobian jac_x.
    lagrange_multipliers : 1D numpy array
        The multipliers for each of the equality 
        constraints
    """
    n_x = len(x)
    n = n_x + len(c_x)
    A = np.zeros((n, n))
    b = np.zeros(n)

    JTJ = jac_x.T.dot(jac_x)
    A[:n_x,:n_x] = JTJ/np.linalg.norm(JTJ)*n*n # includes scaling
    A[:n_x,n_x:] = c_prime.T
    A[n_x:,:n_x] = c_prime
    b[n_x:] = c_x


    luA = lu_factor(A) 
    dx_m = lu_solve(luA, -b) # lu_solve computes the solution of ax = b
    
    x_mod = x + dx_m[:n_x]
    lagrange_multipliers = dx_m[n_x:]
    return (x_mod, lagrange_multipliers) 

def __init__(self, fun, t0, y0, t_bound, max_step=np.inf,
                 rtol=1e-3, atol=1e-6, jac=None, jac_sparsity=None,
                 vectorized=False, **extraneous):
        warn_extraneous(extraneous)
        super(Radau, self).__init__(fun, t0, y0, t_bound, vectorized)
        self.y_old = None
        self.max_step = validate_max_step(max_step)
        self.rtol, self.atol = validate_tol(rtol, atol, self.n)
        self.f = self.fun(self.t, self.y)
        # Select initial step assuming the same order which is used to control
        # the error.
        self.h_abs = select_initial_step(
            self.fun, self.t, self.y, self.f, self.direction,
            3, self.rtol, self.atol)
        self.h_abs_old = None
        self.error_norm_old = None

        self.newton_tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5))
        self.sol = None

        self.jac_factor = None
        self.jac, self.J = self._validate_jac(jac, jac_sparsity)
        if issparse(self.J):
            def lu(A):
                self.nlu += 1
                return splu(A)

            def solve_lu(LU, b):
                return LU.solve(b)

            I = eye(self.n, format='csc')
        else:
            def lu(A):
                self.nlu += 1
                return lu_factor(A, overwrite_a=True)

            def solve_lu(LU, b):
                return lu_solve(LU, b, overwrite_b=True)

            I = np.identity(self.n)

        self.lu = lu
        self.solve_lu = solve_lu
        self.I = I

        self.current_jac = True
        self.LU_real = None
        self.LU_complex = None
        self.Z = None 

def __init__(self, fun, t0, y0, t_bound, max_step=np.inf,
                 rtol=1e-3, atol=1e-6, jac=None, jac_sparsity=None,
                 vectorized=False, **extraneous):
        warn_extraneous(extraneous)
        super(BDF, self).__init__(fun, t0, y0, t_bound, vectorized,
                                  support_complex=True)
        self.max_step = validate_max_step(max_step)
        self.rtol, self.atol = validate_tol(rtol, atol, self.n)
        f = self.fun(self.t, self.y)
        self.h_abs = select_initial_step(self.fun, self.t, self.y, f,
                                         self.direction, 1,
                                         self.rtol, self.atol)
        self.h_abs_old = None
        self.error_norm_old = None

        self.newton_tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5))

        self.jac_factor = None
        self.jac, self.J = self._validate_jac(jac, jac_sparsity)
        if issparse(self.J):
            def lu(A):
                self.nlu += 1
                return splu(A)

            def solve_lu(LU, b):
                return LU.solve(b)

            I = eye(self.n, format='csc', dtype=self.y.dtype)
        else:
            def lu(A):
                self.nlu += 1
                return lu_factor(A, overwrite_a=True)

            def solve_lu(LU, b):
                return lu_solve(LU, b, overwrite_b=True)

            I = np.identity(self.n, dtype=self.y.dtype)

        self.lu = lu
        self.solve_lu = solve_lu
        self.I = I

        kappa = np.array([0, -0.1850, -1/9, -0.0823, -0.0415, 0])
        self.gamma = np.hstack((0, np.cumsum(1 / np.arange(1, MAX_ORDER + 1))))
        self.alpha = (1 - kappa) * self.gamma
        self.error_const = kappa * self.gamma + 1 / np.arange(1, MAX_ORDER + 2)

        D = np.empty((MAX_ORDER + 3, self.n), dtype=self.y.dtype)
        D[0] = self.y
        D[1] = f * self.h_abs * self.direction
        self.D = D

        self.order = 1
        self.n_equal_steps = 0
        self.LU = None 

def sparse_factor_solve_kkt_reg(spH_, A_, spA_, spC_tilde, rx, rs, rz, ry, neq, nineq, nz):
    if neq > 0:
        g_ = torch.cat([rx, rs], 1)
        h_ = torch.cat([rz, ry], 1)
    else:
        g_ = torch.cat([rx, rs], 1)
        h_ = rz
    # XXX Not batched for now
    g_ = g_.squeeze(0).numpy()
    h_ = h_.squeeze(0).numpy()

    H_LU = splu(spH_)

    invH_A_ = H_LU.solve(A_.transpose())
    invH_g_ = H_LU.solve(g_)

    S_ = spA_.dot(invH_A_)  # A H-1 AT
    # A H-1 AT + C_tilde
    S_ -= spC_tilde
    # S_LU = btrifact_hack(S_)
    S_LU = lu_factor(S_)
    # [(H-1 g)T AT]T - h = A H-1 g - h
    t_ = spA_.dot(invH_g_) - h_
    # w = (A H-1 AT + C_tilde)-1 (A H-1 g - h) <= Av - eps I w = h
    # w_ = -t_.btrisolve(*S_LU)
    w_ = -lu_solve(S_LU, t_)
    # XXX Shouldn't it be just g (no minus)?
    # (Doesn't seem to make a difference, though...)
    t_ = -g_ - spA_.transpose().dot(w_)  # -g - AT w
    # v_ = t_.btrisolve(*H_LU)  # v = H-1 (-g - AT w)
    v_ = H_LU.solve(t_)

    dx = v_[:nz]
    ds = v_[nz:]
    dz = w_[:nineq]
    dy = w_[nineq:] if neq > 0 else None

    dx = torch.DoubleTensor(dx).unsqueeze(0)
    ds = torch.DoubleTensor(ds).unsqueeze(0)
    dz = torch.DoubleTensor(dz).unsqueeze(0)
    dy = torch.DoubleTensor(dy).unsqueeze(0) if neq > 0 else None

    return dx, ds, dz, dy 

def get_all_vectors(self):
        """Function to compute all vectors shorter than cutoff distance.
        """

        # Create a matrix of basis vectors.
        basis = np.array(self.input_vectors, dtype=float).T

        # Create ability to invert it.
        det_basis = np.linalg.det(basis)


        if det_basis == 0 or det_basis < 1e-14:
            raise RuntimeError("Vectors are not linearly independent.")

        fac = lu_factor(basis)

        # Compute range of each variable.
        cutoff_distance = np.math.sqrt(self.cutoff_distance_sq)
        step_range = []

        for i in range(3):
            max_disp = 0.0
            for j in range(3):
                max_disp += np.dot(self.input_vectors[i], self.input_vectors[
                    j]) / norm(self.input_vectors[i])
            step_range.append(int(np.math.ceil(max_disp / cutoff_distance)) + 1)

        # Ensure that we have sufficient range to get the cutoff distance
        # away from the origin by checking that we have large enough range to
        #  access a point cutoff distance away along the direction of xy,
        # xz and yz cross products.
        for dir in range(3):
            point = np.cross(self.input_vectors[dir], self.input_vectors[(dir
                                    + 1) % 3])
            point = point * cutoff_distance / norm(point)
            sln = lu_solve(fac, point)
            step_range = [max(step_range[i], int(np.math.ceil(abs(sln[i]))))
                          for i in range(3)]

        # Create the initial vector.
        for x in range(-step_range[0], 1 + step_range[0]):
            for y in range(-step_range[1], 1 + step_range[1]):
                for z in range(-step_range[2], 1 + step_range[2]):
                    a = np.array([x, y, z])
                    l = self.compute_vector(a)
                    dist_sq = l[0] ** 2 + l[1] ** 2 +  l[2] ** 2
                    if dist_sq <= self.cutoff_distance_sq:
                        if not self.include_zero and x == 0 and y == 0 and z \
                                == 0:
                            continue
                        self.super_cells.append(a)
                        self.vectors.append(l) 

def __init__(self, fun, t0, y0, t_bound, max_step=np.inf,
                 rtol=1e-3, atol=1e-6, jac=None, jac_sparsity=None,
                 vectorized=False, **extraneous):
        warn_extraneous(extraneous)
        super(Radau, self).__init__(fun, t0, y0, t_bound, vectorized)
        self.y_old = None
        self.max_step = validate_max_step(max_step)
        self.rtol, self.atol = validate_tol(rtol, atol, self.n)
        self.f = self.fun(self.t, self.y)
        # Select initial step assuming the same order which is used to control
        # the error.
        self.h_abs = select_initial_step(
            self.fun, self.t, self.y, self.f, self.direction,
            3, self.rtol, self.atol)
        self.h_abs_old = None
        self.error_norm_old = None

        self.newton_tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5))
        self.sol = None

        self.jac_factor = None
        self.jac, self.J = self._validate_jac(jac, jac_sparsity)
        if issparse(self.J):
            def lu(A):
                self.nlu += 1
                return splu(A)

            def solve_lu(LU, b):
                return LU.solve(b)

            I = eye(self.n, format='csc')
        else:
            def lu(A):
                self.nlu += 1
                return lu_factor(A, overwrite_a=True)

            def solve_lu(LU, b):
                return lu_solve(LU, b, overwrite_b=True)

            I = np.identity(self.n)

        self.lu = lu
        self.solve_lu = solve_lu
        self.I = I

        self.current_jac = True
        self.LU_real = None
        self.LU_complex = None
        self.Z = None 

def __init__(self, fun, t0, y0, t_bound, max_step=np.inf,
                 rtol=1e-3, atol=1e-6, jac=None, jac_sparsity=None,
                 vectorized=False, **extraneous):
        warn_extraneous(extraneous)
        super(BDF, self).__init__(fun, t0, y0, t_bound, vectorized,
                                  support_complex=True)
        self.max_step = validate_max_step(max_step)
        self.rtol, self.atol = validate_tol(rtol, atol, self.n)
        f = self.fun(self.t, self.y)
        self.h_abs = select_initial_step(self.fun, self.t, self.y, f,
                                         self.direction, 1,
                                         self.rtol, self.atol)
        self.h_abs_old = None
        self.error_norm_old = None

        self.newton_tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5))

        self.jac_factor = None
        self.jac, self.J = self._validate_jac(jac, jac_sparsity)
        if issparse(self.J):
            def lu(A):
                self.nlu += 1
                return splu(A)

            def solve_lu(LU, b):
                return LU.solve(b)

            I = eye(self.n, format='csc', dtype=self.y.dtype)
        else:
            def lu(A):
                self.nlu += 1
                return lu_factor(A, overwrite_a=True)

            def solve_lu(LU, b):
                return lu_solve(LU, b, overwrite_b=True)

            I = np.identity(self.n, dtype=self.y.dtype)

        self.lu = lu
        self.solve_lu = solve_lu
        self.I = I

        kappa = np.array([0, -0.1850, -1/9, -0.0823, -0.0415, 0])
        self.gamma = np.hstack((0, np.cumsum(1 / np.arange(1, MAX_ORDER + 1))))
        self.alpha = (1 - kappa) * self.gamma
        self.error_const = kappa * self.gamma + 1 / np.arange(1, MAX_ORDER + 2)

        D = np.empty((MAX_ORDER + 3, self.n), dtype=self.y.dtype)
        D[0] = self.y
        D[1] = f * self.h_abs * self.direction
        self.D = D

        self.order = 1
        self.n_equal_steps = 0
        self.LU = None