Solve a real triangular system with multiple right-hand sides.
541691
def cublasStrsm(handle, side, uplo, trans, diag, m, n, alpha, A, lda, B, ldb):
    
    
    status = _libcublas.cublasStrsm_v2(handle, 
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       m, n, ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda, int(B), ldb)
    cublasCheckStatus(status)
199118
def solve_kkt(Q_LU, d, G, A, S_LU, rx, rs, rz, ry):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    invQ_rx = rx.btrisolve(*Q_LU)
    if neq > 0:
        h = torch.cat((invQ_rx.unsqueeze(1).bmm(A.transpose(1, 2)).squeeze(1) - ry,
                       invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz), 1)
    else:
        h = invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz

    w = -(h.btrisolve(*S_LU))

    g1 = -rx - w[:, neq:].unsqueeze(1).bmm(G).squeeze(1)
    if neq > 0:
        g1 -= w[:, :neq].unsqueeze(1).bmm(A).squeeze(1)
    g2 = -rs - w[:, neq:]

    dx = g1.btrisolve(*Q_LU)
    ds = g2 / d
    dz = w[:, neq:]
    dy = w[:, :neq] if neq > 0 else None

    return dx, ds, dz, dy
270610
def _brzozowski_algebraic_method_solve(self):
        """"""
        orderedstates = self._bfs_sort(
            sorted(
                self.mma.states,
                key=attrgetter('initial'),
                reverse=True)[0])
        for n in range(len(orderedstates) - 1, 0, -1):
            # print "n:" + repr(n)
            if self.A[
                    orderedstates[n].stateid,
                    orderedstates[n].stateid] != self.empty:
                # B[n] := star(A[n,n]) . B[n]
                if self.B[orderedstates[n].stateid] != self.empty:
                    self.B[orderedstates[n].stateid] = \
                        self.star(self.A[orderedstates[n].stateid, orderedstates[n].stateid]) \
                        + self.B[orderedstates[n].stateid]
                else:
                    self.B[orderedstates[n].stateid] = self.star(
                        self.A[orderedstates[n].stateid, orderedstates[n].stateid])
                for j in range(0, n):
                    # A[n,j] := star(A[n,n]) . A[n,j]
                    if self.A[
                            orderedstates[n].stateid,
                            orderedstates[j].stateid] != self.empty:
                        self.A[
                            orderedstates[n].stateid,
                            orderedstates[j].stateid] = \
                            self.star(self.A[orderedstates[n].stateid, orderedstates[n].stateid]) \
                            + self.A[orderedstates[n].stateid, orderedstates[j].stateid]
                    else:
                        self.A[orderedstates[n].stateid, orderedstates[j].stateid] = self.star(
                            self.A[orderedstates[n].stateid, orderedstates[n].stateid])
            for i in range(0, n):
                # B[i] += A[i,n] . B[n]
                newnode = None
                if self.A[orderedstates[i].stateid, orderedstates[n].stateid] != self.empty \
                        and self.B[orderedstates[n].stateid] != self.empty:
                    newnode = self.A[orderedstates[i].stateid, orderedstates[
                        n].stateid] + self.B[orderedstates[n].stateid]
                elif self.A[orderedstates[i].stateid, orderedstates[n].stateid] != self.empty:
                    newnode = self.A[
                        orderedstates[i].stateid,
                        orderedstates[n].stateid]
                elif self.B[orderedstates[n].stateid] != self.empty:
                    newnode = self.B[orderedstates[n].stateid]
                if self.B[orderedstates[i].stateid] != self.empty:
                    if newnode is not None:
                        self.B[orderedstates[i].stateid] += newnode
                else:
                    self.B[orderedstates[i].stateid] = newnode
                for j in range(0, n):
                    # A[i,j] += A[i,n] . A[n,j]
                    newnode = None
                    if self.A[
                            orderedstates[i].stateid,
                            orderedstates[n].stateid] != self.empty \
                            and self.A[orderedstates[n].stateid, orderedstates[j].stateid] \
                                    != self.empty:
                        newnode = self.A[orderedstates[i].stateid, orderedstates[
                            n].stateid] + self.A[orderedstates[n].stateid, orderedstates[j].stateid]
                    elif self.A[orderedstates[i].stateid, orderedstates[n].stateid] != self.empty:
                        newnode = self.A[
                            orderedstates[i].stateid,
                            orderedstates[n].stateid]
                    elif self.A[orderedstates[n].stateid, orderedstates[j].stateid] != self.empty:
                        newnode = self.A[
                            orderedstates[n].stateid,
                            orderedstates[j].stateid]
                    if self.A[
                            orderedstates[i].stateid,
                            orderedstates[j].stateid] != self.empty:
                        if newnode is not None:
                            self.A[
                                orderedstates[i].stateid,
                                orderedstates[j].stateid] += newnode
                    else:
                        self.A[
                            orderedstates[i].stateid,
                            orderedstates[j].stateid] = newnode
