Matrix-diagonal matrix product for real general matrix.
341621
def cublasSdgmm(handle, mode, m, n, A, lda, x, incx, C, ldc):
    

    status = _libcublas.cublasSdgmm(handle,
                                    _CUBLAS_SIDE[mode],
                                    m, n, 
                                    int(A), lda, 
                                    int(x), incx,
                                    int(C), ldc)
    cublasCheckStatus(status)
320729
function svdPack() {
                var result = lapack$1.sgesvd('A', 'A', this.elements);

                return {
                    U: $M(result.U),
                    S: $M(result.S).column(1).toDiagonalMatrix(),
            	V: $M(result.VT).transpose()
                };
            }
1337662
public static Basic2DMatrix diagonal(int size, double diagonal) {
        double[][] array = new double[size][size];

        for (int i = 0; i < size; i++) {
            array[i][i] = diagonal;
        }

        return new Basic2DMatrix(array);
    }
32165
def xstep(self):
        """

        self.YU[:] = self.Y - self.U
        YUf = sl.rfftn(self.YU, None, self.cri.axisN)
        YUf0 = self.block_sep0(YUf)
        YUf1 = self.block_sep1(YUf)

        b = self.rho * np.sum(np.conj(self.GDf) * YUf1, axis=-1)
        if self.cri.Cd > 1:
            b = np.sum(b, axis=self.cri.axisC, keepdims=True)
        b += self.DSf + self.rho*YUf0

        # Concatenate multiple GDf components on axisC. For
        # single-channel signals, and multi-channel signals with a
        # single-channel dictionary, we end up with sl.solvemdbi_ism
        # solving a linear system of rank dimN+1 (corresponding to the
        # dictionary and a gradient operator per spatial dimension) plus
        # an identity. For multi-channel signals with a multi-channel
        # dictionary, we end up with sl.solvemdbi_ism solving a linear
        # system of rank C.d (dimN+1) (corresponding to the dictionary
        # and a gradient operator per spatial dimension for each
        # channel) plus an identity.

        # The structure of the linear system to be solved depends on the
        # number of channels in the signal and dictionary. Both branches are
        # the same in the single-channel signal case (the choice of handling
        # it via the 'else' branch is somewhat arbitrary).
        if self.cri.C > 1 and self.cri.Cd == 1:
            # Concatenate multiple GDf components on the final axis
            # of GDf (that indexes the number of gradient operators). For
            # multi-channel signals with a single-channel dictionary,
            # sl.solvemdbi_ism has to solve a linear system of rank dimN+1
            # (corresponding to the dictionary and a gradient operator per
            # spatial dimension)
            DfGDf = np.concatenate(
                [self.Df[..., np.newaxis],] +
                [np.sqrt(self.rho)*self.GDf[..., k, np.newaxis] for k
                 in range(self.GDf.shape[-1])], axis=-1)
            self.Xf[:] = sl.solvemdbi_ism(DfGDf, self.rho, b[..., np.newaxis],
                                          self.cri.axisM, -1)[..., 0]
        else:
            # Concatenate multiple GDf components on axisC. For multi-channel
            # signals with a multi-channel dictionary, sl.solvemdbi_ism has
            # to solve a linear system of rank C.d (dimN+1) (corresponding to
            # the dictionary and a gradient operator per spatial dimension
            # for each channel) plus an identity.
            DfGDf = np.concatenate(
                [self.Df,] + [np.sqrt(self.rho)*self.GDf[..., k] for k
                              in range(self.GDf.shape[-1])],
                axis=self.cri.axisC)
            self.Xf[:] = sl.solvemdbi_ism(DfGDf, self.rho, b, self.cri.axisM,
                                          self.cri.axisC)

        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            if self.cri.C > 1 and self.cri.Cd == 1:
                Dop = lambda x: sl.inner(DfGDf, x[..., np.newaxis],
                                         axis=self.cri.axisM)
                DHop = lambda x: sl.inner(np.conj(DfGDf), x, axis=-1)
                ax = DHop(Dop(self.Xf))[..., 0] + self.rho*self.Xf
            else:
                Dop = lambda x: sl.inner(DfGDf, x, axis=self.cri.axisM)
                DHop = lambda x: sl.inner(np.conj(DfGDf), x,
                                          axis=self.cri.axisC)
                ax = DHop(Dop(self.Xf)) + self.rho*self.Xf
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
1174972
def rotation_matrix_2D(th):
    '''
    
    '''
    s = np.sin(th)
    c = np.cos(th)
    return np.array([[c, -s], [s, c]])
543285
def Svd(a, uv, full):
    
    u, s, v = np.linalg.svd(a, full_matrices=full, compute_uv=uv)
    return s, u, v
1300178
def diagonalize(self):
        ''''''
        self.eigvals, self.eigvecs = np.linalg.eig(
            (self.tensor.transpose() + self.tensor) / 2.0)
        self.eigvals = np.diag(np.dot(
            np.dot(self.eigvecs.transpose(), self.tensor), self.eigvecs))
834516
def MatrixInverse(a, adj):
    
    return np.linalg.inv(a if not adj else _adjoint(a)),
341185
def cublasZgemm(handle, transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc):
    

    status = _libcublas.cublasZgemm_v2(handle,
                                       _CUBLAS_OP[transa],
                                       _CUBLAS_OP[transb], m, n, k, 
                                       ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                         alpha.imag)),
                                       int(A), lda, int(B), ldb,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                         beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
1262849
public void householderSymmetric( int row , double gamma )
    {
        int startU = row*N;

        // compute v = -gamma*A*u
        for( int i = row+1; i < N; i++ ) {
            double totalReal = 0;
            double totalImag = 0;

            // the lower triangle is not written to so it needs to traverse upwards
            // to get the information.  Reduces the number of matrix writes need
            // improving large matrix performance
            for( int j = row+1; j < i; j++ ) {
                double realA = QT.data[(j*N+i)*2];
                double imagA = -QT.data[(j*N+i)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            for( int j = i; j < N; j++ ) {
                double realA = QT.data[(i*N+j)*2];
                double imagA = QT.data[(i*N+j)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            w[i*2]   = -gamma*totalReal;
            w[i*2+1] = -gamma*totalImag;
        }

        // alpha = -0.5*gamma*u^T*v
        double realAplha = 0;
        double imageAlpha = 0;

        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = -QT.data[(startU+i)*2+1];

            double realV = w[i*2];
            double imagV = w[i*2+1];

            realAplha += realU*realV - imagU*imagV;
            imageAlpha += realU*imagV + imagU*realV;
        }
        realAplha *= -0.5*gamma;
        imageAlpha *= -0.5*gamma;

        // w = v + alpha*u
        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = QT.data[(startU+i)*2+1];

            w[i*2]   += realAplha*realU - imageAlpha*imagU;
            w[i*2+1] += realAplha*imagU + imageAlpha*realU;
        }

        // A = A + w*u^T + u*w^T
        for( int i = row+1; i < N; i++ ) {

            double realWW = w[i*2];
            double imagWW = w[i*2+1];

            double realUU = QT.data[(startU+i)*2];
            double imagUU = QT.data[(startU+i)*2+1];

            int indA = (i*N+i)*2;
            for( int j = i; j < N; j++ ) {
                // only write to the upper portion of the matrix
                // this reduces the number of cache misses

                double realU = QT.data[(startU+j)*2];
                double imagU = -QT.data[(startU+j)*2+1];

                double realW = w[j*2];
                double imagW = -w[j*2+1];

                QT.data[indA++] += realWW*realU - imagWW*imagU + realW*realUU - imagW*imagUU;
                QT.data[indA++] += realWW*imagU + imagWW*realU + realW*imagUU + imagW*realUU;
            }
        }
    }
204942
def _svd_sym_koopman(K, C00_train, Ctt_train):
    
    from pyemma._ext.variational.solvers.direct import spd_inv_sqrt
    # reweight operator to empirical distribution
    C0t_re = mdot(C00_train, K)
    # symmetrized operator and SVD
    K_sym = mdot(spd_inv_sqrt(C00_train), C0t_re, spd_inv_sqrt(Ctt_train))
    U, S, Vt = np.linalg.svd(K_sym, compute_uv=True, full_matrices=False)
    # projects back to singular functions of K
    U = mdot(spd_inv_sqrt(C00_train), U)
    Vt = mdot(Vt,spd_inv_sqrt(Ctt_train))
    return U, S, Vt.T
1214304
def grad_eigh(ans, x, UPLO='L'):
    """"""
    N = x.shape[-1]
    w, v = ans              # Eigenvalues, eigenvectors.
    def vjp(g):
        wg, vg = g          # Gradient w.r.t. eigenvalues, eigenvectors.
        w_repeated = anp.repeat(w[..., anp.newaxis], N, axis=-1)
        off_diag = anp.ones((N, N)) - anp.eye(N)
        F = off_diag / (T(w_repeated) - w_repeated + anp.eye(N))
        return _dot(v * wg[..., anp.newaxis, :] + _dot(v, F * _dot(T(v), vg)), T(v))
    return vjp
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);
  }
}
613024
def _fill_diagonals(m, diag_indices):
    """"""
    assert m.ndim == 3, "M must have 3 dimensions"
    assert m.shape[1] == m.shape[2], "M must contain square matrices"
    for r in range(m.shape[0]):
        scipy.fill_diagonal(m[r], 0)
        m[r][diag_indices] -= scipy.sum(m[r], axis=1)
1320237
function multiply_matrices(l, r) {
          var al = l[0];
          var cl = l[1];
          var el = l[2];
          var bl = l[3];
          var dl = l[4];
          var fl = l[5];

          var ar = r[0];
          var cr = r[1];
          var er = r[2];
          var br = r[3];
          var dr = r[4];
          var fr = r[5];

          var a = al * ar + cl * br;
          var c = al * cr + cl * dr;
          var e = al * er + cl * fr + el;
          var b = bl * ar + dl * br;
          var d = bl * cr + dl * dr;
          var f = bl * er + dl * fr + fl;

          return [a, c, e, b, d, f];
      }
281035
def mag_pmf(matrix):
    
    nmags, ndists, nlons, nlats, neps = matrix.shape
    mag_pmf = numpy.zeros(nmags)
    for i in range(nmags):
        mag_pmf[i] = numpy.prod(
            [1. - matrix[i, j, k, l, m]
             for j in range(ndists)
             for k in range(nlons)
             for l in range(nlats)
             for m in range(neps)])
    return 1. - mag_pmf
4449
private DoubleMatrix1D rDual(DoubleMatrix1D X, DoubleMatrix1D V, DoubleMatrix1D gradX) {
		if(getA()==null){
			return gradX;
		}
		//return getAT().zMult(V, gradX.copy(), 1., 1, false);
		return ColtUtils.zMult(getAT(), V, gradX, 1);
	}