Rank-k operation on real symmetric matrix.
776553
def cublasDsyrk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasDsyrk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       n, k, ctypes.byref(cuda.cuFloatComplex(alpha.real,     
                                                                              alpha.imag)),
                                       int(A), lda, 
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
683347
def K(self, parm):
        
        return OU_K_matrix(self.X, parm) + np.identity(self.X.shape[0])*(10**-10)
371693
private void _symmetricTql2 ()
  {
    // This is derived from the Algol procedures tql2, by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    for (int i = 1; i < m_nDim; i++)
      m_aEVe[i - 1] = m_aEVe[i];
    m_aEVe[m_nDim - 1] = 0.0;

    double f = 0.0;
    double tst1 = 0.0;
    for (int l = 0; l < m_nDim; l++)
    {
      // Find small subdiagonal element
      tst1 = Math.max (tst1, MathHelper.abs (m_aEVd[l]) + MathHelper.abs (m_aEVe[l]));
      int m = l;
      while (m < m_nDim)
      {
        if (MathHelper.abs (m_aEVe[m]) <= EPSILON * tst1)
          break;
        m++;
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.

      if (m > l)
      {
        do
        {
          // Compute implicit shift
          double g = m_aEVd[l];
          double p = (m_aEVd[l + 1] - g) / (2.0 * m_aEVe[l]);
          double r = MathHelper.hypot (p, 1.0);
          if (p < 0)
          {
            r = -r;
          }
          m_aEVd[l] = m_aEVe[l] / (p + r);
          m_aEVd[l + 1] = m_aEVe[l] * (p + r);
          final double dl1 = m_aEVd[l + 1];
          double h = g - m_aEVd[l];
          for (int i = l + 2; i < m_nDim; i++)
          {
            m_aEVd[i] -= h;
          }
          f += h;

          // Implicit QL transformation.
          p = m_aEVd[m];
          double c = 1.0;
          double c2 = c;
          double c3 = c;
          final double el1 = m_aEVe[l + 1];
          double s = 0.0;
          double s2 = 0.0;
          for (int i = m - 1; i >= l; i--)
          {
            c3 = c2;
            c2 = c;
            s2 = s;
            g = c * m_aEVe[i];
            h = c * p;
            r = MathHelper.hypot (p, m_aEVe[i]);
            m_aEVe[i + 1] = s * r;
            s = m_aEVe[i] / r;
            c = p / r;
            p = c * m_aEVd[i] - s * g;
            m_aEVd[i + 1] = h + s * (c * g + s * m_aEVd[i]);

            // Accumulate transformation.

            for (int k = 0; k < m_nDim; k++)
            {
              h = m_aEigenVector[k][i + 1];
              m_aEigenVector[k][i + 1] = s * m_aEigenVector[k][i] + c * h;
              m_aEigenVector[k][i] = c * m_aEigenVector[k][i] - s * h;
            }
          }
          p = -s * s2 * c3 * el1 * m_aEVe[l] / dl1;
          m_aEVe[l] = s * p;
          m_aEVd[l] = c * p;

          // Check for convergence.

        } while (MathHelper.abs (m_aEVe[l]) > EPSILON * tst1);
      }
      m_aEVd[l] = m_aEVd[l] + f;
      m_aEVe[l] = 0.0;
    }

    // Sort eigenvalues and corresponding vectors.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      int k = i;
      double p = m_aEVd[i];
      for (int j = i + 1; j < m_nDim; j++)
      {
        if (m_aEVd[j] < p)
        {
          k = j;
          p = m_aEVd[j];
        }
      }
      if (k != i)
      {
        m_aEVd[k] = m_aEVd[i];
        m_aEVd[i] = p;
        for (int j = 0; j < m_nDim; j++)
        {
          p = m_aEigenVector[j][i];
          m_aEigenVector[j][i] = m_aEigenVector[j][k];
          m_aEigenVector[j][k] = p;
        }
      }
    }
  }