657230
def quadraticEval(a, b, c, x):
    """"""
    return a*(x**2) + b*x + c
300356
def cublasDsymm(handle, side, uplo, m, n, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasDsymm_v2(handle,
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo],
                                       m, n, ctypes.byref(ctypes.c_double(alpha)),
                                       int(A), lda, int(B), ldb, 
                                       ctypes.byref(ctypes.c_double(beta)), 
                                       int(C), ldc)
    cublasCheckStatus(status)
213036
public void solve(double[] b, double[] x) {
        if (!full) {
            throw new IllegalStateException("This is not a FULL singular value decomposition.");
        }

        if (b.length != m || x.length != n) {
            throw new IllegalArgumentException("Dimensions do not agree.");
        }

        double[] tmp = new double[n];
        for (int j = 0; j < n; j++) {
            double r = 0.0;
            if (s[j] > tol) {
                for (int i = 0; i < m; i++) {
                    r += U.get(i, j) * b[i];
                }
                r /= s[j];
            }
            tmp[j] = r;
        }

        for (int j = 0; j < n; j++) {
            double r = 0.0;
            for (int jj = 0; jj < n; jj++) {
                r += V.get(j, jj) * tmp[jj];
            }
            x[j] = r;
        }
    }
86176
public Vector4d hermite(Vector4dc t0, Vector4dc v1, Vector4dc t1, double t, Vector4d dest) {
        double t2 = t * t;
        double t3 = t2 * t;
        dest.x = (x + x - v1.x() - v1.x() + t1.x() + t0.x()) * t3 + (3.0 * v1.x() - 3.0 * x - t0.x() - t0.x() - t1.x()) * t2 + x * t + x;
        dest.y = (y + y - v1.y() - v1.y() + t1.y() + t0.y()) * t3 + (3.0 * v1.y() - 3.0 * y - t0.y() - t0.y() - t1.y()) * t2 + y * t + y;
        dest.z = (z + z - v1.z() - v1.z() + t1.z() + t0.z()) * t3 + (3.0 * v1.z() - 3.0 * z - t0.z() - t0.z() - t1.z()) * t2 + z * t + z;
        dest.w = (w + w - v1.w() - v1.w() + t1.w() + t0.w()) * t3 + (3.0 * v1.w() - 3.0 * w - t0.w() - t0.w() - t1.w()) * t2 + w * t + w;
        return dest;
    }
