CUSOLVERRF (forward and backward triangular) solves
763289
public static int cusolverRfSolve(
        cusolverRfHandle handle, 
        Pointer P, 
        Pointer Q, 
        int nrhs, //only nrhs=1 is supported
        Pointer Temp, //of size ldt*nrhs (ldt>=n)
        int ldt, 
        /** Input/Output (in the device memory) */
        Pointer XF, 
        /** Input */
        int ldxf)
    {
        return checkResult(cusolverRfSolveNative(handle, P, Q, nrhs, Temp, ldt, XF, ldxf));
    }
594543
public static int cusolverRfBatchResetValues(
        int batchSize, 
        int n, 
        int nnzA, 
        Pointer csrRowPtrA, 
        Pointer csrColIndA, 
        Pointer csrValA_array, 
        Pointer P, 
        Pointer Q, 
        /** Output */
        cusolverRfHandle handle)
    {
        return checkResult(cusolverRfBatchResetValuesNative(batchSize, n, nnzA, csrRowPtrA, csrColIndA, csrValA_array, P, Q, handle));
    }
845218
public static int cusolverSpScsrlsvluHost(
        cusolverSpHandle handle, 
        int n, 
        int nnzA, 
        cusparseMatDescr descrA, 
        Pointer csrValA, 
        Pointer csrRowPtrA, 
        Pointer csrColIndA, 
        Pointer b, 
        float tol, 
        int reorder, 
        Pointer x, 
        int[] singularity)
    {
        return checkResult(cusolverSpScsrlsvluHostNative(handle, n, nnzA, descrA, csrValA, csrRowPtrA, csrColIndA, b, tol, reorder, x, singularity));
    }
699930
def _multicomplex2(f, fx, x, h):
        """"""
        n = len(x)
        ee = np.diag(h)
        hess = np.outer(h, h)
        cmplx_wrap = Bicomplex.__array_wrap__
        for i in range(n):
            for j in range(i, n):
                zph = Bicomplex(x + 1j * ee[i, :], ee[j, :])
                hess[i, j] = cmplx_wrap(f(zph)).imag12 / hess[j, i]
                hess[j, i] = hess[i, j]
        return hess
1355444
def solve_dp(self, lam):
        ''''''
        cur_converge = self.converge+1
        step = 0
        # Get an initial estimate using the GFL
        self.solve_gfl(lam)
        beta2 = np.copy(self.beta)
        while cur_converge > self.converge and step < self.max_dp_steps:
            # Weight each edge differently
            u = lam / (1 + np.abs(self.beta[self.trails[::2]] - self.beta[self.trails[1::2]]))
            # Swap the beta buffers
            temp = self.beta
            self.beta = beta2
            beta2 = temp
            # Solve the edge-weighted GFL problem, which updates beta
            self.solve_gfl(u)
            # Check for convergence
            cur_converge = np.sqrt(((self.beta - beta2)**2).sum())
            step += 1
        self.steps.append(step)
        return self.beta
714193
private void ssMergeForward(int PA, int first, int middle, int last, int buf,
                                int depth) {
        // PA, first, middle, last, buf are pointers to SA
        int a, b, c, bufend;// pointers to SA
        int t, r;

        bufend = buf + (middle - first) - 1;
        ssBlockSwap(buf, first, middle - first);

        for (t = SA[a = first], b = buf, c = middle; ; ) {
            r = ssCompare(PA + SA[b], PA + SA[c], depth);
            if (r < 0) {
                do {
                    SA[a++] = SA[b];
                    if (bufend <= b) {
                        SA[bufend] = t;
                        return;
                    }
                    SA[b++] = SA[a];
                }
                while (SA[b] < 0);
            } else if (r > 0) {
                do {
                    SA[a++] = SA[c];
                    SA[c++] = SA[a];
                    if (last <= c) {
                        while (b < bufend) {
                            SA[a++] = SA[b];
                            SA[b++] = SA[a];
                        }
                        SA[a] = SA[b];
                        SA[b] = t;
                        return;
                    }
                }
                while (SA[c] < 0);
            } else {
                SA[c] = ~SA[c];
                do {
                    SA[a++] = SA[b];
                    if (bufend <= b) {
                        SA[bufend] = t;
                        return;
                    }
                    SA[b++] = SA[a];
                }
                while (SA[b] < 0);

                do {
                    SA[a++] = SA[c];
                    SA[c++] = SA[a];
                    if (last <= c) {
                        while (b < bufend) {
                            SA[a++] = SA[b];
                            SA[b++] = SA[a];
                        }
                        SA[a] = SA[b];
                        SA[b] = t;
                        return;
                    }
                }
                while (SA[c] < 0);
            }
        }

    }
