Matrix triangular solve op.
979000
def MatrixTriangularSolve(a, rhs, lower, adj):
    
    trans = 0 if not adj else 2

    r = np.empty(rhs.shape).astype(a.dtype)
    for coord in np.ndindex(a.shape[:-2]):
        pos = coord + (Ellipsis,)
        r[pos] = sp.linalg.solve_triangular(a[pos] if not adj else np.conj(a[pos]), rhs[pos],
                                            trans=trans, lower=lower)

    return r,
1347819
function triaUpSolve(A, b) {
    var size = A[0].length;
    var x = jStat.zeros(1, size)[0];
    var parts;
    var matrix_mode = false;

    if (b[0].length != undefined) {
      b = b.map(function(i){ return i[0] });
      matrix_mode = true;
    }

    jStat.arange(size - 1, -1, -1).forEach(function(i) {
      parts = jStat.arange(i + 1, size).map(function(j) {
        return x[j] * A[i][j];
      });
      x[i] = (b[i] - jStat.sum(parts)) / A[i][i];
    });

    if (matrix_mode)
      return x.map(function(i){ return [i] });
    return x;
  }
897702
public static FloatMatrix solveSymmetric(FloatMatrix A, FloatMatrix B) {
		A.assertSquare();
		FloatMatrix X = B.dup();
		int[] ipiv = new int[B.rows];
		SimpleBlas.sysv('U', A.dup(), ipiv, X);
		return X;
	}
970521
public function solve($B)
    {
        if ($B->getRowDimension() == $this->m) {
            if ($this->isNonsingular()) {
                // Copy right hand side with pivoting
                $nx = $B->getColumnDimension();
                $X  = $B->getMatrix($this->piv, 0, $nx-1);
                // Solve L*Y = B(piv,:)
                for ($k = 0; $k < $this->n; $k++) {
                    for ($i = $k+1; $i < $this->n; $i++) {
                        for ($j = 0; $j < $nx; $j++) {
                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                        }
                    }
                }
                // Solve U*X = Y;
                for ($k = $this->n-1; $k >= 0; $k--) {
                    for ($j = 0; $j < $nx; $j++) {
                        $X->A[$k][$j] /= $this->LU[$k][$k];
                    }
                    for ($i = 0; $i < $k; $i++) {
                        for ($j = 0; $j < $nx; $j++) {
                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                        }
                    }
                }
                return $X;
            } else {
                trigger_error(MatrixSingularException, ERROR);
            }
        } else {
            trigger_error(MatrixSquareException, ERROR);
        }
    }
378501
def _solve(self, x0, A, l, u, xmin, xmax):
        
        # Indexes of constrained lines.
        il = [i for i,ln in enumerate(self._ln) if 0.0 < ln.rate_a < 1e10]
        nl2 = len(il)

        neqnln = 2 * self._nb # no. of non-linear equality constraints
        niqnln = 2 * len(il)  # no. of lines with constraints

        user_data = {"A": A, "neqnln": neqnln, "niqnln": niqnln}

        self._f(x0)
        Jdata = self._dg(x0, False, user_data)
#        Hdata = self._h(x0, ones(neqnln + niqnln), None, False, user_data)

        lmbda = {"eqnonlin": ones(neqnln),
                 "ineqnonlin": ones(niqnln)}
        H = tril(self._hessfcn(x0, lmbda), format="coo")
        self._Hrow, self._Hcol = H.row, H.col

        n = len(x0) # the number of variables
        xl = xmin
        xu = xmax
        gl = r_[zeros(2 * self._nb), -Inf * ones(2 * nl2), l]
        gu = r_[zeros(2 * self._nb),       zeros(2 * nl2), u]
        m = len(gl) # the number of constraints
        nnzj = len(Jdata) # the number of nonzeros in Jacobian matrix
        nnzh = 0#len(H.data) # the number of non-zeros in Hessian matrix

        f_fcn, df_fcn, g_fcn, dg_fcn, h_fcn = \
            self._f, self._df, self._g, self._dg, self._h

        nlp = pyipopt.create(n, xl, xu, m, gl, gu, nnzj, nnzh,
                             f_fcn, df_fcn, g_fcn, dg_fcn)#, h_fcn)

#        print dir(nlp)
#        nlp.str_option("print_options_documentation", "yes")
#        nlp.int_option("max_iter", 10)

#        x, zl, zu, obj = nlp.solve(x0)
        success = nlp.solve(x0, user_data)
        nlp.close()
