code: def gemv_c_code(
      y, A, x, z, alpha, beta, fail, force_init_beta=False, params=None
  ):
      """
      z <- beta * y + alpha * dot(A, x)
  
      where A is a matrix, y and x are vectors (ergo z is vector)
      """
      code = """
  
      int elemsize;
      float fbeta;
      double dbeta;
  
      if (PyArray_DIMS(%(A)s)[0] != PyArray_DIMS(%(y)s)[0])
      {
          PyErr_SetString(PyExc_ValueError,
                          "Shape mismatch: A.shape[0] != y.shape[0]");
          %(fail)s;
      }
      if (PyArray_DIMS(%(A)s)[1] != PyArray_DIMS(%(x)s)[0])
      {
          PyErr_SetString(PyExc_ValueError,
                          "Shape mismatch: A.shape[1] != x.shape[0]");
          %(fail)s;
      }
  
      if  (PyArray_DESCR(%(y)s)->type_num == NPY_DOUBLE) { elemsize = 8; }
      else if (PyArray_DESCR(%(y)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 y if not destructive
      if (!%(params)s->inplace)
      {
          if ((NULL == %(z)s)
              || (PyArray_DIMS(%(z)s)[0] != PyArray_DIMS(%(y)s)[0]))
          {
              Py_XDECREF(%(z)s);
              %(z)s = (PyArrayObject*)PyArray_SimpleNew(1,
                  PyArray_DIMS(%(y)s), PyArray_TYPE(%(y)s));
              if(!%(z)s) {
                  PyErr_SetString(PyExc_MemoryError,
                                  "failed to alloc gemv output");
                  %(fail)s
              }
          }
          if (dbeta != 0)
          {
              if (PyArray_CopyInto(%(z)s, %(y)s) != 0) {
                  %(fail)s
              }
          }
          else if (%(force_init_beta)d)
          {
              PyObject *zero = PyFloat_FromDouble(0.);
              if (zero == NULL) %(fail)s;
              if (PyArray_FillWithScalar(%(z)s, zero) != 0) %(fail)s;
              Py_DECREF(zero);
          }
      }
      else
      {
          if (%(z)s != %(y)s)
          {
              Py_XDECREF(%(z)s);
              %(z)s = %(y)s;
              Py_INCREF(%(z)s);
          }
      }
      {
          char TRANS = 'T';
          char NOTRANS = 'N';
          int NA0 = PyArray_DIMS(%(A)s)[0];
          int NA1 = PyArray_DIMS(%(A)s)[1];
          /* This formula is needed in the case where A 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 SA0 = (NA0 > 1) ? (PyArray_STRIDES(%(A)s)[0] / elemsize) : (NA1 + 1);
          int SA1 = (NA1 > 1) ? (PyArray_STRIDES(%(A)s)[1] / elemsize) : (NA0 + 1);
          int Sz = PyArray_STRIDES(%(z)s)[0] / elemsize;
          int Sx = PyArray_STRIDES(%(x)s)[0] / elemsize;
  
          dtype_%(x)s* x_data = (dtype_%(x)s*) PyArray_DATA(%(x)s);
          dtype_%(z)s* z_data = (dtype_%(z)s*) PyArray_DATA(%(z)s);
          // gemv expects pointers to the beginning of memory arrays,
          // but numpy provides a pointer to the first element,
          // so when the stride is negative, we need to get the last one.
          if (Sx < 0)
              x_data += (NA1 - 1) * Sx;
          if (Sz < 0)
              z_data += (NA0 - 1) * Sz;
  
          if (NA0 * NA1)
          {
              // If A 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(%(A)s)[0] < 0)
                  || (PyArray_STRIDES(%(A)s)[1] < 0)
                  || ((PyArray_STRIDES(%(A)s)[0] != elemsize)
                      && (PyArray_STRIDES(%(A)s)[1] != elemsize)))
              {
                  npy_intp dims[2];
                  dims[0] = NA0;
                  dims[1] = NA1;
  
                  PyArrayObject * A_copy = (PyArrayObject *) PyArray_Copy(
                                                                     %(A)s);
                  if (!A_copy)
                      %(fail)s
                  Py_XDECREF(%(A)s);
                  %(A)s = A_copy;
                  SA0 = (NA0 > 1) ? (PyArray_STRIDES(%(A)s)[0] / elemsize) : (NA1 + 1);
                  SA1 = (NA1 > 1) ? (PyArray_STRIDES(%(A)s)[1] / elemsize) : (NA0 + 1);
              }
  
              if (PyArray_STRIDES(%(A)s)[0] == elemsize)
              {
                  if (PyArray_DESCR(%(A)s)->type_num == NPY_FLOAT)
                  {
                      float alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];
                      sgemv_(&NOTRANS, &NA0, &NA1,
                          &alpha,
                          (float*)(PyArray_DATA(%(A)s)), &SA1,
                          (float*)x_data, &Sx,
                          &fbeta,
                          (float*)z_data, &Sz);
                  }
                  else if (PyArray_DESCR(%(A)s)->type_num == NPY_DOUBLE)
                  {
                      double alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];
                      dgemv_(&NOTRANS, &NA0, &NA1,
                          &alpha,
                          (double*)(PyArray_DATA(%(A)s)), &SA1,
                          (double*)x_data, &Sx,
                          &dbeta,
                          (double*)z_data, &Sz);
                  }
                  else
                  {
                      PyErr_SetString(PyExc_AssertionError,
                                      "neither float nor double dtype");
                      %(fail)s
                  }
              }
              else if (PyArray_STRIDES(%(A)s)[1] == elemsize)
              {
                  if (PyArray_DESCR(%(A)s)->type_num == NPY_FLOAT)
                  {
                      float alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];
  
                      // Check for vector-vector dot (NA0 == 1). The code may work
                      // for SA1 != 1 as well, but has not been tested for this case,
                      // so SA1 == 1 is required for safety.
                      if (NA0 == 1 && SA1 == 1)
                      {
                          if (fbeta != 0.f) {
                            z_data[0] = fbeta*z_data[0];
                          } else {
                            z_data[0] = 0.f;
                          }
                          z_data[0] += alpha*sdot_(&NA1,
                                (float*)(PyArray_DATA(%(A)s)), &SA1,
                                (float*)x_data, &Sx);
                      }
                      else
                      {
                          sgemv_(&TRANS, &NA1, &NA0,
                              &alpha,
                              (float*)(PyArray_DATA(%(A)s)), &SA0,
                              (float*)x_data, &Sx,
                              &fbeta,
                              (float*)z_data, &Sz);
                      }
                  }
                  else if (PyArray_DESCR(%(A)s)->type_num == NPY_DOUBLE)
                  {
                      double alpha = ((dtype_%(alpha)s*)PyArray_DATA(%(alpha)s))[0];
  
                      // Check for vector-vector dot (NA0 == 1). The code may work
                      // for SA1 != 1 as well, but has not been tested for this case,
                      // so SA1 == 1 is required for safety.
                      if (NA0 == 1 && SA1 == 1)
                      {
                          if (dbeta != 0.) {
                            z_data[0] = dbeta*z_data[0];
                          } else {
                            z_data[0] = 0.;
                          }
                          z_data[0] += alpha*ddot_(&NA1,
                                (double*)(PyArray_DATA(%(A)s)), &SA1,
                                (double*)x_data, &Sx);
                      }
                      else
                      {
                          dgemv_(&TRANS, &NA1, &NA0,
                              &alpha,
                              (double*)(PyArray_DATA(%(A)s)), &SA0,
                              (double*)x_data, &Sx,
                              &dbeta,
                              (double*)z_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 A.
              dtype_%(z)s * zptr = (dtype_%(z)s*)(PyArray_DATA(%(z)s));
              for (int i = 0; i < NA0; ++i)
              {
                  zptr[i * Sz] = (dbeta == 0.0 ? 0.0 : zptr[i * Sz] * dbeta);
              }
          }
      }
      """
      return code % locals()
  
imports: from pytensor.link.c.op import COp
  from pytensor.link.c.params_type import ParamsType
  from pytensor.scalar import bool
  from pytensor.tensor.blas import Gemv, Ger, blas_header_text, blas_header_version, ldflags
label_desc: Use image node_id 2 for calling a global function with example usage: gemv_c_code(y, A, x, z, alpha, beta, fail, force_init_beta, params) and returns: unknown