1047141
def xstep(self):
        
        """

        self.X = np.asarray(sl.cho_solve_ATAI(
            self.D, self.rho, self.DTS + self.rho * (self.Y - self.U),
            self.lu, self.piv), dtype=self.dtype)

        if self.opt['LinSolveCheck']:
            b = self.DTS + self.rho * (self.Y - self.U)
            ax = self.D.T.dot(self.D.dot(self.X)) + self.rho*self.X
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
212810
def cublasZtrsm(handle, side, uplo, transa, diag, m, n, alpha, A, lda, B, ldb):
    
    
    status = _libcublas.cublasZtrsm_v2(handle, 
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       m, n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,                    
                                                                               alpha.imag)),
                                       int(A), lda, int(B), ldb)
    cublasCheckStatus(status)
1289444
def cublasSsymm(handle, side, uplo, m, n, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasSsymm_v2(handle,
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       m, n, ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda, int(B), ldb, 
                                       ctypes.byref(ctypes.c_float(beta)), 
                                       int(C), ldc)
    cublasCheckStatus(status)
444449
def forward(inputs_i, Q, G, A, b, h, U_Q, U_S, R, verbose=False):
    
    nineq, nz, neq, _ = get_sizes(G, A)

    # find initial values
    d = torch.ones(nineq).type_as(Q)
    nb = -b if b is not None else None
    factor_kkt(U_S, R, d)
    x, s, z, y = solve_kkt(
        U_Q, d, G, A, U_S,
        inputs_i, torch.zeros(nineq).type_as(Q), -h, nb)
    # x1, s1, z1, y1 = factor_solve_kkt(Q, torch.eye(nineq).type_as(Q), G, A, inputs_i,
    # torch.zeros(nineq).type_as(Q), -h, nb)

    if torch.min(s) < 0:
        s -= torch.min(s) - 1
    if torch.min(z) < 0:
        z -= torch.min(z) - 1

    prev_resid = None
    for i in range(20):
        # affine scaling direction
        rx = (torch.mv(A.t(), y) if neq > 0 else 0.) + \
            torch.mv(G.t(), z) + torch.mv(Q, x) + inputs_i
        rs = z
        rz = torch.mv(G, x) + s - h
        ry = torch.mv(A, x) - b if neq > 0 else torch.Tensor([0.])
        mu = torch.dot(s, z) / nineq
        pri_resid = torch.norm(ry) + torch.norm(rz)
        dual_resid = torch.norm(rx)
        resid = pri_resid + dual_resid + nineq * mu
        d = z / s
        if verbose:
            print(("primal_res = {0:.5g}, dual_res = {1:.5g}, " +
                   "gap = {2:.5g}, kappa(d) = {3:.5g}").format(
                pri_resid, dual_resid, mu, min(d) / max(d)))
        # if (pri_resid < 5e-4 and dual_resid < 5e-4 and mu < 4e-4):
        improved = (prev_resid is None) or (resid < prev_resid + 1e-6)
        if not improved or resid < 1e-6:
            return x, y, z
        prev_resid = resid

        factor_kkt(U_S, R, d)
        dx_aff, ds_aff, dz_aff, dy_aff = solve_kkt(U_Q, d, G, A, U_S,
                                                   rx, rs, rz, ry)
        # D = torch.diag((z/s).cpu()).type_as(Q)
        # dx_aff1, ds_aff1, dz_aff1, dy_aff1 = factor_solve_kkt(Q, D, G, A, rx, rs, rz, ry)

        # compute centering directions
        alpha = min(min(get_step(z, dz_aff), get_step(s, ds_aff)), 1.0)
        sig = (torch.dot(s + alpha * ds_aff, z +
                         alpha * dz_aff) / (torch.dot(s, z)))**3
        dx_cor, ds_cor, dz_cor, dy_cor = solve_kkt(
            U_Q, d, G, A, U_S, torch.zeros(nz).type_as(Q),
            (-mu * sig * torch.ones(nineq).type_as(Q) + ds_aff * dz_aff) / s,
            torch.zeros(nineq).type_as(Q), torch.zeros(neq).type_as(Q))
        # dx_cor, ds_cor, dz_cor, dy_cor = factor_solve_kkt(Q, D, G, A,
        #     torch.zeros(nz).type_as(Q),
        #     (-mu*sig*torch.ones(nineq).type_as(Q) + ds_aff*dz_aff)/s,
        #     torch.zeros(nineq).type_as(Q), torch.zeros(neq).type_as(Q))

        dx = dx_aff + dx_cor
        ds = ds_aff + ds_cor
        dz = dz_aff + dz_cor
        dy = dy_aff + dy_cor if neq > 0 else None
        alpha = min(1.0, 0.999 * min(get_step(s, ds), get_step(z, dz)))
        dx_norm = torch.norm(dx)
        dz_norm = torch.norm(dz)
        if np.isnan(dx_norm) or dx_norm > 1e5 or dz_norm > 1e5:
            # Overflow, return early
            return x, y, z

        x += alpha * dx
        s += alpha * ds
        z += alpha * dz
        y = y + alpha * dy if neq > 0 else None

    return x, y, z
471662
def scipy_solve_symm_block_tridiag(H_diag, H_upper_diag, v, ab=None):
    
    from scipy.linalg import solveh_banded
    ab = convert_block_tridiag_to_banded(H_diag, H_upper_diag) \
        if ab is None else ab
    x = solveh_banded(ab, v.ravel(), lower=True)
    return x.reshape(v.shape)
24831
def solve(self):
        
        case = self.om.case
        self._base_mva = case.base_mva
        # TODO: Find an explanation for this value.
        self.opt["cost_mult"] = 1e-4

        # Unpack the OPF model.
        self._bs, self._ln, self._gn, _ = self._unpack_model(self.om)
        # Compute problem dimensions.
        self._ipol, _, self._nb, self._nl, _, self._ny, self._nxyz = \
            self._dimension_data(self._bs, self._ln, self._gn)

        # Compute problem dimensions.
        self._ng = len(self._gn)
#        gpol = [g for g in gn if g.pcost_model == POLYNOMIAL]

        # Linear constraints (l <= A*x <= u).
        A, l, u = self.om.linear_constraints()
#        AA, bb = self._linear_constraints(self.om)

        _, xmin, xmax = self._var_bounds()

        # Select an interior initial point for interior point solver.
        x0 = self._initial_interior_point(self._bs, self._gn, xmin, xmax, self._ny)

        # Build admittance matrices.
        self._Ybus, self._Yf, self._Yt = case.Y

        # Optimisation variables.

        self._Pg = self.om.get_var("Pg")
        self._Qg = self.om.get_var("Qg")
        self._Va = self.om.get_var("Va")
        self._Vm = self.om.get_var("Vm")

        # Adds a constraint on the reference bus angles.
#        xmin, xmax = self._ref_bus_angle_constraint(bs, Va, xmin, xmax)

        # Solve using Python Interior Point Solver (PIPS).
        s = self._solve(x0, A, l, u, xmin, xmax)

        Vang, Vmag, Pgen, Qgen = self._update_solution_data(s)

        self._update_case(self._bs, self._ln, self._gn, self._base_mva,
            self._Yf, self._Yt, Vang, Vmag, Pgen, Qgen, s["lmbda"])

        return s
892155
def rsolve(A, b, epsilon=_epsilon):
    
    """
    A = asarray(A, float)
    b = asarray(b, float)
    if A.shape[0] == 0:
        return zeros((A.shape[1],))
    if A.shape[1] == 0:
        return zeros((0,))
    try:
        x = lstsq(A, b, rcond=epsilon)
        r = sum(x[3] > epsilon)
        if r == 0:
            return zeros(A.shape[1])
        return x[0]
    except (ValueError, LinAlgError) as e:
        warnings.warn(str(e), RuntimeWarning)
    return solve(A, b)
