PyObject* multi_axpy(PyObject *self, PyObject *args)
{
  PyArrayObject* alpha;
  PyArrayObject* x;
  PyArrayObject* y;
  if (!PyArg_ParseTuple(args, "OOO", &alpha, &x, &y)) 
    return NULL;
  int n0 = PyArray_DIMS(x)[0];
  int n = PyArray_DIMS(x)[1];
  for (int d = 2; d < PyArray_NDIM(x); d++)
    n *= PyArray_DIMS(x)[d];
  int incx = 1;
  int incy = 1;

   if (PyArray_DESCR(alpha)->type_num == NPY_DOUBLE)
    {
      if (PyArray_DESCR(x)->type_num == NPY_CDOUBLE)
        n *= 2;
      double *ap = DOUBLEP(alpha);
      double *xp = DOUBLEP(x);
      double *yp = DOUBLEP(y);
      for (int i = 0; i < n0; i++)
        {
          daxpy_(&n, &ap[i], 
                 (void*)xp, &incx,
                 (void*)yp, &incy);
          xp += n;
          yp += n;
        }
    }
  else
    {
      double_complex *ap = COMPLEXP(alpha);
      double_complex *xp = COMPLEXP(x);
      double_complex *yp = COMPLEXP(y);
      for (int i = 0; i < n0; i++)
        {
          zaxpy_(&n, (void*)(&ap[i]), 
                 (void*)xp, &incx,
                 (void*)yp, &incy);
          xp += n;
          yp += n;
        }
    }
  Py_RETURN_NONE;
}

int NCFormat_from_spMatrix(SuperMatrix * A, int m, int n, int nnz,
			   PyArrayObject * nzvals, PyArrayObject * rowind,
			   PyArrayObject * colptr, int typenum)
{
    int ok = 0;

    ok = (PyArray_EquivTypenums(PyArray_DESCR(nzvals)->type_num, typenum) &&
          PyArray_EquivTypenums(PyArray_DESCR(rowind)->type_num, NPY_INT) &&
          PyArray_EquivTypenums(PyArray_DESCR(colptr)->type_num, NPY_INT) &&
          PyArray_NDIM(nzvals) == 1 &&
          PyArray_NDIM(rowind) == 1 &&
          PyArray_NDIM(colptr) == 1 &&
          PyArray_IS_C_CONTIGUOUS(nzvals) &&
          PyArray_IS_C_CONTIGUOUS(rowind) &&
          PyArray_IS_C_CONTIGUOUS(colptr) &&
          nnz <= PyArray_DIM(nzvals, 0) &&
          nnz <= PyArray_DIM(rowind, 0) &&
          n+1 <= PyArray_DIM(colptr, 0));
    if (!ok) {
	PyErr_SetString(PyExc_ValueError,
			"sparse matrix arrays must be 1-D C-contiguous and of proper "
                        "sizes and types");
	return -1;
    }


    if (setjmp(_superlu_py_jmpbuf))
	return -1;
    else {
	if (!CHECK_SLU_TYPE(nzvals->descr->type_num)) {
	    PyErr_SetString(PyExc_TypeError, "Invalid type for array.");
	    return -1;
	}
	Create_CompCol_Matrix(nzvals->descr->type_num,
			      A, m, n, nnz, nzvals->data,
			      (int *) rowind->data, (int *) colptr->data,
			      SLU_NC,
			      NPY_TYPECODE_TO_SLU(nzvals->descr->type_num),
			      SLU_GE);
    }

    return 0;
}

/* evaluates if given array is non-empty 1D NumPy with double elements */
int is1DDoubleArray(PyArrayObject* x)
{
    if (!(PyArray_NDIM(x) == 1 && PyArray_DIM(x, 0) > 0 &&
          PyArray_DESCR(x)->type_num == NPY_DOUBLE)) {
        PyErr_SetString(PyExc_ValueError,
                        "Array(s) must be non-empty 1D array of type Float");
        return 0;
    } else
        return 1;
}

