Python scipy.sparse.linalg.spsolve() Examples

def test_set_up_tms(self, tms_sphere):
        m, cond, dAdt, E_analytical = tms_sphere
        S = fem.FEMSystem.tms(m, cond)
        b = S.assemble_tms_rhs(dAdt)
        x = spalg.spsolve(S.A, b)
        v = mesh_io.NodeData(x, 'FEM', mesh=m)
        E = -v.gradient().value * 1e3 - dAdt.node_data2elm_data().value
        m.elmdata = [mesh_io.ElementData(E_analytical, 'analytical'),
                     mesh_io.ElementData(E, 'E_FEM'),
                     mesh_io.ElementData(E_analytical + dAdt.node_data2elm_data().value, 'grad_analytical'),
                     dAdt]

        m.nodedata = [mesh_io.NodeData(x, 'FEM')]
        #mesh_io.write_msh(m, '~/Tests/fem.msh')
        assert rdm(E, E_analytical) < .2
        assert np.abs(mag(E, E_analytical)) < np.log(1.1) 

def solve_direct_lagrange(k_lg, f):
    """Solves a linear system of equations (Ku = f) using the direct solver method and the
    Lagrangian multiplier method.

    :param k: (N+1) x (N+1) Lagrangian multiplier matrix of the linear system
    :type k: :class:`scipy.sparse.csc_matrix`
    :param f: N x 1 right hand side of the linear system
    :type f: :class:`numpy.ndarray`

    :return: The solution vector to the linear system of equations
    :rtype: :class:`numpy.ndarray`

    :raises RuntimeError: If the Lagrangian multiplier method exceeds a tolerance of 1e-5
    """

    u = spsolve(k_lg, np.append(f, 0))

    # compute error
    err = u[-1] / max(np.absolute(u))

    if err > 1e-5:
        err = "Lagrangian multiplier method error exceeds tolerance of 1e-5."
        raise RuntimeError(err)

    return u[:-1] 

def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)
        me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A, rhs.values)
        return me 

def solve_system_1(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(u0)
        me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A, rhs.values)
        return me 

def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(u0)
        me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A, rhs.values)
        return me 

def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)
        self.uext.values[0] = 0.0
        self.uext.values[-1] = 0.0
        self.uext.values[1:-1] = rhs.values[:]
        me.values = spsolve(sp.eye(self.params.nvars + 2, format='csc') - factor * self.A, self.uext.values)[1:-1]
        # me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A[1:-1, 1:-1], rhs.values)
        return me 

def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the nonlinear system
            factor (float): abbrev. for the node-to-node stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver (not used here so far)
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        M = self.Id - factor * self.A

        b = np.concatenate((rhs.values[0, :], rhs.values[1, :]))

        sol = spsolve(M, b)

        me = self.dtype_u(self.init)
        me.values[0, :], me.values[1, :] = np.split(sol, 2)

        return me 

def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)

        if self.params.direct_solver:
            me.values = spsolve(self.Id - factor * self.A, rhs.values.flatten())
        else:
            me.values = gmres(self.Id - factor * self.A, rhs.values.flatten(), x0=u0.values.flatten(),
                              tol=self.params.lintol, maxiter=self.params.liniter)[0]
        me.values = me.values.reshape(self.params.nvars)
        return me 

def timestep(self, u0, dt):

        # Solve for stages
        for i in range(0, self.nstages):

            # Construct RHS
            rhs = np.copy(u0)
            for j in range(0, i):
                rhs += dt * self.A_hat[i, j] * (self.f_slow(self.stages[j, :])) + dt * self.A[i, j] * \
                    (self.f_fast(self.stages[j, :]))

            # Solve for stage i
            if self.A[i, i] == 0:
                # Avoid call to spsolve with identity matrix
                self.stages[i, :] = np.copy(rhs)
            else:
                self.stages[i, :] = self.f_fast_solve(rhs, dt * self.A[i, i])

        # Update
        for i in range(0, self.nstages):
            u0 += dt * self.b_hat[i] * (self.f_slow(self.stages[i, :])) + dt * self.b[i] * \
                (self.f_fast(self.stages[i, :]))

        return u0 

