Python scipy.linalg.lu_solve() Examples

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 solve(self, vec):
        """Solve the impedance matrix for a source vector. Caches the
        factorised matrix for efficiently solving multiple vectors"""
        if self.part_o != self.part_s:
            raise ValueError("Can only invert a self-impedance matrix")

        Z_lu = self.factored()
        if isinstance(vec, LookupArray):
            vec = vec.simple_view()

        lookup = (self.unknowns, (self.part_s, self.basis_container))

        if len(vec.shape) > 1:
            lookup = lookup+(vec.shape[1],)

        I = LookupArray(lookup, dtype=np.complex128)
        I_simp = I.simple_view()
        I_simp[:] = la.lu_solve(Z_lu, vec)
        return I 

def lu_compute_mtx_obj(Tbeta_lst, num_bands, K, lu_R_GtGinv_Rt_lst):
    """
    compute the matrix (M) in the objective function:
        min   c^H M c
        s.t.  c0^H c = 1

    Parameters
    ----------
    GtG_lst: 
        list of G^H * G
    Tbeta_lst: 
        list of Teoplitz matrices for beta-s
    Rc0: 
        right dual matrix for the annihilating filter (same of each block -> not a list)
    """
    mtx = np.zeros((K + 1, K + 1), dtype=float)  # <= assume G, Tbeta and Rc0 are real-valued

    for loop in range(num_bands):
        Tbeta_loop = Tbeta_lst[loop]
        mtx += np.dot(Tbeta_loop.T,
                      linalg.lu_solve(lu_R_GtGinv_Rt_lst[loop],
                                      Tbeta_loop, check_finite=False)
                      )

    return mtx 

def compute_obj_val(GtG_inv_lst, Tbeta_lst, Rc0, c_ri_half, num_bands, K):
    """
    compute the fitting error.
    CAUTION: Here we assume use_lu = True
    """
    mtx = np.zeros((K + 1, K + 1), dtype=float)  # <= assume G, Tbeta and Rc0 are real-valued
    fitting_error = 0
    for band_count in range(num_bands):
        #GtG_loop = GtG_lst[band_count]
        Tbeta_loop = Tbeta_lst[band_count]

        mtx += np.dot(Tbeta_loop.T,
                      linalg.solve(
                          np.dot(Rc0, 
                              #linalg.lu_solve(GtG_loop, Rc0.T, check_finite=False)
                              np.dot(GtG_inv_lst[band_count], Rc0.T)
                              ),
                          Tbeta_loop, check_finite=False
                      )
                      )

    fitting_error += np.sqrt(np.dot(c_ri_half.T, np.dot(mtx, c_ri_half)).real)

    return fitting_error, mtx 

def solve_storing_lu(A, b):
    LOG.debug(f"Solve with LU decomposition of {A}.")
    return sl.lu_solve(lu_decomp(A), b)


# ITERATIVE SOLVER 

def lu_solve(lu_A, b, trans=0):
    """
    LU solve wrapper.

    Computes the solution to

    .. math:: \mathbf{Ax} = \mathbf{b}

    or

    .. math:: \mathbf{A}^T\mathbf{x} = \mathbf{b}

    if ``trans=1``.

    It uses the ``SuperLU.solve()`` method if the input is a ``SuperLU`` or else will revert to the dense methods
    in scipy.

    Args:
        lu_A (SuperLU or tuple): object or tuple containing the information of the LU factorisation
        b (np.ndarray): Right hand side vector to solve
        trans (int): ``0`` or ``1`` for either solution option.

    Returns:
        np.ndarray: Solution to the system.

    """
    transpose_mode_dict = {0: 'N', 1: 'T'}
    if type(lu_A) == scsp.linalg.SuperLU:
        return lu_A.solve(b, trans=transpose_mode_dict[trans])
    else:
        return sclalg.lu_solve(lu_A, b, trans=trans) 

def mimo_block_arnoldi(self, frequency, r):

        n = self.ss.states

        A = self.ss.A
        B = self.ss.B
        C = self.ss.C

        for i in range(self.nfreq):

            if self.frequency[i] == np.inf:
                F = A
                G = B
            else:
                lu_a = krylovutils.lu_factor(frequency[i], A)
                F = krylovutils.lu_solve(lu_a, np.eye(n))
                G = krylovutils.lu_solve(lu_a, B)

            if i == 0:
                V = krylovutils.block_arnoldi_krylov(r, F, G)
            else:
                Vi = krylovutils.block_arnoldi_krylov(r, F, G)
                V = np.block([V, Vi])

        self.V = V

        Ar = V.T.dot(A.dot(V))
        Br = V.T.dot(B)
        Cr = C.dot(V)

        return Ar, Br, Cr 

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_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 _right_matrix_multiply(self, other):
        return sla.lu_solve(
            self._inv_lu_and_piv, other.T, not self._inv_lu_transposed,
            check_finite=False).T 

def _left_matrix_multiply(self, other):
        return sla.lu_solve(
            self._inv_lu_and_piv, other, self._inv_lu_transposed, 
            check_finite=False) 

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 _matvec(self, x):
        return lu_solve(self.M_lu, x) 

def solve_system(self, rhs):
        # To do preconditioning below, we will need to scale each column of A elementwise by the entries of some vector
        col_scale = lambda A, scale: (A.T * scale).T  # Uses the trick that A*x scales the 0th column of A by x[0], etc.

        if self.factorisation_current:
            # A(preconditioned) = diag(left_scaling) * A(original) * diag(right_scaling)
            # Solve A(original)\rhs
            return col_scale(LA.lu_solve((self.lu, self.piv), col_scale(rhs, self.left_scaling)), self.right_scaling)
        else:
            if self.do_logging:
                logging.warning("model.solve_system not using factorisation")
            A, left_scaling, right_scaling = self.interpolation_matrix()
            return col_scale(LA.solve(A, col_scale(rhs, left_scaling)), right_scaling) 

def _matvec(self, x):
        return lu_solve(self.M_lu, x) 

def solve_linear(self, d_outputs, d_residuals, mode):
        if mode == 'fwd':
            d_outputs['circulations'] = lu_solve(self.lu, d_residuals['circulations'], trans=0)
        else:
            d_residuals['circulations'] = lu_solve(self.lu, d_outputs['circulations'], trans=1) 

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

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

def _matvec(self, x):
        return lu_solve(self.M_lu, x) 

def _matvec(self, x):
        return lu_solve(self.M_lu, x) 

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 __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 __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 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(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 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)