NPY_NO_EXPORT PyObject *
array_int(PyArrayObject *v)
{
    PyObject *pv, *pv2;
    if (PyArray_SIZE(v) != 1) {
        PyErr_SetString(PyExc_TypeError, "only length-1 arrays can be"\
                        " converted to Python scalars");
        return NULL;
    }
    pv = PyArray_DESCR(v)->f->getitem(PyArray_DATA(v), v);
    if (pv == NULL) {
        return NULL;
    }
    if (Py_TYPE(pv)->tp_as_number == 0) {
        PyErr_SetString(PyExc_TypeError, "cannot convert to an int; "\
                        "scalar object is not a number");
        Py_DECREF(pv);
        return NULL;
    }
    if (Py_TYPE(pv)->tp_as_number->nb_int == 0) {
        PyErr_SetString(PyExc_TypeError, "don't know how to convert "\
                        "scalar number to int");
        Py_DECREF(pv);
        return NULL;
    }
    /*
     * If we still got an array which can hold references, stop
     * because it could point back at 'v'.
     */
    if (PyArray_Check(pv) &&
                PyDataType_REFCHK(PyArray_DESCR((PyArrayObject *)pv))) {
        PyErr_SetString(PyExc_TypeError,
                "object array may be self-referencing");
        Py_DECREF(pv);
        return NULL;
    }

    pv2 = Py_TYPE(pv)->tp_as_number->nb_int(pv);
    Py_DECREF(pv);
    return pv2;
}

static PyObject *
array_repr_builtin(PyArrayObject *self, int repr)
{
    PyObject *ret;
    char *string;
    /* max_n initial value is arbitrary, dump_data will extend it */
    Py_ssize_t n = 0, max_n = PyArray_NBYTES(self) * 4 + 7;

    if ((string = PyArray_malloc(max_n)) == NULL) {
        return PyErr_NoMemory();
    }

    if (dump_data(&string, &n, &max_n, PyArray_DATA(self),
                  PyArray_NDIM(self), PyArray_DIMS(self),
                  PyArray_STRIDES(self), self) < 0) {
        PyArray_free(string);
        return NULL;
    }

    if (repr) {
        if (PyArray_ISEXTENDED(self)) {
            ret = PyUString_FromFormat("array(%s, '%c%d')",
                                       string,
                                       PyArray_DESCR(self)->type,
                                       PyArray_DESCR(self)->elsize);
        }
        else {
            ret = PyUString_FromFormat("array(%s, '%c')",
                                       string,
                                       PyArray_DESCR(self)->type);
        }
    }
    else {
        ret = PyUString_FromStringAndSize(string, n);
    }

    PyArray_free(string);
    return ret;
}

//initialization of BallTree object
// argument is a single array of size [D,N]
static int
BallTree_init(BallTreeObject *self, PyObject *args, PyObject *kwds){
  //we use goto statements : all variables should be declared up front
    PyObject *arg=NULL;
    PyObject *arr=NULL;
    long int leaf_size=20;

    static char *kwlist[] = {"x", "leafsize", NULL};

    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|l", kwlist,
        &arg,&leaf_size) )
        goto fail;

    if(leaf_size <= 0){
        PyErr_SetString(PyExc_ValueError,
            "BallTree : leaf size must be greater than zero");
        goto fail;
    }

  //view this object as an array of doubles,
    arr = PyArray_FROM_OTF(arg,NPY_DOUBLE,0);
    if(arr==NULL)
        goto fail;

  //check that it is 2D
    if( PyArray_NDIM(arr) != 2)
        goto fail;

    if(self != NULL){
    //create the list of points
        self->size = PyArray_DIMS(arr)[0];
        self->dim = PyArray_DIMS(arr)[1];

        int inc = PyArray_STRIDES(arr)[1]/PyArray_DESCR(arr)->elsize;
        self->Points = new std::vector<BallTree_Point*>(self->size);
        for(int i=0;i<self->size;i++)
            self->Points->at(i) = new BallTree_Point(arr,
            (double*)PyArray_GETPTR2(arr,i,0),
            inc, PyArray_DIM(arr,1));
        self->tree = new BallTree<BallTree_Point>(*(self->Points),
            leaf_size);
    }
    self->data = arr;
  //Py_DECREF(arr);
    return 0;

    fail:
    Py_XDECREF(arr);
    return -1;
}

