Index: numpy/core/src/multiarray/methods.c
===================================================================
--- numpy/core/src/multiarray/methods.c	(revision 7128)
+++ numpy/core/src/multiarray/methods.c	(working copy)
@@ -756,6 +756,40 @@
 
 
 static PyObject *
+array_preparearray(PyArrayObject *self, PyObject *args)
+{
+    PyObject *arr;
+    PyObject *ret;
+
+    if (PyTuple_Size(args) < 1) {
+        PyErr_SetString(PyExc_TypeError,
+                        "only accepts 1 argument");
+        return NULL;
+    }
+    arr = PyTuple_GET_ITEM(args, 0);
+    if (!PyArray_Check(arr)) {
+        PyErr_SetString(PyExc_TypeError,
+                        "can only be called with ndarray object");
+        return NULL;
+    }
+
+    Py_INCREF(PyArray_DESCR(arr));
+    ret = PyArray_NewFromDescr(self->ob_type,
+                               PyArray_DESCR(arr),
+                               PyArray_NDIM(arr),
+                               PyArray_DIMS(arr),
+                               PyArray_STRIDES(arr), PyArray_DATA(arr),
+                               PyArray_FLAGS(arr), (PyObject *)self);
+    if (ret == NULL) {
+        return NULL;
+    }
+    Py_INCREF(arr);
+    PyArray_BASE(ret) = arr;
+    return ret;
+}
+
+
+static PyObject *
 array_wraparray(PyArrayObject *self, PyObject *args)
 {
     PyObject *arr;
@@ -2031,6 +2065,8 @@
     /* for subtypes */
     {"__array__", (PyCFunction)array_getarray,
          METH_VARARGS, NULL},
+    {"__array_prepare__", (PyCFunction)array_preparearray,
+         METH_VARARGS, NULL},
     {"__array_wrap__", (PyCFunction)array_wraparray,
          METH_VARARGS, NULL},
 
Index: numpy/core/src/umath/ufunc_object.c
===================================================================
--- numpy/core/src/umath/ufunc_object.c	(revision 7128)
+++ numpy/core/src/umath/ufunc_object.c	(working copy)
@@ -234,6 +234,122 @@
     "and can't coerce safely to supported types";
 
 /*
+ * This function analyzes the input arguments
+ * and determines an appropriate __array_prepare__ function to call
+ * for the outputs.
+ *
+ * If an output argument is provided, then it is wrapped
+ * with its own __array_prepare__ not with the one determined by
+ * the input arguments.
+ *
+ * if the provided output argument is already an ndarray,
+ * the wrapping function is None (which means no wrapping will
+ * be done --- not even PyArray_Return).
+ *
+ * A NULL is placed in output_wrap for outputs that
+ * should just have PyArray_Return called.
+ */
+static void
+_find_array_prepare(PyObject *args, PyObject **output_wrap, int nin, int nout)
+{
+    Py_ssize_t nargs;
+    int i;
+    int np = 0;
+    PyObject *with_wrap[NPY_MAXARGS], *wraps[NPY_MAXARGS];
+    PyObject *obj, *wrap = NULL;
+
+    nargs = PyTuple_GET_SIZE(args);
+    for (i = 0; i < nin; i++) {
+        obj = PyTuple_GET_ITEM(args, i);
+        if (PyArray_CheckExact(obj) || PyArray_IsAnyScalar(obj)) {
+            continue;
+        }
+        wrap = PyObject_GetAttrString(obj, "__array_prepare__");
+        if (wrap) {
+            if (PyCallable_Check(wrap)) {
+                with_wrap[np] = obj;
+                wraps[np] = wrap;
+                ++np;
+            }
+            else {
+                Py_DECREF(wrap);
+                wrap = NULL;
+            }
+        }
+        else {
+            PyErr_Clear();
+        }
+    }
+    if (np > 0) {
+        /* If we have some wraps defined, find the one of highest priority */
+        wrap = wraps[0];
+        if (np > 1) {
+            double maxpriority = PyArray_GetPriority(with_wrap[0],
+                        PyArray_SUBTYPE_PRIORITY);
+            for (i = 1; i < np; ++i) {
+                double priority = PyArray_GetPriority(with_wrap[i],
+                            PyArray_SUBTYPE_PRIORITY);
+                if (priority > maxpriority) {
+                    maxpriority = priority;
+                    Py_DECREF(wrap);
+                    wrap = wraps[i];
+                }
+                else {
+                    Py_DECREF(wraps[i]);
+                }
+            }
+        }
+    }
+
+    /*
+     * Here wrap is the wrapping function determined from the
+     * input arrays (could be NULL).
+     *
+     * For all the output arrays decide what to do.
+     *
+     * 1) Use the wrap function determined from the input arrays
+     * This is the default if the output array is not
+     * passed in.
+     *
+     * 2) Use the __array_prepare__ method of the output object.
+     * This is special cased for
+     * exact ndarray so that no PyArray_Return is
+     * done in that case.
+     */
+    for (i = 0; i < nout; i++) {
+        int j = nin + i;
+        int incref = 1;
+        output_wrap[i] = wrap;
+        if (j < nargs) {
+            obj = PyTuple_GET_ITEM(args, j);
+            if (obj == Py_None) {
+                continue;
+            }
+            if (PyArray_CheckExact(obj)) {
+                output_wrap[i] = Py_None;
+            }
+            else {
+                PyObject *owrap = PyObject_GetAttrString(obj,
+                            "__array_prepare__");
+                incref = 0;
+                if (!(owrap) || !(PyCallable_Check(owrap))) {
+                    Py_XDECREF(owrap);
+                    owrap = wrap;
+                    incref = 1;
+                    PyErr_Clear();
+                }
+                output_wrap[i] = owrap;
+            }
+        }
+        if (incref) {
+            Py_XINCREF(output_wrap[i]);
+        }
+    }
+    Py_XDECREF(wrap);
+    return;
+}
+
+/*
  * Called for non-NULL user-defined functions.
  * The object should be a CObject pointing to a linked-list of functions
  * storing the function, data, and signature of all user-defined functions.
@@ -1054,6 +1169,7 @@
     npy_intp temp_dims[NPY_MAXDIMS];
     npy_intp *out_dims;
     int out_nd;
+    PyObject *wraparr[NPY_MAXARGS];
 
     /* Check number of arguments */
     nargs = PyTuple_Size(args);
@@ -1332,16 +1448,59 @@
             return -1;
         }
 