368902
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        else if( B.numRows != numRows || B.numCols != X.numCols )
            throw new IllegalArgumentException("Unexpected dimensions for B");

        U = decomposer.getR(U,true);
        final double gammas[] = decomposer.getGammas();
        final double dataQR[] = QR.data;

        final int BnumCols = B.numCols;

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                a[i*2]   = B.data[indexB];
                a[i*2+1] = B.data[indexB+1];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^H)*b = b - u*(gamma*U^H*b)
            for( int n = 0; n < numCols; n++ ) {
                int indexU = (n*numRows + n + 1)*2;

                double realUb = a[n*2];
                double imagUb = a[n*2+1];

                // U^H*b
                for( int i = n+1; i < numRows; i++ ) {
                    double realU = dataQR[indexU++];
                    double imagU = -dataQR[indexU++];

                    double realB = a[i*2];
                    double imagB = a[i*2+1];

                    realUb += realU*realB - imagU*imagB;
                    imagUb += realU*imagB + imagU*realB;
                }

                // gamma*U^T*b
                realUb *= gammas[n];
                imagUb *= gammas[n];

                a[n*2]   -= realUb;
                a[n*2+1] -= imagUb;

                indexU = (n*numRows + n + 1)*2;
                for( int i = n+1; i < numRows; i++) {
                    double realU = dataQR[indexU++];
                    double imagU = dataQR[indexU++];

                    a[i*2]   -= realU*realUb - imagU*imagUb;
                    a[i*2+1] -= realU*imagUb + imagU*realUb;
                }
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_ZDRM.solveU(U.data, a, numCols);

            // save the results

            for( int i = 0; i < numCols; i++ ) {
                int indexX = (i*X.numCols+colB)*2;

                X.data[indexX]   = a[i*2];
                X.data[indexX+1] = a[i*2+1];
            }
        }
    }
1353711
def solve(self):
        
        result = Formula(path_actions(a_star_search(
            ({f: self.cube[f] for f in "LUFDRB"}, 
             self.cube.select_type("edge") & self.cube.has_colour(self.cube["D"].colour)), 
            self.cross_successors, 
            self.cross_state_value, 
            self.cross_goal, 
            )))
        self.cube(result)
        return result
856425
def solve_t(self, Mt):
        
        if len(Mt.shape)==2:    _Mt = Mt[:, :, sp.newaxis]
        else:                   _Mt = Mt
        LMt = vei_CoR_veX(_Mt, R=self.Lr(), C=self.Lc())
        DLMt = self.D()[:, :, sp.newaxis] * LMt
        RV = vei_CoR_veX(DLMt, R=self.Lr().T, C=self.Lc().T)
        if len(Mt.shape)==2:    RV = RV[:, :, 0]
        return RV