/*
 * check if in "alhs @[email protected] orhs" that alhs is a temporary (refcnt == 1) so we
 * can do in-place operations instead of creating a new temporary
 * "cannot" is set to true if it cannot be done even with swapped arguments
 */
static int
can_elide_temp(PyArrayObject * alhs, PyObject * orhs, int * cannot)
{
    /*
     * to be a candidate the array needs to have reference count 1, be an exact
     * array of a basic type, own its data and size larger than threshold
     */
    if (Py_REFCNT(alhs) != 1 || !PyArray_CheckExact(alhs) ||
            !PyArray_ISNUMBER(alhs) ||
            !PyArray_CHKFLAGS(alhs, NPY_ARRAY_OWNDATA) ||
            !PyArray_ISWRITEABLE(alhs) ||
            PyArray_CHKFLAGS(alhs, NPY_ARRAY_UPDATEIFCOPY) ||
            PyArray_CHKFLAGS(alhs, NPY_ARRAY_WRITEBACKIFCOPY) ||
            PyArray_NBYTES(alhs) < NPY_MIN_ELIDE_BYTES) {
        return 0;
    }
    if (PyArray_CheckExact(orhs) ||
        PyArray_CheckAnyScalar(orhs)) {
        PyArrayObject * arhs;

        /* create array from right hand side */
        Py_INCREF(orhs);
        arhs = (PyArrayObject *)PyArray_EnsureArray(orhs);
        if (arhs == NULL) {
            return 0;
        }

        /*
         * if rhs is not a scalar dimensions must match
         * TODO: one could allow broadcasting on equal types
         */
        if (!(PyArray_NDIM(arhs) == 0 ||
              (PyArray_NDIM(arhs) == PyArray_NDIM(alhs) &&
               PyArray_CompareLists(PyArray_DIMS(alhs), PyArray_DIMS(arhs),
                                    PyArray_NDIM(arhs))))) {
                Py_DECREF(arhs);
                return 0;
        }

        /* must be safe to cast (checks values for scalar in rhs) */
        if (PyArray_CanCastArrayTo(arhs, PyArray_DESCR(alhs),
                                   NPY_SAFE_CASTING)) {
            Py_DECREF(arhs);
            return check_callers(cannot);
        }
        Py_DECREF(arhs);
    }

    return 0;
}

static PyObject * merge(PyObject * self,
     PyObject * args, PyObject * kwds) {

    static char * kwlist[] = {
        "data", "A", "B", "out", NULL
    };

    PyArrayObject * data, * A, * B, * out;
    if(!PyArg_ParseTupleAndKeywords(args, kwds, 
        "O!O!O!O!", kwlist,
        &PyArray_Type, &data, 
        &PyArray_Type, &A, 
        &PyArray_Type, &B, 
        &PyArray_Type, &out)) return NULL;

    npy_intp i = 0;

    int (*compare)(void *, void *, void*) = PyArray_DESCR(data)->f->compare;

    size_t sizeA = PyArray_SIZE(A);
    size_t sizeB = PyArray_SIZE(B);
    npy_intp * Aptr = PyArray_DATA(A);
    npy_intp * Aend = Aptr + sizeA;

    npy_intp * Bptr = PyArray_DATA(B);
    npy_intp * Bend = Bptr + sizeB;
    npy_intp * Optr = PyArray_DATA(out);

    #define VA PyArray_GETPTR1(data, *Aptr)
    #define VB PyArray_GETPTR1(data, (*Bptr) + sizeA)

    Py_BEGIN_ALLOW_THREADS
    while(Aptr < Aend|| Bptr < Bend) {
        while(Aptr < Aend && (Bptr == Bend || compare(VA, VB, data) <= 0)) {
            *Optr = *Aptr;
            Aptr++;
            Optr++;
            //printf("adding from A, i = %ld, k = %ld v = %ld\n", i, k, v);
        }
        while(Bptr < Bend && (Aptr == Aend || compare(VA, VB, data) >= 0)) {
            *Optr = *Bptr + sizeA;
            Bptr++;
            Optr++;
            //printf("adding from B, j = %ld, k = %ld v = %ld\n", j, k, v);
        }
    }
    Py_END_ALLOW_THREADS
    Py_INCREF(Py_None);
    return Py_None;
}

