r"""
        Solves a triangular system of equations with multiple righthand
        sides. Computes

        .. math::

            B &:= L^{-1} B  \text{ if trans is 'N'}

           B &:= L^{-T} B  \text{ if trans is 'T'}
558672
def trsm(self,B,trans='N'):
         
        """
        
        if trans=='N':
            cp.trsm(self._L0,B)
            pftrsm(self._V,self._L,self._B,B,trans='N')
        elif trans=='T':
            pftrsm(self._V,self._L,self._B,B,trans='T')
            cp.trsm(self._L0,B,trans='T')
        elif type(trans) is str:
            raise ValueError("trans must be 'N' or 'T'")
        else:
            raise TypeError("trans must be 'N' or 'T'")
        return
352578
def solve(A, b):
    
    """
    A = asarray(A, float)
    b = asarray(b, float)
    if A.shape[0] == 1:

        with errstate(divide="ignore"):
            A_ = array([[1.0 / A[0, 0]]])

        if not isfinite(A_[0, 0]):
            raise LinAlgError("Division error.")

        return dot(A_, b)
    elif A.shape[0] == 2:
        a = A[0, 0]
        b_ = A[0, 1]
        c = A[1, 0]
        d = A[1, 1]
        A_ = array([[d, -b_], [-c, a]])

        with errstate(divide="ignore"):
            A_ /= a * d - b_ * c

        if not npy_all(isfinite(A_)):
            raise LinAlgError("Division error.")

        return dot(A_, b)
    return _solve(A, b)
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)
1269769
public static DoubleMatrix solveSymmetric(DoubleMatrix A, DoubleMatrix B) {
		A.assertSquare();
		DoubleMatrix X = B.dup();
		int[] ipiv = new int[B.rows];
		SimpleBlas.sysv('U', A.dup(), ipiv, X);
		return X;
	}
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
787006
public function solve($B)
    {
        if ($B->getRowDimension() == $this->m) {
            if ($this->isFullRank()) {
                // Copy right hand side
                $nx = $B->getColumnDimension();
                $X  = $B->getArrayCopy();
                // Compute Y = transpose(Q)*B
                for ($k = 0; $k < $this->n; $k++) {
                    for ($j = 0; $j < $nx; $j++) {
                        $s = 0.0;
                        for ($i = $k; $i < $this->m; $i++) {
                            $s += $this->QR[$i][$k] * $X[$i][$j];
                        }
                        $s = -$s/$this->QR[$k][$k];
                        for ($i = $k; $i < $this->m; $i++) {
                            $X[$i][$j] += $s * $this->QR[$i][$k];
                        }
                    }
                }
                // Solve R*X = Y;
                for ($k = $this->n-1; $k >= 0; $k--) {
                    for ($j = 0; $j < $nx; $j++) {
                        $X[$k][$j] /= $this->Rdiag[$k];
                    }
                    for ($i = 0; $i < $k; $i++) {
                        for ($j = 0; $j < $nx; $j++) {
                            $X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
                        }
                    }
                }
                $X = new Matrix($X);
                return ($X->getMatrix(0, $this->n-1, 0, $nx));
            } else {
                trigger_error(MatrixRankException, ERROR);
            }
        } else {
            trigger_error(MatrixDimensionException, ERROR);
        }
    }
237726
def solve(self,b):
        

        x = b.copy()
        self.mumps.set_rhs(x)
        self.mumps.run(job=3)

        return x
190026
private double[][] solveLtransposed(double[][] X) {
    for(int k = L.length - 1; k >= 0; k--) {
      final double[] Lk = L[k], Xk = X[k];
      timesEquals(Xk, 1. / Lk[k]);
      for(int i = 0; i < k; i++) {
        minusTimesEquals(X[i], Xk, Lk[i]);
      }
    }
    return X;
  }
2344
def solver(AA, N_max, symNx = 2, throw_out_modes=False):
    

    # Find all integer component n_vectors which lie within sphere of radius N_max
    # Here we have assumed that the potential is symmetric x->-x, y->-y, z->-z
    # This can be relaxed by changing symN to 1
    # Additionally due to time reversal symmetry S_n = -S_-n so we only consider
    # "half" of the n-vector-space

    angs = unroll_angles(AA.T[3:].T,np.ones(3))

    symNz = 2
    NNx = range(-N_max, N_max+1, symNx)
    NNy = range(-N_max, N_max+1, symNz)
    NNz = range(-N_max, N_max+1, symNz)
    n_vectors = np.array([[i,j,k] for (i,j,k) in product(NNx,NNy,NNz)
                          if(not(i==0 and j==0 and k==0)            # exclude zero vector
                             and (k>0                               # northern hemisphere
                                  or (k==0 and j>0)                 # half of x-y plane
                                  or (k==0 and j==0 and i>0))       # half of x axis
                             and np.sqrt(i*i+j*j+k*k)<=N_max)])     # inside sphere

    xxx = check_each_direction(n_vectors,angs)

    if(throw_out_modes):
        n_vectors = np.delete(n_vectors,check_each_direction(n_vectors,angs),axis=0)

    n = len(n_vectors)+3
    b = np.zeros(shape=(n, ))
    a = np.zeros(shape=(n,n))

    a[:3,:3]=len(AA)*np.identity(3)

    for i in AA:
        a[:3,3:]+=2.*n_vectors.T[:3]*np.cos(np.dot(n_vectors,i[3:]))
        a[3:,3:]+=4.*np.dot(n_vectors,n_vectors.T)*np.outer(np.cos(np.dot(n_vectors,i[3:])),np.cos(np.dot(n_vectors,i[3:])))
        b[:3]+=i[:3]
        b[3:]+=2.*np.dot(n_vectors,i[:3])*np.cos(np.dot(n_vectors,i[3:]))

    a[3:,:3]=a[:3,3:].T

    return np.array(solve(a,b)), n_vectors
760169
def factorize_and_solve(self, A, b):
        

        self.factorize(A)

        return self.solve(b)
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;
        }
    }
1318529
function(eqns, solve_for, solutions) { 
        if(eqns instanceof Equation) {
            //if it's zero then we're done
            if(eqns.isZero())
                return [new Symbol(0)];
            //if the lhs = x then we're done
            if(eqns.LHS.equals(solve_for) && !eqns.RHS.contains(solve_for)) {
                return [eqns.RHS];
            }
            //if the rhs = x then we're done
            if(eqns.RHS.equals(solve_for) && !eqns.LHS.contains(solve_for)) {
                return [eqns.LHS];
            }
        }
        //unwrap the vector since what we want are the elements
        if(eqns instanceof core.Vector)
            eqns = eqns.elements;
        solve_for = solve_for || 'x'; //assumes x by default
        //If it's an array then solve it as a system of equations
        if(isArray(eqns)) {
            return sys_solve.apply(undefined, arguments);
        }
        //parse out functions. Fix for issue #300
        //eqns = core.Utils.evaluate(eqns);
        
        solutions = solutions || [];
        //maybe we get lucky
        if(eqns.group === S && eqns.contains(solve_for)) {
            solutions.push(new Symbol(0));
            return solutions;
        }  
        if(eqns.group === CB) {
            var sf = String(solve_for); //everything else belongs to the coeff
            eqns.each(function(x) {
                if(x.contains(sf))
                    solve(x, solve_for, solutions);
            });
 
            return solutions;
        }
        var existing = {}, //mark existing solutions as not to have duplicates
            add_to_result = function(r, has_trig) {
                var r_is_symbol = isSymbol(r);
                if(r === undefined || typeof r === 'number' && isNaN(r))
                    return;
                if(isArray(r)) 
                    solutions = solutions.concat(r);
                else { 
                    if(r.valueOf() !== 'null') {
                        if(!r_is_symbol)
                            r = _.parse(r);
                        //try to convert the number to multiples of pi
                        if(core.Settings.make_pi_conversions && has_trig) {
                            var temp = _.divide(r.clone(), new Symbol(Math.PI)),
                                m = temp.multiplier,
                                a = Math.abs(m.num),
                                b = Math.abs(m.den);
                            if(a < 10 && b < 10)
                                r = _.multiply(temp, new Symbol('pi'));
                        }
                        //convert to a string so we can mark it as a known solution
                        var r_str = r.toString();
                        if(!existing[r_str])
                            solutions.push(r);
                        //mark the answer as seen
                        existing[r_str] = true;
                    }
                }
            };
        //gets points around which to solve. It does that because it builds on the principle that if
        //the sign changes over an interval then there must be a zero on that interval
        
        
        var attempt_Newton = function(symbol) { 
            var has_trig = symbol.hasTrig();
            // we get all the points where a possible zero might exist
            var points1 = get_points(symbol, 0.1);
            var points2 = get_points(symbol, 0.05);
            var points3 = get_points(symbol, 0.01);
            var points = core.Utils.arrayUnique(points1.concat(points2).concat(points3)),
                l = points.length;
            //compile the function and the derivative of the function
            var f = build(symbol.clone()),
                fp = build(_C.diff(symbol.clone()));
            for(var i=0; i<l; i++) {
                var point = points[i];
                add_to_result(Newton(point, f, fp), has_trig);
            }
            solutions.sort();
        };

        var eq = core.Utils.isSymbol(eqns) ? eqns : toLHS(eqns),
            vars = core.Utils.variables(eq),//get a list of all the variables
            numvars = vars.length;//how many variables are we dealing with
        //if we're dealing with a single variable then we first check if it's a 
        //polynomial (including rationals).If it is then we use the Jenkins-Traubb algorithm.     
        //Don't waste time
        if(eq.group === S || eq.group === CB && eq.contains(solve_for))
            return [new Symbol(0)];
        //force to polynomial. We go through each and then we look at what it would 
        //take for its power to be an integer
        //if the power is a fractional we divide by the fractional power
        var fractionals = {},
            cfact;
        var correct_denom = function(symbol) { 
            var original = symbol.clone(); //preserve the original
            if(symbol.symbols) {
                for(var x in symbol.symbols) { 
                    var sym = symbol.symbols[x];
                    var parts = explode(sym, solve_for);
                    var is_sqrt = parts[1].fname === core.Settings.SQRT;
                    var v = Symbol.unwrapSQRT(parts[1]);
                    var p = v.power.clone();
                    //circular logic with sqrt. Since sqrt(x) becomes x^(1/2) which then becomes sqrt(x), this continues forever
                    //this needs to be terminated if p = 1/2
                    if(!isSymbol(p) && !p.equals(1/2)) {
                        if(p.den.gt(1)) { 
                            if(is_sqrt) {
                                symbol = _.subtract(symbol, sym.clone());
                                symbol = _.add(symbol, _.multiply(parts[0].clone(), v));
                                return correct_denom(symbol);
                            }
                            var c = fractionals[p.den];
                            fractionals[p.den] = c ? c++ : 1;
                        }
                        else if(p.sign() === -1){ 
                            var factor = _.parse(solve_for+'^'+Math.abs(p)); //this
                            //unwrap the symbol's denoniator
                            symbol.each(function(y, index) {
                                if(y.contains(solve_for)) {
                                    symbol.symbols[index] = _.multiply(y, factor.clone());
                                } 
                            });
                            fractionals = {};
                            return correct_denom(_.parse(symbol));
                        }
                        else if(sym.group === PL) {
                            var min_p = core.Utils.arrayMin(core.Utils.keys(sym.symbols));
                            if(min_p < 0) {
                                var factor = _.parse(solve_for+'^'+Math.abs(min_p));
                                var corrected = new Symbol(0);
                                original.each(function(x) {
                                    corrected = _.add(corrected, _.multiply(x.clone(), factor.clone()));
                                }, true);
                                return corrected;
                            }
                        }
                    }
                }
            }
            return symbol;
        };
        //rewrites equations/expression in simpler form
        var rewrite = function(rhs, lhs) { 
            lhs = lhs || new Symbol(0);
            rhs = Symbol.unwrapSQRT(_.expand(rhs)); //expand the term expression go get rid of quotients when possible
            var c = 0, //a counter to see if we have all terms with the variable
                l = rhs.length;
            //try to rewrite the whole thing
            if(rhs.group === CP && rhs.contains(solve_for) && rhs.isLinear()) {
                rhs.distributeMultiplier();
                var t = new Symbol(0);
                //first bring all the terms containing the variable to the lhs
                rhs.each(function(x) {
                    if(x.contains(solve_for)) {
                        c++;
                        t = _.add(t, x.clone());
                    }
                    else
                        lhs = _.subtract(lhs, x.clone());
                });
                rhs = t;
                
                //if not all the terms contain the variable so it's in the form
                //a*x^2+x
                if(c !== l)
                    return rewrite(rhs, lhs);
                else { 
                    return [rhs, lhs];
                }
            }
            else if(rhs.group === CB && rhs.contains(solve_for) && rhs.isLinear()) {
                if(rhs.multiplier.lessThan(0)) {
                    rhs.multiplier = rhs.multiplier.multiply(new core.Frac(-1));
                    lhs.multiplier = lhs.multiplier.multiply(new core.Frac(-1));
                }
                if(lhs.equals(0))
                    return new Symbol(0);
                else {
                    var t = new Symbol(1);
                    rhs.each(function(x) { 
                        if(x.contains(solve_for)) 
                            t = _.multiply(t, x.clone());
                        else 
                            lhs = _.divide(lhs, x.clone());
                    });
                    rhs = t;
                    return rewrite(rhs, lhs);
                    
                }
            }   
            else if(!rhs.isLinear() && rhs.contains(solve_for)) { 
                var p = _.parse(rhs.power.clone().invert());
                rhs = _.pow(rhs, p.clone());
                lhs = _.pow(_.expand(lhs), p.clone());
                return rewrite(rhs, lhs);
            }
            else if(rhs.group === FN || rhs.group === S || rhs.group === PL) {
                return [rhs, lhs];
            }
        };

        //separate the equation
        var separate = function(eq) {
            var lhs = new Symbol(0),
                rhs = new Symbol(0);
            eq.each(function(x) {
                if(x.contains(solve_for, true))
                    lhs = _.add(lhs, x.clone());
                else
                    rhs = _.subtract(rhs, x.clone());
            });
            return [lhs, rhs];
        };
        
        //first remove any denominators
        eq = correct_denom(eq);  

        if(eq.equals(0))
            return [eq];
        //correct fractionals. I can only handle one type right now
        var fkeys = core.Utils.keys(fractionals);
        if(fkeys.length === 1) { 
            //make a note of the factor
            cfact = fkeys[0];
            eq.each(function(x, index) {
                if(x.contains(solve_for)) {
                    var parts = explode(x, solve_for);
                    var v = parts[1];
                    var p = v.power;
                    if(p.den.gt(1)) { 
                        v.power = p.multiply(new core.Frac(cfact));
                        eq.symbols[index] = _.multiply(v, parts[0]);
                    }
                } 
            });
            eq = _.parse(eq);
        }
        //polynomial single variable
        if(numvars === 1) { 
            if(eq.isPoly(true)) { 
                var coeffs = core.Utils.getCoeffs(eq, solve_for),
                    deg = coeffs.length - 1,
                    was_calculated = false;
                if(vars[0] === solve_for) {
                    //check to see if all the coefficients are constant
                    if(checkAll(coeffs, function(x) {
                        return x.group !== core.groups.N;
                    })) {
                        var roots = core.Algebra.proots(eq);
                        //if all the roots are integers then return those
                        if(checkAll(roots, function(x) {
                            return !core.Utils.isInt(x);
                        })) {
                            //roots have been calculates
                            was_calculated = true;
                            roots.map(function(x) {
                                add_to_result(new Symbol(x));
                            });
                        }
                    }
                    
                    if(!was_calculated) {
                        //we can solve algebraically for degrees 1, 2, 3. The remainder we switch to Jenkins-
                        if(deg === 1) 
                            add_to_result(_.divide(coeffs[0], coeffs[1].negate()));
                        else if(deg === 2) { 
                            add_to_result(_.expand(quad.apply(undefined, coeffs)));
                            coeffs.push('-');
                            add_to_result(_.expand(quad.apply(undefined, coeffs)));
                        }
                        else if(deg === 3)
                            add_to_result(cubic.apply(undefined, coeffs));
                        else {
                            /*
                            var sym_roots = csolve(eq, solve_for); 
                            if(sym_roots.length === 0)
                                sym_roots = divnconsolve(eq, solve_for);
                            if(sym_roots.length > 0) 
                                add_to_result(sym_roots);
                            else
                            */
                            _A.proots(eq).map(add_to_result);
                        }
                    }
                            
                }
            }
            else { 
                //since it's not a polynomial then we'll try to look for a solution using Newton's method
                //this is not a very broad search but takes the positions that something is better than nothing
                attempt_Newton(eq);
            }
        }
        else { 
            //The idea here is to go through the equation and collect the coefficients
            //place them in an array and call the quad or cubic function to get the results
            if(!eq.hasFunc(solve_for) && eq.isComposite()) { 
                try {
                    var coeffs = core.Utils.getCoeffs(eq, solve_for);
                    
                    var l = coeffs.length,
                        deg = l-1; //the degree of the polynomial
                    //handle the problem based on the degree
                    switch(deg) {
                        case 0:
                            var separated = separate(eq);
                            var lhs = separated[0],
                                rhs = separated[1];
                            if(lhs.group === core.groups.EX) {
                                add_to_result(_.parse(core.Utils.format('log(({0})/({2}))/log({1})', rhs, lhs.value, lhs.multiplier)));
                            }
                            break;
                        case 1:
                            //nothing to do but to return the quotient of the constant and the LT
                            //e.g. 2*x-1
                            add_to_result(_.divide(coeffs[0], coeffs[1].negate()));
                            break;
                        case 2:
                            add_to_result(quad.apply(undefined, coeffs));
                            coeffs.push('-');
                            add_to_result(quad.apply(undefined, coeffs));
                            break;
                        case 3:
                            add_to_result(cubic.apply(undefined, coeffs));
                            break;
                        case 4:
                            add_to_result(quartic.apply(undefined, coeffs));
                            break;
                        default:
                            add_to_result(csolve(eq, solve_for));
                            if(solutions.length === 0)
                                add_to_result(divnconsolve(eq, solve_for));
                    }
                }
                catch(e) { /*something went wrong. EXITING*/; } 
            }
            else {
                try {
                    var rw = rewrite(eq);
                    var lhs = rw[0];
                    var rhs = rw[1];
                    if(lhs.group === FN) {
                        if(lhs.fname === 'abs') {
                            solutions.push(rhs.clone());
                            solutions.push(rhs.negate());
                        }
                        else
                            solutions.push(_.subtract(lhs, rhs));
                    }
                }
                catch(error) {
                    //Let's try this another way
                    try { 
                        //1. if the symbol is in the form a*b*c*... then the solution is zero if 
                        //either a or b or c is zero.
                        if(eq.group === CB)
                            solutions.push(0);
                        else if(eq.group === CP) {
                            var separated = separate(eq);
                            var lhs = separated[0],
                                rhs = separated[1];
                            
                            //reduce the equation
                            if(lhs.group === core.groups.EX && lhs.value === solve_for) {
                                //change the base of both sides
                                var p = lhs.power.clone().invert();
                                solutions.push(_.pow(rhs, p));
                            }
                        }
                    }
                    catch(error){;}
                }
            }
        }
        
        if(cfact) {
            solutions = solutions.map(function(x) {
                return _.pow(x, new Symbol(cfact));
            });
        }
        
        return solutions;
    }
279061
def B(g,i):
    """"""
    g1 = g&(2**i)
    if i:
        nA = Awidth(i)
        nB = Bwidth(i)
        i=i-1
        Bn = B(g,i)
        if g1:
            return Bn            << (nA+nB) | int('1'*nA,2) << nB | Bn
        else:
            return int('1'*nB,2) << (nA+nB) | A(g,i)      << nB | Bn
    else:
        if g1:
            return 1
        else:
            return 0
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
679061
public static void multTranA( DMatrix3x3 A , DMatrix3x3 B , DMatrix3x3 C , DMatrix3x3 output )
	{
		double t11 = A.a11*B.a11 + A.a21*B.a21 + A.a31*B.a31;
		double t12 = A.a11*B.a12 + A.a21*B.a22 + A.a31*B.a32;
		double t13 = A.a11*B.a13 + A.a21*B.a23 + A.a31*B.a33;

		double t21 = A.a12*B.a11 + A.a22*B.a21 + A.a32*B.a31;
		double t22 = A.a12*B.a12 + A.a22*B.a22 + A.a32*B.a32;
		double t23 = A.a12*B.a13 + A.a22*B.a23 + A.a32*B.a33;

		double t31 = A.a13*B.a11 + A.a23*B.a21 + A.a33*B.a31;
		double t32 = A.a13*B.a12 + A.a23*B.a22 + A.a33*B.a32;
		double t33 = A.a13*B.a13 + A.a23*B.a23 + A.a33*B.a33;

		output.a11 = t11*C.a11 + t12*C.a21 + t13*C.a31;
		output.a12 = t11*C.a12 + t12*C.a22 + t13*C.a32;
		output.a13 = t11*C.a13 + t12*C.a23 + t13*C.a33;

		output.a21 = t21*C.a11 + t22*C.a21 + t23*C.a31;
		output.a22 = t21*C.a12 + t22*C.a22 + t23*C.a32;
		output.a23 = t21*C.a13 + t22*C.a23 + t23*C.a33;

		output.a31 = t31*C.a11 + t32*C.a21 + t33*C.a31;
		output.a32 = t31*C.a12 + t32*C.a22 + t33*C.a32;
		output.a33 = t31*C.a13 + t32*C.a23 + t33*C.a33;
	}
524879
def cublasDtrsm(handle, side, uplo, trans, diag, m, n, alpha, A, lda, B, ldb):
    
    
    status = _libcublas.cublasDtrsm_v2(handle, 
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       m, n, ctypes.byref(ctypes.c_double(alpha)),
                                       int(A), lda, int(B), ldb)
    cublasCheckStatus(status)
1307818
def quadraticSolver(a, b, c):
    """"""
    if a == 0:
        if b == 0:  return [] # attempting to solve an equation with infinite (0=0) or impossible (0=3) solutions for x
        else:       return [-c / float(b)]
    else:
        d = b**2.0 - 4*a*c # calculate the discriminant
        if d < 0:   return [] # would require math.sqrt a negative number, yielding a complex number
        solvedPart0 = 1/(2.0*a)
        solvedPart1 = -b * solvedPart0
        if d == 0:  return [ solvedPart1 ] # two values at same 'x'
        else: # i.e. d > 0   apply full quadratic formula
            solvedPart2 = math.sqrt(d)*solvedPart0
            return [ solvedPart1+solvedPart2, solvedPart1-solvedPart2 ]