1199464
def dot_solve(self, y):
        

        """
        return np.dot(y.T, cho_solve(self._factor, y))
898384
def _cho_solve_AATI(A, rho, b, c, lwr, check_finite=True):
    """"""

    N, M = A.shape
    if N >= M:
        x = (b - _cho_solve((c, lwr), b.dot(A).T,
            check_finite=check_finite).T.dot(A.T)) / rho
    else:
        x = _cho_solve((c, lwr), b.T, check_finite=check_finite).T
    return x
908829
@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        if( B.numRows != numRows )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numRows);
        X.reshape(numCols,B.numCols);

        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++ ) {
                a.data[i] = B.data[i*BnumCols + colB];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                double []u = QR[n];

                double vv = u[n];
                u[n] = 1;
                QrHelperFunctions_DDRM.rank1UpdateMultR(a, u, gammas[n], 0, n, numRows, temp.data);
                u[n] = vv;
            }

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

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.data[i*X.numCols+colB] = a.data[i];
            }
        }
    }
1354194
def correct(self, z):
        ''''''
        c = self.linear_system.Ml*(
            self.linear_system.b - self.linear_system.A*z)
        c = utils.inner(self.W, c, ip_B=self.ip_B)
        if self.Q is not None and self.R is not None:
            c = scipy.linalg.solve_triangular(self.R, self.Q.T.conj().dot(c))
        if self.WR is not self.VR:
            c = self.WR.dot(scipy.linalg.solve_triangular(self.VR, c))
        return z + self.W.dot(c)
886371
def xstep(self):
        
        """

        YU = self.Y - self.U
        self.X = np.asarray(sl.cho_solve_ATAI(
            self.D, 1.0, self.block_sep0(YU) +
            self.D.T.dot(self.block_sep1(YU)), self.lu, self.piv),
            dtype=self.dtype)
205648
def MatrixSolve(a, rhs, adj):
    
    return np.linalg.solve(a if not adj else _adjoint(a), rhs),
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
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)
1035795
def solve(self):
        
        base_mva = self.om.case.base_mva
        Bf = self.om._Bf
        Pfinj = self.om._Pfinj
        # Unpack the OPF model.
        bs, ln, gn, cp = self._unpack_model(self.om)
        # Compute problem dimensions.
        ipol, ipwl, nb, nl, nw, ny, nxyz = self._dimension_data(bs, ln, gn)
        # Split the constraints in equality and inequality.
        AA, ll, uu = self._linear_constraints(self.om)
        # Piece-wise linear components of the objective function.
        Npwl, Hpwl, Cpwl, fparm_pwl, any_pwl = self._pwl_costs(ny, nxyz, ipwl)
        # Quadratic components of the objective function.
        Npol, Hpol, Cpol, fparm_pol, polycf, npol = \
            self._quadratic_costs(gn, ipol, nxyz, base_mva)
        # Combine pwl, poly and user costs.
        NN, HHw, CCw, ffparm = \
            self._combine_costs(Npwl, Hpwl, Cpwl, fparm_pwl, any_pwl,
                                Npol, Hpol, Cpol, fparm_pol, npol, nw)
        # Transform quadratic coefficients for w into coefficients for X.
        HH, CC, C0 = self._transform_coefficients(NN, HHw, CCw, ffparm, polycf,
                                                  any_pwl, npol, nw)
        # Bounds on the optimisation variables.
        _, xmin, xmax = self._var_bounds()

        # Select an interior initial point for interior point solver.
        x0 = self._initial_interior_point(bs, gn, xmin, xmax, ny)

        # Call the quadratic/linear solver.
        s = self._run_opf(HH, CC, AA, ll, uu, xmin, xmax, x0, self.opt)

        # Compute the objective function value.
        Va, Pg = self._update_solution_data(s, HH, CC, C0)

        # Set case result attributes.
        self._update_case(bs, ln, gn, base_mva, Bf, Pfinj, Va, Pg, s["lmbda"])

        return s
702712
def solve_kkt(U_Q, d, G, A, U_S, rx, rs, rz, ry, dbg=False):
    
    nineq, nz, neq, _ = get_sizes(G, A)

    invQ_rx = torch.potrs(rx.view(-1, 1), U_Q).view(-1)
    if neq > 0:
        h = torch.cat([torch.mv(A, invQ_rx) - ry,
                       torch.mv(G, invQ_rx) + rs / d - rz], 0)
    else:
        h = torch.mv(G, invQ_rx) + rs / d - rz

    w = -torch.potrs(h.view(-1, 1), U_S).view(-1)

    g1 = -rx - torch.mv(G.t(), w[neq:])
    if neq > 0:
        g1 -= torch.mv(A.t(), w[:neq])
    g2 = -rs - w[neq:]

    dx = torch.potrs(g1.view(-1, 1), U_Q).view(-1)
    ds = g2 / d
    dz = w[neq:]
    dy = w[:neq] if neq > 0 else None

    # if np.all(np.array([x.norm() for x in [rx, rs, rz, ry]]) != 0):
    if dbg:
        import IPython
        import sys
        IPython.embed()
        sys.exit(-1)

    # if rs.norm() > 0: import IPython, sys; IPython.embed(); sys.exit(-1)
    return dx, ds, dz, dy
243161
def solve_factorized_aug(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """RHS for condensed system"""
    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))
    dxnu = mysolve(LU, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
71508
def inv_posdef(mat):
    """"""
    # See:
    # - http://www.seas.ucla.edu/~vandenbe/103/lectures/chol.pdf
    # - http://scicomp.stackexchange.com/questions/3188
    chol = np.linalg.cholesky(mat)
    ident = np.eye(mat.shape[0])
    res = solve_triangular(chol, ident, lower=True, overwrite_b=True)
    return np.transpose(res).dot(res)