PyObject* gemm(PyObject *self, PyObject *args)
{
  Py_complex alpha;
  PyArrayObject* a;
  PyArrayObject* b;
  Py_complex beta;
  PyArrayObject* c;
  char transa = 'n';
  if (!PyArg_ParseTuple(args, "DOODO|c", &alpha, &a, &b, &beta, &c, &transa))
    return NULL;
  int m, k, lda, ldb, ldc;
  if (transa == 'n')
    {
      m = PyArray_DIMS(a)[1];
      for (int i = 2; i < PyArray_NDIM(a); i++)
        m *= PyArray_DIMS(a)[i];
      k = PyArray_DIMS(a)[0];
      lda = MAX(1, PyArray_STRIDES(a)[0] / PyArray_STRIDES(a)[PyArray_NDIM(a) - 1]);
      ldb = MAX(1, PyArray_STRIDES(b)[0] / PyArray_STRIDES(b)[1]);
      ldc = MAX(1, PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[PyArray_NDIM(c) - 1]);
    }
  else
    {
      k = PyArray_DIMS(a)[1];
      for (int i = 2; i < PyArray_NDIM(a); i++)
        k *= PyArray_DIMS(a)[i];
      m = PyArray_DIMS(a)[0];
      lda = MAX(1, k);
      ldb = MAX(1, PyArray_STRIDES(b)[0] / PyArray_STRIDES(b)[PyArray_NDIM(b) - 1]);
      ldc = MAX(1, PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[1]);

    }
  int n = PyArray_DIMS(b)[0];
  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    dgemm_(&transa, "n", &m, &n, &k,
           &(alpha.real),
           DOUBLEP(a), &lda,
           DOUBLEP(b), &ldb,
           &(beta.real),
           DOUBLEP(c), &ldc);
  else
    zgemm_(&transa, "n", &m, &n, &k,
           &alpha,
           (void*)COMPLEXP(a), &lda,
           (void*)COMPLEXP(b), &ldb,
           &beta,
           (void*)COMPLEXP(c), &ldc);
  Py_RETURN_NONE;
}

/*NUMPY_API
 * Min
 */
NPY_NO_EXPORT PyObject *
PyArray_Min(PyArrayObject *ap, int axis, PyArrayObject *out)
{
    PyArrayObject *arr;
    PyObject *ret;

    arr=(PyArrayObject *)PyArray_CheckAxis(ap, &axis, 0);
    if (arr == NULL) {
        return NULL;
    }
    ret = PyArray_GenericReduceFunction(arr, n_ops.minimum, axis,
                                        PyArray_DESCR(arr)->type_num, out);
    Py_DECREF(arr);
    return ret;
}

int checkArray(PyObject *check, const char* name) {

  if(check==NULL) {
    PyErr_SetString(PyExc_ValueError, "Unspecified array in kdtree");
    return 1;
  }

  PyArray_Descr *descr = PyArray_DESCR(check);
  if(descr==NULL || descr->kind!='f' || descr->elsize!=sizeof(T)) {
    PyErr_Format(PyExc_TypeError, "Incorrect numpy data type for %s passed to kdtree - must match C %s",name,c_name<T>());
    return 1;
  }
  return 0;

}