650218
def solve_factorized_schur(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)

    """Assemble right hand side of augmented system"""
    r1 = rd + mydot(G.T, mydot(SIG, rp2)) - mydot(G.T, rc/s)
    r2 = rp1

    """Unpack LU-factors"""
    LU_S, LU_H = LU

    """Assemble right hand side for normal equation"""
    b = r2 - mydot(A, mysolve(LU_H, r1))  

    """Solve for dnu"""
    dnu = mysolve(LU_S, b)
       
    """Solve for dx"""
    dx = mysolve(LU_H, -(r1 + mydot(A.T, dnu)))    
    
    """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
466222
def kdists(matrix, k=7, ix=None):
    

    ix = ix or kindex(matrix, k)
    return matrix[ix][np.newaxis].T
738950
static void computePseudo(DMatrixRMaj A , DMatrixRMaj output ) {
		final int N = A.numRows;

		DMatrix2x2 m = new DMatrix2x2();
		DMatrix2x2 m_inv = new DMatrix2x2();

		for (int i = 0, index=0; i < N; i++) {
			double a_i1 = A.data[index++];
			double a_i2 = A.data[index++];
			m.a11 += a_i1*a_i1;
			m.a12 += a_i1*a_i2;
			m.a22 += a_i2*a_i2;
		}
		m.a21 = m.a12;
		CommonOps_DDF2.invert(m,m_inv);

		output.reshape(2,N);
		for (int i = 0,index=0; i < N; i++) {
			output.data[i] = A.data[index++]*m_inv.a11 + A.data[index++]*m_inv.a12;
		}
		int end = 2*N;
		for (int i = N,A_index=0; i < end; i++) {
			output.data[i] = A.data[A_index++]*m_inv.a21 + A.data[A_index++]*m_inv.a22;
		}
	}
386102
def _lmder1_linear_full_rank(n, m, factor, target_fnorm1, target_fnorm2):
    """"""

    def func(params, vec):
        s = params.sum()
        temp = 2. * s / m + 1
        vec[:] = -temp
        vec[:params.size] += params

    def jac(params, jac):
        # jac.shape = (n, m) by LMDER standards
        jac.fill(-2. / m)
        for i in range(n):
            jac[i,i] += 1

    guess = np.ones(n) * factor

    #_lmder1_test(m, func, jac, guess)
    _lmder1_driver(m, func, jac, guess,
                   target_fnorm1, target_fnorm2,
                   [-1] * n)
490009
def precompute_sharp_round(nxk, nyk, xc, yc):
    

    # Create arrays for the two- and four-fold symmetry computations:
    s4m = np.ones((nyk,nxk),dtype=np.int16)
    s4m[yc, xc] = 0

    s2m = np.ones((nyk,nxk),dtype=np.int16)
    s2m[yc, xc] = 0
    s2m[yc:nyk, 0:xc]     = -1;
    s2m[0:yc+1, xc+1:nxk] = -1;

    return s2m, s4m
1320422
function computeHOAMatrices(matrix) {
  // We start by computing the 2nd-order matrix from the 1st-order matrix.
  for (let i = 2; i <= matrix.length; i++) {
    computeBandRotation(matrix, i);
  }
}
32760
def _kak_decomposition_to_operations(q0: ops.Qid,
                                     q1: ops.Qid,
                                     kak: linalg.KakDecomposition,
                                     allow_partial_czs: bool,
                                     atol: float = 1e-8
                                     ) -> List[ops.Operation]:
    """"""
    b0, b1 = kak.single_qubit_operations_before
    pre = [_do_single_on(b0, q0, atol=atol), _do_single_on(b1, q1, atol=atol)]
    a0, a1 = kak.single_qubit_operations_after
    post = [_do_single_on(a0, q0, atol=atol), _do_single_on(a1, q1, atol=atol)]

    return list(cast(Iterable[ops.Operation], ops.flatten_op_tree([
        pre,
        _non_local_part(q0,
                        q1,
                        kak.interaction_coefficients,
                        allow_partial_czs,
                        atol=atol),
        post,
    ])))
1237858
def cublasStbmv(handle, uplo, trans, diag, n, k, A, lda, x, incx):
    
    
    status = _libcublas.cublasStbmv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, k, int(A), lda, int(x), incx)
    cublasCheckStatus(status)
395013
def skew_matrix(w):
    ''''''
    return np.array([[0, -w[2], w[1]],
                     [w[2], 0, -w[0]],
                     [-w[1], w[0], 0]])
163316
def Ak(self, k, m_pred, P_pred):  # returns state iteration matrix
        

        return self.A[:, :, self.index[self.A_time_var_index, k]]
258913
def K(self, parm):
        
        return RQ_K_matrix(self.X, parm) + np.identity(self.X.shape[0])*(10**-10)
786792
function improveRankOrdering(g, fixed, movable) {
    var weights = rankWeights(g, fixed, movable);

    var result = [];

    var layer = movable.slice(0);
    
    // Move fixed nodes into the result array first
    for (var i = 0; i < layer.length; ++i) {
      var u = layer[i];
      if (weights[u] === -1) {
        result[i] = u;
        layer[i] = null;
      }
    }

    layer.sort(function(x, y) { return (x ? weights[x] : -1) - (y ? weights[y] : -1); });

    var nextIdx = 0;
    for (var i = 0; i < layer.length; ++i) {
      if (layer[i] !== null) {
        while (result[nextIdx] !== undefined) {
          ++nextIdx;
        }
        result[nextIdx] = layer[i];
      }
    }

    return result;
  }
906378
def K(self,X,X2,target):
        """"""
        AX = np.dot(X,self.transform)
        if X2 is None:
            X2 = X
            AX2 = AX
        else:
            AX2 = np.dot(X2, self.transform)
        self.k.K(X,X2,target)
        self.k.K(AX,X2,target)
        self.k.K(X,AX2,target)
        self.k.K(AX,AX2,target)
1357713
def cublasCgeru(handle, m, n, alpha, x, incx, y, incy, A, lda):
    

    status = _libcublas.cublasCgeru_v2(handle,
                                       m, n, ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                              alpha.imag)),
                                       int(x), incx, int(y), incy, int(A), lda)
    cublasCheckStatus(status)