Matrix-vector product for complex general banded matrix.
486859
def cublasZgbmv(handle, trans, m, n, kl, ku, alpha, A, lda, 
                x, incx, beta, y, incy):
    

    status = _libcublas.cublasZgbmv_v2(handle,
                                       trans, m, n, kl, ku,
                                       ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                         alpha.imag)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                         beta.imag)),
                              int(y), incy)
    cublasCheckStatus(status)
30541
def _embedding_classical_mds(matrix, dimensions=3, additive_correct=False):
    
    if additive_correct:
        dbc = double_centre(_additive_correct(matrix))
    else:
        dbc = double_centre(matrix)
    decomp = eigen(dbc)
    lambda_ = np.diag(np.sqrt(np.abs(decomp.vals[:dimensions])))
    evecs = decomp.vecs[:, :dimensions]
    coords = evecs.dot(lambda_)
    return coords
443947
def C(w, Xs):
    '''
    '''
    n = len(Xs)
    P = projection_matrix(w)
    Ys = [np.dot(P, X) for X in Xs]
    A = calc_A(Ys)
    A_hat = calc_A_hat(A, skew_matrix(w))

    return np.dot(A_hat, sum(np.dot(Y, Y) * Y for Y in Ys)) / np.trace(np.dot(A_hat, A))
910674
def cublasZgerc(handle, m, n, alpha, x, incx, y, incy, A, lda):
    

    status = _libcublas.cublasZgerc_v2(handle,
                                       m, n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                               alpha.imag)),
                                       int(x), incx, int(y), incy, int(A), lda)
    cublasCheckStatus(status)
521521
def cublasSgbmv(handle, trans, m, n, kl, ku, alpha, A, lda,
                x, incx, beta, y, incy):
    

    status = _libcublas.cublasSgbmv_v2(handle,
                                       trans, m, n, kl, ku,
                                       ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda,
                                       int(x), incx,
                                       ctypes.byref(ctypes.c_float(beta)),
                                       int(y), incy)
    cublasCheckStatus(status)
1028992
def numbaGaussian2d(psf, sy, sx):
    '''
    
    '''
    ps0, ps1 = psf.shape
    c0,c1 = ps0//2, ps1//2
    ssx = 2*sx**2
    ssy = 2*sy**2
    for i in range(ps0):
        for j in range(ps1):
            psf[i,j]=exp( -( (i-c0)**2/ssy
                            +(j-c1)**2/ssx) )
    psf/=psf.sum()
1289444
def cublasSsymm(handle, side, uplo, m, n, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasSsymm_v2(handle,
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       m, n, ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda, int(B), ldb, 
                                       ctypes.byref(ctypes.c_float(beta)), 
                                       int(C), ldc)
    cublasCheckStatus(status)
500618
function svdJs() {
                var A = this;
                var V = Matrix$b.I(A.rows());
                var S = A.transpose();
                var U = Matrix$b.I(A.cols());
                var err = Number.MAX_VALUE;
                var i = 0;
                var maxLoop = 100;

                while(err > 2.2737e-13 && i < maxLoop) {
                    var qr = S.transpose().qrJs();
                    S = qr.R;
                    V = V.x(qr.Q);
                    qr = S.transpose().qrJs();
                    U = U.x(qr.Q);
                    S = qr.R;

                    var e = S.triu(1).unroll().norm();
                    var f = S.diagonal().norm();

                    if(f == 0)
                        f = 1;

                    err = e / f;

                    i++;
                }

                var ss = S.diagonal();
                var s = [];

                for(var i = 1; i <= ss.cols(); i++) {
                    var ssn = ss.e(i);
                    s.push(Math.abs(ssn));

                    if(ssn < 0) {
                        for(var j = 0; j < U.rows(); j++) {
                            V.elements[j][i - 1] = -(V.elements[j][i - 1]);
                        }
                    }
                }

                return {U: U, S: $V(s).toDiagonalMatrix(), V: V};
            }
1299717
public void mul(Complex a, IMatrix result) {
        if ((result.rows != rows) || (result.columns != columns)) result.reshape(rows, columns);

        int i, j;
        for (i = 0; i < rows; i++)
            for (j = 0; j < columns; j++) {
                result.realmatrix[i][j] = realmatrix[i][j] * a.real - imagmatrix[i][j] * a.imag;
                result.imagmatrix[i][j] = realmatrix[i][j] * a.imag + imagmatrix[i][j] * a.real;
            }
    }
611029
def decompose_covariance(c):
    

    # make it a kickass copy of the original
    c = _n.array(c)

    # first get the error vector
    e = []
    for n in range(0, len(c[0])): e.append(_n.sqrt(c[n][n]))

    # now cycle through the matrix, dividing by e[1]*e[2]
    for n in range(0, len(c[0])):
        for m in range(0, len(c[0])):
            c[n][m] = c[n][m] / (e[n]*e[m])

    return [_n.array(e), _n.array(c)]
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
417743
def cublasCsymm(handle, side, uplo, m, n, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasCsymm_v2(handle, 
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       m, n, 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)
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)
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)
6969
def cublasZhpmv(handle, uplo, n, alpha, AP, x, incx, beta, y, incy):
    
    
    status = _libcublas.cublasZhpmv_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                            alpha.imag)),
                                       int(AP), int(x), incx,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real, 
                                                                         beta.imag)),
                                       int(y), incy)
    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;
            }
        }
    }
964528
def cublasZsymv(handle, uplo, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasZsymv_v2(handle, 
                                        _CUBLAS_FILL_MODE[uplo], n, 
                                        ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                            alpha.imag)), 
                                        int(A), lda, int(x), incx, 
                                        ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                            beta.imag)), 
                                        int(y), incy)
    cublasCheckStatus(status)