int getBitDepth(PyObject *check) {

  if(check==NULL) {
    return 0;
  }

  PyArray_Descr *descr = PyArray_DESCR(check);
  if(descr!=NULL && descr->kind=='f' && descr->elsize==sizeof(float))
      return 32;
  else if(descr!=NULL && descr->kind=='f' && descr->elsize==sizeof(double))
      return 64;
  else
      return 0;

}

static PyObject *
array_index(PyArrayObject *v)
{
    if (!PyArray_ISINTEGER(v) || PyArray_SIZE(v) != 1) {
        PyErr_SetString(PyExc_TypeError, "only integer arrays with "     \
                        "one element can be converted to an index");
        return NULL;
    }
    if (PyArray_NDIM(v) != 0) {
        if (DEPRECATE("converting an array with ndim > 0 to an index"
                      " will result in an error in the future") < 0) {
            return NULL;
        }
    }
    return PyArray_DESCR(v)->f->getitem(PyArray_DATA(v), v);
}

static PyObject *sps_putdatarow(PyObject *self, PyObject *args)
{
  char *spec_version, *array_name;
  int ptype, stype;
  PyObject *in_src;
  PyArrayObject *src;
  int no_items;
  int in_row;


  if (!PyArg_ParseTuple(args, "ssiO", &spec_version, &array_name, &in_row,
            &in_src)) {
    return NULL;
  }

  if (!(src = (PyArrayObject*) PyArray_ContiguousFromObject(in_src,
                NPY_NOTYPE, 1, 1))) {
    struct module_state *st = GETSTATE(self);
    PyErr_SetString(st->SPSError, "Input Array is not a 1 dim array");
    return NULL;
  }

  ptype = PyArray_DESCR(src)->type_num;
  stype = sps_py2type(ptype);

  if (ptype == -1) {
    struct module_state *st = GETSTATE(self);
    PyErr_SetString(st->SPSError, "Type of data in shared memory not supported");
    Py_DECREF(src);
    return NULL;
  }

  no_items = (int) (PyArray_DIMS(src)[0]);

  if (SPS_CopyRowToShared(spec_version, array_name, PyArray_DATA(src), stype,
              in_row, no_items, NULL)
      == -1) {
    struct module_state *st = GETSTATE(self);
    PyErr_SetString(st->SPSError, "Error copying data to shared memory");
    Py_DECREF(src);
    return NULL;
  }else
   Py_DECREF(src);

  Py_INCREF(Py_None);
  return Py_None;
}

static void
_default_copyswapn(void *dst, npy_intp dstride, void *src,
                   npy_intp sstride, npy_intp n, int swap, void *arr)
{
    npy_intp i;
    PyArray_CopySwapFunc *copyswap;
    char *dstptr = dst;
    char *srcptr = src;

    copyswap = PyArray_DESCR(arr)->f->copyswap;

    for (i = 0; i < n; i++) {
        copyswap(dstptr, srcptr, swap, arr);
        dstptr += dstride;
        srcptr += sstride;
    }
}

void _copyTransform(const PyArrayObject *transMat, RDGeom::Transform3D &trans) {
  unsigned int nrows = PyArray_DIM(transMat, 0);
  unsigned int ncols = PyArray_DIM(transMat, 1);
  if ((nrows != 4) || (ncols != 4)) {
    throw_value_error("The transform has to be square matrix, of size 4x4");
  }
  if (PyArray_DESCR(const_cast<PyArrayObject *>(transMat))->type_num !=
      NPY_DOUBLE)
    throw_value_error("Only double arrays allowed for transform object ");

  unsigned int dSize = nrows * nrows;
  const double *inData = reinterpret_cast<const double *>(
      PyArray_DATA(const_cast<PyArrayObject *>(transMat)));
  double *tData = trans.getData();
  memcpy(static_cast<void *>(tData), static_cast<const void *>(inData),
         dSize * sizeof(double));
}