112739
public static void solveTranL(DMatrixSparseCSC L , double []x )
    {
        final int N = L.numCols;

        for (int j = N-1; j >= 0; j--) {
            int idx0 = L.col_idx[j];
            int idx1 = L.col_idx[j+1];

            for (int p = idx0+1; p < idx1; p++) {
                x[j] -= L.nz_values[p]*x[L.nz_rows[p]];
            }
            x[j] /= L.nz_values[idx0];
        }
    }
929959
def cublasZsyrk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasZsyrk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       n, k, ctypes.byref(cuda.cuDoubleComplex(alpha.real,    
                                                                               alpha.imag)),
                                       int(A), lda,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                         beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
1028630
public static int cusolverSpXcsrsymamdHost(
        cusolverSpHandle handle, 
        int n, 
        int nnzA, 
        cusparseMatDescr descrA, 
        Pointer csrRowPtrA, 
        Pointer csrColIndA, 
        Pointer p)
    {
        return checkResult(cusolverSpXcsrsymamdHostNative(handle, n, nnzA, descrA, csrRowPtrA, csrColIndA, p));
    }
988274
def cublasCher2k(handle, uplo, trans, n, k, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasCher2k_v2(handle, 
                                        _CUBLAS_FILL_MODE[uplo], 
                                        _CUBLAS_OP[trans], 
                                        n, k, ctypes.byref(cuda.cuFloatComplex(alpha.real,                 
                                                                               alpha.imag)),
                                        int(A), lda, int(B), ldb, 
                                        ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                         beta.imag)),
                                        int(C), ldc)
    cublasCheckStatus(status)
336215
public static int cusolverSpScsrlsvqrHost(
        cusolverSpHandle handle, 
        int m, 
        int nnz, 
        cusparseMatDescr descrA, 
        Pointer csrValA, 
        Pointer csrRowPtrA, 
        Pointer csrColIndA, 
        Pointer b, 
        float tol, 
        int reorder, 
        Pointer x, 
        int[] singularity)
    {
        return checkResult(cusolverSpScsrlsvqrHostNative(handle, m, nnz, descrA, csrValA, csrRowPtrA, csrColIndA, b, tol, reorder, x, singularity));
    }
1165685
private void solveInternalL() {
        // This takes advantage of the diagonal elements always being real numbers

        // solve L*y=b storing y in x
        TriangularSolver_ZDRM.solveL_diagReal(t, vv, n);

        // solve L^T*x=y
        TriangularSolver_ZDRM.solveConjTranL_diagReal(t, vv, n);
    }
600200
def solve(self):
        '''
        
        '''
        self.defUtilityFuncs()
        self.DiscFacEff = self.DiscFac*self.LivPrb
        self.makePFcFunc()
        self.defValueFuncs()
        solution = ConsumerSolution(cFunc=self.cFunc, vFunc=self.vFunc, vPfunc=self.vPfunc,
                                    mNrmMin=self.mNrmMin, hNrm=self.hNrmNow,
                                    MPCmin=self.MPC, MPCmax=self.MPC)
        #solution = self.addSSmNrm(solution)
        return solution
552428
def cublasCgemv(handle, trans, m, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasCgemv_v2(handle,
                                       _CUBLAS_OP[trans], m, n,
                                       ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)
124553
def cublasCsyr(handle, uplo, n, alpha, x, incx, A, lda):
    

    status = _libcublas.cublasCsyr_v2(handle,
                                        _CUBLAS_FILL_MODE[uplo], n, 
                                        ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                        int(x), incx, int(A), lda)
    cublasCheckStatus(status)
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)
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];
            }
        }
    }