-	/* Recover mps[i]. */
-	if (self->core_enabled) {
-	    PyArrayObject *ao = mps[i];
-	    mps[i] = (PyArrayObject *)mps[i]->base;
-	    Py_DECREF(ao);
-	}
+        /* Recover mps[i]. */
+        if (self->core_enabled) {
+            PyArrayObject *ao = mps[i];
+            mps[i] = (PyArrayObject *)mps[i]->base;
+            Py_DECREF(ao);
+        }
 
     }
 
     /*
+     * Use __array_prepare__ on all outputs
+     * if present on one of the input arguments.
+     * If present for multiple inputs:
+     * use __array_prepare__ of input object with largest
+     * __array_priority__ (default = 0.0)
+     *
+     * Exception:  we should not wrap outputs for items already
+     * passed in as output-arguments.  These items should either
+     * be left unwrapped or wrapped by calling their own __array_prepare__
+     * routine.
+     *
+     * For each output argument, wrap will be either
+     * NULL --- call PyArray_Return() -- default if no output arguments given
+     * None --- array-object passed in don't call PyArray_Return
+     * method --- the __array_prepare__ method to call.
+     */
+    _find_array_prepare(args, wraparr, loop->ufunc->nin, loop->ufunc->nout);
+
+    /* wrap outputs */
+    for (i = 0; i < loop->ufunc->nout; i++) {
+        int j = loop->ufunc->nin+i;
+        PyObject *wrap;
+        wrap = wraparr[i];
+        if (wrap != NULL) {
+            if (wrap == Py_None) {
+                Py_DECREF(wrap);
+                continue;
+            }
+            PyObject *res = PyObject_CallFunction(wrap, "O(OOi)",
+                        mps[j], loop->ufunc, args, i);
+            Py_DECREF(wrap);
+            if ((res == NULL) || (res == Py_None)) {
+                if (!PyErr_Occurred()){
+                    PyErr_SetString(PyExc_TypeError,
+                            "__array_prepare__ must return an ndarray or subclass thereof");
+                }
+                return -1;
+            }
+            mps[j] = (PyArrayObject *)res;
+        }
+    }
+
+    /*
      * If any of different type, or misaligned or swapped
      * then must use buffers
      */