/* Fill in the info structure */
static _buffer_info_t*
_buffer_info_new(PyArrayObject *arr)
{
    _buffer_info_t *info;
    _tmp_string_t fmt = {NULL, 0, 0};
    int k;

    info = malloc(sizeof(_buffer_info_t));
    if (info == NULL) {
        goto fail;
    }

    /* Fill in format */
    if (_buffer_format_string(PyArray_DESCR(arr), &fmt, arr, NULL, NULL) != 0) {
        free(fmt.s);
        goto fail;
    }
    _append_char(&fmt, '\0');
    info->format = fmt.s;

    /* Fill in shape and strides */
    info->ndim = PyArray_NDIM(arr);

    if (info->ndim == 0) {
        info->shape = NULL;
        info->strides = NULL;
    }
    else {
        info->shape = malloc(sizeof(Py_ssize_t) * PyArray_NDIM(arr) * 2 + 1);
        if (info->shape == NULL) {
            goto fail;
        }
        info->strides = info->shape + PyArray_NDIM(arr);
        for (k = 0; k < PyArray_NDIM(arr); ++k) {
            info->shape[k] = PyArray_DIMS(arr)[k];
            info->strides[k] = PyArray_STRIDES(arr)[k];
        }
    }

    return info;

fail:
    free(info);
    return NULL;
}

PyObject* scal(PyObject *self, PyObject *args)
{
  Py_complex alpha;
  PyArrayObject* x;
  if (!PyArg_ParseTuple(args, "DO", &alpha, &x))
    return NULL;
  int n = PyArray_DIMS(x)[0];
  for (int d = 1; d < PyArray_NDIM(x); d++)
    n *= PyArray_DIMS(x)[d];
  int incx = 1;

  if (PyArray_DESCR(x)->type_num == NPY_DOUBLE)
    dscal_(&n, &(alpha.real), DOUBLEP(x), &incx);
  else
    zscal_(&n, &alpha, (void*)COMPLEXP(x), &incx);

  Py_RETURN_NONE;
}

NPY_NO_EXPORT int
_zerofill(PyArrayObject *ret)
{
    if (PyDataType_REFCHK(PyArray_DESCR(ret))) {
        PyObject *zero = PyInt_FromLong(0);
        PyArray_FillObjectArray(ret, zero);
        Py_DECREF(zero);
        if (PyErr_Occurred()) {
            Py_DECREF(ret);
            return -1;
        }
    }
    else {
        npy_intp n = PyArray_NBYTES(ret);
        memset(PyArray_DATA(ret), 0, n);
    }
    return 0;
}

PyObject* multi_dotu(PyObject *self, PyObject *args)
{
  PyArrayObject* a;
  PyArrayObject* b;
  PyArrayObject* c;
  if (!PyArg_ParseTuple(args, "OOO", &a, &b, &c)) 
    return NULL;
  int n0 = PyArray_DIMS(a)[0];
  int n = PyArray_DIMS(a)[1];
  for (int i = 2; i < PyArray_NDIM(a); i++)
    n *= PyArray_DIMS(a)[i];
  int incx = 1;
  int incy = 1;
  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    {
      double *ap = DOUBLEP(a);
      double *bp = DOUBLEP(b);
      double *cp = DOUBLEP(c);

      for (int i = 0; i < n0; i++)
        {
          cp[i] = ddot_(&n, (void*)ap, 
             &incx, (void*)bp, &incy);
          ap += n;
          bp += n;
        }
    }
  else
    {
      double_complex* ap = COMPLEXP(a);
      double_complex* bp = COMPLEXP(b);
      double_complex* cp = COMPLEXP(c);
      for (int i = 0; i < n0; i++)
        {
          cp[i] = 0.0;
          for (int j = 0; j < n; j++)
              cp[i] += ap[j] * bp[j];
          ap += n;
          bp += n;
        }
    }
  Py_RETURN_NONE;
}