def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)

        if self.params.direct_solver:
            me.values = spsolve(self.Id - factor * self.A, rhs.values.flatten())
        else:
            me.values = gmres(self.Id - factor * self.A, rhs.values.flatten(), x0=u0.values.flatten(),
                              tol=self.params.lintol, maxiter=self.params.liniter)[0]
        me.values = me.values.reshape(self.params.nvars)
        return me 

def solve_system_jacobian(self, dfdu, rhs, factor, u0, t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            dfdu: the Jacobian of the RHS of the ODE
            rhs: right-hand side for the linear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        me = self.dtype_u(self.init)
        me.values = spsolve(sp.eye(self.params.nvars) - factor * dfdu, rhs.values)
        return me 

def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        me = mesh(self.nvars)
        me.values = LA.spsolve(sp.eye(self.nvars)-factor*self.A,rhs.values)
        return me 

def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """
        
        M  = self.Id - factor*self.A
        
        b = np.concatenate( (rhs.values[0,:], rhs.values[1,:]) )
        
        sol = LA.spsolve(M, b)

        me = mesh(self.nvars)
        me.values[0,:], me.values[1,:] = np.split(sol, 2)
        
        return me 

def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """
        
        M  = self.Id - factor*self.A
        
        b = np.concatenate( (rhs.values[0,:], rhs.values[1,:]) )
        
        sol = LA.spsolve(M, b)

        me = mesh(self.nvars)
        me.values[0,:], me.values[1,:] = np.split(sol, 2)
        
        return me 

def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        me = mesh(self.nvars)
        me.values = LA.spsolve(sp.eye(self.nvars)-factor*self.A,rhs.values)

        return me 

def test_dirichlet_problem_cube(self, cube_lr):
        m = cube_lr.crop_mesh(5)
        cond = np.ones(len(m.elm.tetrahedra))
        top = m.nodes.node_number[m.nodes.node_coord[:, 2] > 49]
        bottom = m.nodes.node_number[m.nodes.node_coord[:, 2] < -49]
        bcs = [fem.DirichletBC(top, np.ones(len(top))),
               fem.DirichletBC(bottom, -np.ones(len(bottom)))]
        S = fem.FEMSystem(m, cond)
        A = S.A
        b = np.zeros(m.nodes.nr)
        dof_map = S.dof_map
        for bc in bcs:
            A, b, dof_map = bc.apply(A, b, dof_map)
        x = spalg.spsolve(A, b)
        for bc in bcs:
            x, dof_map = bc.apply_to_solution(x, dof_map)
        order = dof_map.inverse.argsort()
        x = x[order]
        sol = m.nodes.node_coord[:, 2]/50.
        assert np.allclose(sol, x.T) 

def solve(self):

        shape = self.bcs.shape
        if self._L is None:
            self._L = self.lhs.matrix(shape) # expensive operation, so cache it

        L = sparse.lil_matrix(self._L)
        f = self.rhs.reshape(-1, 1)

        nz = list(self.bcs.row_inds())

        L[nz, :] = self.bcs.lhs[nz, :]
        f[nz] = np.array(self.bcs.rhs[nz].toarray()).reshape(-1, 1)

        L = sparse.csr_matrix(L)
        return spsolve(L, f).reshape(shape) 

def solve ( self, A, b ):

    if len(A.shape) == 2:
      C = self.getConstraintsMatrix()

      a = zeros(len(self))
      a[self.constrainedDofs] = self.constrainedFac * array(self.constrainedVals)

      A_constrained = C.transpose() * (A * C )
      b_constrained = C.transpose()* ( b - A * a )

      x_constrained = spsolve( A_constrained, b_constrained )

      x = C * x_constrained

      x[self.constrainedDofs] = self.constrainedFac * array(self.constrainedVals)
    
    elif len(A.shape) == 1:
      x = b / A

      x[self.constrainedDofs] = self.constrainedFac * array(self.constrainedVals)
   
    return x 

def static_sol(mat, rhs):
    """Solve a static problem [mat]{u_sol} = {rhs}
    
    Parameters
    ----------
    mat : array
        Array with the system of equations. It can be stored in
        dense or sparse scheme.
    rhs : array
        Array with right-hand-side of the system of equations.
    
    Returns
    -------
    u_sol : array
        Solution of the system of equations.

    Raises
    ------
    
        

    """
    if type(mat) is csr_matrix:
        u_sol = spsolve(mat, rhs)
    elif type(mat) is ndarray:
        u_sol = solve(mat, rhs)
    else:
        raise TypeError("Matrix should be numpy array or csr_matrix.")

    return u_sol 

def calculate_PTDF(sub_network,skip_pre=False):
    """
    Calculate the Power Transfer Distribution Factor (PTDF) for
    sub_network.

    Sets sub_network.PTDF as a (dense) numpy array.

    Parameters
    ----------
    sub_network : pypsa.SubNetwork
    skip_pre : bool, default False
        Skip the preliminary steps of computing topology, calculating dependent values,
        finding bus controls and computing B and H.

    """

    if not skip_pre:
        calculate_B_H(sub_network)

    #calculate inverse of B with slack removed

    n_pvpq = len(sub_network.pvpqs)
    index = np.r_[:n_pvpq]

    I = csc_matrix((np.ones(n_pvpq), (index, index)))

    B_inverse = spsolve(sub_network.B[1:, 1:],I)

    #exception for two-node networks, where B_inverse is a 1d array
    if issparse(B_inverse):
        B_inverse = B_inverse.toarray()
    elif B_inverse.shape == (1,):
        B_inverse = B_inverse.reshape((1,1))

    #add back in zeroes for slack
    B_inverse = np.hstack((np.zeros((n_pvpq,1)),B_inverse))
    B_inverse = np.vstack((np.zeros(n_pvpq+1),B_inverse))

    sub_network.PTDF = sub_network.H*B_inverse 

def solve_one_iter(self):
        """Solve one iteration of Gauss-Newton."""
        # precision * dx = information
        precision, information, cost = self._get_precision_information_and_cost()
        dx = splinalg.spsolve(precision, information)

        # Backtrack line search
        if self.options.linesearch_max_iters > 0:
            best_step_size, best_cost = self._do_line_search(dx)
        else:
            best_step_size, best_cost = 1., cost

        return best_step_size * dx, best_cost 

def inv(A):
    """
    Compute the inverse of a sparse matrix

    Parameters
    ----------
    A : (M,M) ndarray or sparse matrix
        square matrix to be inverted

    Returns
    -------
    Ainv : (M,M) ndarray or sparse matrix
        inverse of `A`

    Notes
    -----
    This computes the sparse inverse of `A`.  If the inverse of `A` is expected
    to be non-sparse, it will likely be faster to convert `A` to dense and use
    scipy.linalg.inv.

    .. versionadded:: 0.12.0

    """
    I = speye(A.shape[0], A.shape[1], dtype=A.dtype, format=A.format)
    Ainv = spsolve(A, I)
    return Ainv 

def _solve_P_Q(U, V, structure=None):
    """
    A helper function for expm_2009.

    Parameters
    ----------
    U : ndarray
        Pade numerator.
    V : ndarray
        Pade denominator.
    structure : str, optional
        A string describing the structure of both matrices `U` and `V`.
        Only `upper_triangular` is currently supported.

    Notes
    -----
    The `structure` argument is inspired by similar args
    for theano and cvxopt functions.

    """
    P = U + V
    Q = -U + V
    if isspmatrix(U):
        return spsolve(Q, P)
    elif structure is None:
        return solve(Q, P)
    elif structure == UPPER_TRIANGULAR:
        return solve_triangular(Q, P)
    else:
        raise ValueError('unsupported matrix structure: ' + str(structure)) 

def _solve_direct(A, b):
    """ Direct solver """

    if HAS_MKL:
        # prefered method using MKL. Much faster (on Mac at least)
        pSolve = pardisoSolver(A, mtype=13)
        pSolve.factor()
        x = pSolve.solve(b)
        pSolve.clear()
        return x
    else:
        # scipy solver.
        return spl.spsolve(A, b) 

def test_basic_spsolve_vector():
    ps.remove_stored_factorization()
    ps.free_memory()
    A, b = create_test_A_b_rand()
    xpp = spsolve(A,b)
    xscipy = scipyspsolve(A,b)
    np.testing.assert_array_almost_equal(xpp, xscipy) 

def test_basic_spsolve_matrix():
    ps.remove_stored_factorization()
    ps.free_memory()
    A, b = create_test_A_b_rand(matrix=True)
    xpp = spsolve(A,b)
    xscipy = scipyspsolve(A,b)
    np.testing.assert_array_almost_equal(xpp, xscipy) 

def solve_direct(k, f):
    """Solves a linear system of equations (Ku = f) using the direct solver method.

    :param k: N x N matrix of the linear system
    :type k: :class:`scipy.sparse.csc_matrix`
    :param f: N x 1 right hand side of the linear system
    :type f: :class:`numpy.ndarray`

    :return: The solution vector to the linear system of equations
    :rtype: :class:`numpy.ndarray`
    """

    return spsolve(k, f) 

def solve_linear(model):
    logger.info('solving problem with %d DOFs...'%model.DOF)
    K_,f_=model.K_,model.f_
#    M_x = lambda x: sl.spsolve(P, x)
#    M = sl.LinearOperator((n, n), M_x)
    #print(sl.spsolve(K_,f_))
    delta,info=sl.lgmres(K_,f_.toarray())
    model.is_solved=True
    logger.info('Done!')
    model.d_=delta.reshape((model.node_count*6,1))
    model.r_=model.K*model.d_ 

def newton_raphson_sparse(f, guess, dfdx, x_tol=1e-10, lim_iter=100, distribute_slack=False, slack_weights=None):
    """Solve f(x) = 0 with initial guess for x and dfdx(x). dfdx(x) should
    return a sparse Jacobian.  Terminate if error on norm of f(x) is <
    x_tol or there were more than lim_iter iterations.

    """

    slack_args = {"distribute_slack": distribute_slack,
                  "slack_weights": slack_weights}
    converged = False
    n_iter = 0
    F = f(guess, **slack_args)
    diff = norm(F,np.Inf)

    logger.debug("Error at iteration %d: %f", n_iter, diff)

    while diff > x_tol and n_iter < lim_iter:

        n_iter +=1

        guess = guess - spsolve(dfdx(guess, **slack_args),F)

        F = f(guess, **slack_args)
        diff = norm(F,np.Inf)

        logger.debug("Error at iteration %d: %f", n_iter, diff)

    if diff > x_tol:
        logger.warning("Warning, we didn't reach the required tolerance within %d iterations, error is at %f. See the section \"Troubleshooting\" in the documentation for tips to fix this. ", n_iter, diff)
    elif not np.isnan(diff):
        converged = True

    return guess, n_iter, diff, converged 

def regionfillLaplace(I, mask, maskPerimeter):
    height, width = I.shape
    rightSide = formRightSide(I, maskPerimeter)

    # Location of mask pixels
    maskIdx = np.where(mask)

    # Only keep values for pixels that are in the mask
    rightSide = rightSide[maskIdx]

    # Number the mask pixels in a grid matrix
    grid = -np.ones((height, width))
    grid[maskIdx] = range(0, maskIdx[0].size)
    # Pad with zeros to avoid "index out of bounds" errors in the for loop
    grid = padMatrix(grid)
    gridIdx = np.where(grid >= 0)

    # Form the connectivity matrix D=sparse(i,j,s)
    # Connect each mask pixel to itself
    i = np.arange(0, maskIdx[0].size)
    j = np.arange(0, maskIdx[0].size)
    # The coefficient is the number of neighbors over which we average
    numNeighbors = computeNumberOfNeighbors(height, width)
    s = numNeighbors[maskIdx]
    # Now connect the N,E,S,W neighbors if they exist
    for direction in ((-1, 0), (0, 1), (1, 0), (0, -1)):
        # Possible neighbors in the current direction
        neighbors = grid[gridIdx[0] + direction[0], gridIdx[1] + direction[1]]
        # ConDnect mask points to neighbors with -1's
        index = (neighbors >= 0)
        i = np.concatenate((i, grid[gridIdx[0][index], gridIdx[1][index]]))
        j = np.concatenate((j, neighbors[index]))
        s = np.concatenate((s, -np.ones(np.count_nonzero(index))))

    D = sparse.coo_matrix((s, (i.astype(int), j.astype(int)))).tocsr()
    sol = spsolve(D, rightSide)
    I[maskIdx] = sol
    return I