def gemv_c_code(aa, xx, yy, zz, alpha, beta, destructive, fail,
                force_init_beta=False):
    code = """

    int elemsize ;
    float fbeta;
    double dbeta;

    if (PyArray_NDIM(%(aa)s) != 1)
    {
        PyErr_SetString(PyExc_NotImplementedError, "Gemv: rank(aa) != 1");
        %(fail)s;
    }
    if (PyArray_NDIM(%(xx)s) != 2)
    {
        PyErr_SetString(PyExc_NotImplementedError, "Gemv: rank(xx) != 2");
        %(fail)s;
    }
    if (PyArray_NDIM(%(yy)s) != 1)
    {
        PyErr_SetString(PyExc_NotImplementedError, "Gemv: rank(yy) != 1");
        %(fail)s;
    }
    if (PyArray_NDIM(%(alpha)s) != 0)
    {
        PyErr_SetString(PyExc_NotImplementedError, "Gemv: rank(alpha) != 0");
        %(fail)s;
    }
    if (PyArray_NDIM(%(beta)s) != 0)
    {
        PyErr_SetString(PyExc_NotImplementedError, "Gemv: rank(beta) != 0");
        %(fail)s;
    }

    if (PyArray_DESCR(%(aa)s)->type_num != PyArray_DESCR(%(xx)s)->type_num)
    { PyErr_SetString(PyExc_TypeError, "Gemv: aa vs. xx"); %(fail)s; }
    if (PyArray_DESCR(%(aa)s)->type_num != PyArray_DESCR(%(yy)s)->type_num)
    { PyErr_SetString(PyExc_TypeError, "Gemv: aa vs. yy"); %(fail)s; }

    if (PyArray_DIMS(%(xx)s)[0] != PyArray_DIMS(%(aa)s)[0])
    {
        PyErr_SetString(PyExc_ValueError,
                        "Shape mismatch: A.shape[0] != x.shape[0]");
        %(fail)s;
    }
    if (PyArray_DIMS(%(xx)s)[1] != PyArray_DIMS(%(yy)s)[0])
    {
        PyErr_SetString(PyExc_ValueError,
                        "Shape mismatch: A.shape[1] != y.shape[0]");
        %(fail)s;
    }

    if  (PyArray_DESCR(%(aa)s)->type_num == NPY_DOUBLE) { elemsize = 8; }
    else if (PyArray_DESCR(%(aa)s)->type_num == NPY_FLOAT) { elemsize = 4;}
    else {
        PyErr_SetString(PyExc_NotImplementedError, "complex Gemv");
        %(fail)s;
    }

    fbeta = dbeta = ((dtype_%(beta)s*)PyArray_DATA(%(beta)s))[0];

    // copy aa if not destructive
    if (!%(destructive)s)
    {
        if ((NULL == %(zz)s)
            || (PyArray_DIMS(%(zz)s)[0] != PyArray_DIMS(%(aa)s)[0]))
        {
            Py_XDECREF(%(zz)s);
            %(zz)s = (PyArrayObject*)PyArray_SimpleNew(1,
                PyArray_DIMS(%(aa)s), PyArray_TYPE(%(aa)s));
            if(!%(zz)s) {
                PyErr_SetString(PyExc_MemoryError,
                                "failed to alloc gemv output");
                %(fail)s
            }
        }
        if (%(zz)s == %(aa)s)
        {
            PyErr_SetString(PyExc_AssertionError, "%(zz)s != %(aa)s");
            %(fail)s
        }
        if (dbeta != 0)
        {
            if (PyArray_DESCR(%(zz)s)->type_num == NPY_FLOAT)
            {
                float * zoutdata = (float*)PyArray_DATA(%(zz)s);
                const float * zdata = (float*)PyArray_DATA(%(aa)s);
                int Ai = PyArray_STRIDES(%(aa)s)[0]/sizeof(float);
                int Zi = PyArray_STRIDES(%(zz)s)[0]/sizeof(float);
                for (int i = 0; i < PyArray_DIMS(%(aa)s)[0]; ++i)
                {
                    zoutdata[Zi*i] = fbeta * zdata[Ai*i];
                }
            }
            else if (PyArray_DESCR(%(zz)s)->type_num == NPY_DOUBLE)
            {
                double * zoutdata = (double*) PyArray_DATA(%(zz)s);
                const double * zdata = (double*)PyArray_DATA(%(aa)s);
                int Ai = PyArray_STRIDES(%(aa)s)[0]/sizeof(double);
                int Zi = PyArray_STRIDES(%(zz)s)[0]/sizeof(double);
                for (int i = 0; i < PyArray_DIMS(%(aa)s)[0]; ++i)
                {
                    zoutdata[Zi*i] = dbeta * zdata[Ai*i];
                }
            }
            else
            {
                PyErr_SetString(PyExc_AssertionError,
                                "neither float nor double dtype");
                %(fail)s
            }
            fbeta = dbeta = 1.0;
        }
        else if (%(force_init_beta)d)
        {
            if (PyArray_CHKFLAGS(%(zz)s, NPY_ARRAY_C_CONTIGUOUS))
            {
                memset((void *)PyArray_DATA(%(zz)s), 0, PyArray_SIZE(%(zz)s)*PyArray_ITEMSIZE(%(zz)s));
            }
            else
            {
                if (PyArray_DESCR(%(zz)s)->type_num == NPY_FLOAT)
                {
                    float *zoutdata = (float *)PyArray_DATA(%(zz)s);
                    int Zi = PyArray_STRIDES(%(zz)s)[0]/sizeof(float);
                    for (int i = 0; i < PyArray_DIMS(%(aa)s)[0]; ++i)
                    {
                        zoutdata[Zi*i] = 0.0f;
                    }
                }
                else if (PyArray_DESCR(%(zz)s)->type_num == NPY_DOUBLE)
                {
                    double *zoutdata = (double *)PyArray_DATA(%(zz)s);
                    int Zi = PyArray_STRIDES(%(zz)s)[0]/sizeof(double);
                    for (int i = 0; i < PyArray_DIMS(%(aa)s)[0]; ++i)
                    {
                        zoutdata[Zi*i] = 0.0;
                    }
                }
                else
                {
                    PyErr_SetString(PyExc_AssertionError,
                                    "neither float nor double dtype");
                    %(fail)s
                }
            }
        }
    }
    else
    {
        //fprintf(stderr, "Gemv working in-place \\n");
        if (%(zz)s != %(aa)s)
        {
            if (%(zz)s) { Py_DECREF(%(zz)s); }
            %(zz)s = %(aa)s;
            Py_INCREF(%(zz)s);
        }
    }
    {
        char TRANS = 'T';
        char NOTRANS = 'N';
        int Nx0 = PyArray_DIMS(%(xx)s)[0];
        int Nx1 = PyArray_DIMS(%(xx)s)[1];
        /* This formula is needed in the case where xx is actually a row or
         * column matrix, because BLAS sometimes insists that the strides:
         *  - are not smaller than the number of elements in the array
         *  - are not 0.
         */
        int Sx0 = (Nx0 > 1) ? (PyArray_STRIDES(%(xx)s)[0] / elemsize) : (Nx1 + 1);
        int Sx1 = (Nx1 > 1) ? (PyArray_STRIDES(%(xx)s)[1] / elemsize) : (Nx0 + 1);
        int Sz = PyArray_STRIDES(%(zz)s)[0] / elemsize;
        int Sy = PyArray_STRIDES(%(yy)s)[0] / elemsize;

        dtype_%(yy)s* yy_data = (dtype_%(yy)s*) PyArray_DATA(%(yy)s);
        dtype_%(zz)s* zz_data = (dtype_%(zz)s*) PyArray_DATA(%(zz)s);
        // gemv expects pointers to the beginning of memory arrays,
        // but numpy provides provides a pointer to the first element,
        // so when the stride is negative, we need to get the last one.
        if (Sy < 0)
            yy_data += (Nx1 - 1) * Sy;
        if (Sz < 0)
            zz_data += (Nx0 - 1) * Sz;

        if (Nx0 * Nx1)
        {
            // If xx is neither C- nor F-contiguous, we make a copy.
            // TODO:
            // - if one stride is equal to "- elemsize", we can still call
            //   gemv on reversed matrix and vectors
            // - if the copy is too long, maybe call vector/vector dot on
            //   each row instead
            if ((PyArray_STRIDES(%(xx)s)[0] < 0)
                || (PyArray_STRIDES(%(xx)s)[1] < 0)
                || ((PyArray_STRIDES(%(xx)s)[0] != elemsize)
                    && (PyArray_STRIDES(%(xx)s)[1] != elemsize)))
            {
                npy_intp dims[2];
                dims[0] = Nx0;
                dims[1] = Nx1;

                PyArrayObject * xx_copy = (PyArrayObject *) PyArray_Copy(
                                                                    %(xx)s);
                if (!xx_copy)
                    %(fail)s
                Py_XDECREF(%(xx)s);
                %(xx)s = xx_copy;
                Sx0 = (Nx0 > 1) ? (PyArray_STRIDES(%(xx)s)[0] / elemsize) : (Nx1 + 1);
                Sx1 = (Nx1 > 1) ? (PyArray_STRIDES(%(xx)s)[1] / elemsize) : (Nx0 + 1);
            }

            if (PyArray_STRIDES(%(xx)s)[0] == elemsize)
            {
                if (PyArray_DESCR(%(xx)s)->type_num == NPY_FLOAT)
                {
                    //fprintf(stderr, "A\\n");
                    float alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];
                    sgemv_(&NOTRANS, &Nx0, &Nx1,
                        &alpha,
                        (float*)(PyArray_DATA(%(xx)s)), &Sx1,
                        (float*)yy_data, &Sy,
                        &fbeta,
                        (float*)zz_data, &Sz);
                }
                else if (PyArray_DESCR(%(xx)s)->type_num == NPY_DOUBLE)
                {
                    double alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];
                    dgemv_(&NOTRANS, &Nx0, &Nx1,
                        &alpha,
                        (double*)(PyArray_DATA(%(xx)s)), &Sx1,
                        (double*)yy_data, &Sy,
                        &dbeta,
                        (double*)zz_data, &Sz);
                }
                else
                {
                    PyErr_SetString(PyExc_AssertionError,
                                    "neither float nor double dtype");
                    %(fail)s
                }
            }
            else if (PyArray_STRIDES(%(xx)s)[1] == elemsize)
            {
                if (PyArray_DESCR(%(xx)s)->type_num == NPY_FLOAT)
                {
                    float alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];

                    // Check for vector-vector dot (Nx0 == 1). The code may work
                    // for Sx1 != 1 as well, but has not been tested for this case,
                    // so Sx1 == 1 is required for safety.
                    if (Nx0 == 1 && Sx1 == 1)
                    {
                        zz_data[0] = fbeta*zz_data[0] + alpha*sdot_(&Nx1,
                            (float*)(PyArray_DATA(%(xx)s)), &Sx1,
                            (float*)yy_data, &Sy);
                    }
                    else
                    {
                        sgemv_(&TRANS, &Nx1, &Nx0,
                            &alpha,
                            (float*)(PyArray_DATA(%(xx)s)), &Sx0,
                            (float*)yy_data, &Sy,
                            &fbeta,
                            (float*)zz_data, &Sz);
                    }
                }
                else if (PyArray_DESCR(%(xx)s)->type_num == NPY_DOUBLE)
                {
                    double alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];

                    // Check for vector-vector dot (Nx0 == 1). The code may work
                    // for Sx1 != 1 as well, but has not been tested for this case,
                    // so Sx1 == 1 is required for safety.
                    if (Nx0 == 1 && Sx1 == 1)
                    {
                        zz_data[0] = dbeta*zz_data[0] + alpha*ddot_(&Nx1,
                              (double*)(PyArray_DATA(%(xx)s)), &Sx1,
                              (double*)yy_data, &Sy);
                    }
                    else
                    {
                        dgemv_(&TRANS, &Nx1, &Nx0,
                            &alpha,
                            (double*)(PyArray_DATA(%(xx)s)), &Sx0,
                            (double*)yy_data, &Sy,
                            &dbeta,
                            (double*)zz_data, &Sz);
                    }
                }
                else
                {
                    PyErr_SetString(PyExc_AssertionError,
                                    "neither float nor double dtype");
                    %(fail)s
                }
            }
            else
            {
                PyErr_SetString(PyExc_AssertionError,
                    "xx is a double-strided matrix, and should have been "
                    "copied into a memory-contiguous one.");
                %(fail)s
            }
        }
        else if (dbeta != 1.0)
        {
            // the matrix has at least one dim of length 0
            // so we do this loop, which either iterates over 0 elements
            // or else it does the right thing for length-0 x.
            dtype_%(zz)s * zptr = (dtype_%(zz)s*)(PyArray_DATA(%(zz)s));
            for (int i = 0; i < Nx0; ++i)
            {
                zptr[i * Sz] *= dbeta;
            }
        }
    }

    """
    return code % locals()
Get the iterator flag for the given method.
dataset/ETHPy150Open rizar/attention-lvcsr/libs/Theano/theano/tensor/blas_c.py/gemv_c_code--crossover with--dataset/ETHPy150Open cloudmatrix/esky/esky/sudo/__init__.py/_get_sudo_iterator