Index: doc/source/reference/arrays.classes.rst
===================================================================
--- doc/source/reference/arrays.classes.rst	(revision 7180)
+++ doc/source/reference/arrays.classes.rst	(working copy)
@@ -47,15 +47,30 @@
    update meta-information from the "parent." Subclasses inherit a
    default implementation of this method that does nothing.
 
-.. function:: __array_wrap__(array)
+.. function:: __array_prepare__(array, context=None)
 
-   This method should return an instance of the subclass from the
-   :class:`ndarray` object passed in. For example, this is called
-   after every :ref:`ufunc <ufuncs.output-type>` for the object with
-   the highest array priority. The ufunc-computed array object is
-   passed in and whatever is returned is passed to the
-   user. Subclasses inherit a default implementation of this method.
+   At the beginning of every :ref:`ufunc <ufuncs.output-type>`, this
+   method is called on the input object with the highest array
+   priority, or the output object if one was specified. The output
+   array is passed in and whatever is returned is passed to the ufunc.
+   Subclasses inherit a default implementation of this method which
+   simply returns the output array unmodified. Subclasses may opt to
+   use this method to transform the output array into an instance of
+   the subclass and update metadata before returning the array to the
+   ufunc for computation.
 
+.. function:: __array_wrap__(array, context=None)
+
+   At the end of every :ref:`ufunc <ufuncs.output-type>`, this method
+   is called on the input object with the highest array priority, or
+   the output object if one was specified. The ufunc-computed array
+   is passed in and whatever is returned is passed to the user.
+   Subclasses inherit a default implementation of this method, which
+   transforms the array into a new instance of the object's class. Subclasses
+   may opt to use this method to transform the output array into an
+   instance of the subclass and update metadata before returning the
+   array to the user.
+
 .. data:: __array_priority__
 
    The value of this attribute is used to determine what type of
Index: doc/source/reference/ufuncs.rst
===================================================================
--- doc/source/reference/ufuncs.rst	(revision 7180)
+++ doc/source/reference/ufuncs.rst	(working copy)
@@ -102,19 +102,24 @@
 The output of the ufunc (and its methods) is not necessarily an
 :class:`ndarray`, if all input arguments are not :class:`ndarrays <ndarray>`.
 
-All output arrays will be passed to the :obj:`__array_wrap__`
-method of the input (besides :class:`ndarrays <ndarray>`, and scalars)
-that defines it **and** has the highest :obj:`__array_priority__` of
-any other input to the universal function. The default
-:obj:`__array_priority__` of the ndarray is 0.0, and the default
-:obj:`__array_priority__` of a subtype is 1.0. Matrices have
-:obj:`__array_priority__` equal to 10.0.
+All output arrays will be passed to the :obj:`__array_prepare__` and
+:obj:`__array_wrap__` methods of the input (besides
+:class:`ndarrays <ndarray>`, and scalars) that defines it **and** has
+the highest :obj:`__array_priority__` of any other input to the
+universal function. The default :obj:`__array_priority__` of the
+ndarray is 0.0, and the default :obj:`__array_priority__` of a subtype
+is 1.0. Matrices have :obj:`__array_priority__` equal to 10.0.
 
 The ufuncs can also all take output arguments. The output will be cast
 if necessary to the provided output array. If a class with an
 :obj:`__array__` method is used for the output, results will be
 written to the object returned by :obj:`__array__`. Then, if the class
-also has an :obj:`__array_wrap__` method, the returned
+also has an :obj:`__array_prepare__` method, it is called so metadata
+may be determined based on the context of the ufunc (the context
+consisting of the ufunc itself, the arguments passed to the ufunc, and
+the ufunc domain.) The array object returned by
+:obj:`__array_prepare__` is passed to the ufunc for computation.
+Finally, if the class also has an :obj:`__array_wrap__` method, the returned
 :class:`ndarray` result will be passed to that method just before
 passing control back to the caller.
 
Index: numpy/core/src/multiarray/methods.c
===================================================================
--- numpy/core/src/multiarray/methods.c	(revision 7180)
+++ numpy/core/src/multiarray/methods.c	(working copy)
@@ -756,6 +756,31 @@
 
 
 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 (arr == NULL) {
+        return NULL;
+    }
+    if (!PyArray_Check(arr)) {
+        PyErr_SetString(PyExc_TypeError,
+                        "can only be called with ndarray object");
+        return NULL;
+    }
+    Py_INCREF(arr);
+    return arr;
+}
+
+
+static PyObject *
 array_wraparray(PyArrayObject *self, PyObject *args)
 {
     PyObject *arr;
@@ -2031,6 +2056,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 7180)
+++ 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 +1170,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 +1449,60 @@
             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;
+            }
+            Py_DECREF(mps[j]);
+            mps[j] = (PyArrayObject *)res;
+        }
+    }
+
+    /*
      * If any of different type, or misaligned or swapped
      * then must use buffers
      */
Index: numpy/core/tests/test_umath.py
===================================================================
--- numpy/core/tests/test_umath.py	(revision 7180)
+++ numpy/core/tests/test_umath.py	(working copy)
@@ -484,6 +484,38 @@
         a = A()
         self.failUnlessRaises(RuntimeError, ncu.maximum, a, a)
 
+    def test_default_prepare(self):
+        class with_wrap(object):
+            __array_priority__ = 10
+            def __array__(self):
+                return np.zeros(1)
+            def __array_wrap__(self, arr, context):
+                return arr
+        a = with_wrap()
+        x = ncu.minimum(a, a)
+        assert_equal(x, np.zeros(1))
+        assert_equal(type(x), np.ndarray)
+
+    def test_prepare(self):
+        class with_prepare(np.ndarray):
+            __array_priority__ = 10
+            def __array_prepare__(self, arr, context):
+                # make sure we can return a new 
+                return np.array(arr).view(type=with_prepare)
+        a = np.array(1).view(type=with_prepare)
+        x = np.add(a, a)
+        assert_equal(x, np.array(2))
+        assert_equal(type(x), with_prepare)
+
+    def test_failing_prepare(self):
+        class A(object):
+            def __array__(self):
+                return np.zeros(1)
+            def __array_prepare__(self, arr, context=None):
+                raise RuntimeError
+        a = A()
+        self.failUnlessRaises(RuntimeError, ncu.maximum, a, a)
+
     def test_array_with_context(self):
         class A(object):
             def __array__(self, dtype=None, context=None):
Index: numpy/add_newdocs.py
===================================================================
--- numpy/add_newdocs.py	(revision 7180)
+++ numpy/add_newdocs.py	(working copy)
@@ -1645,8 +1645,14 @@
     """))
 
 
+add_newdoc('numpy.core.multiarray', 'ndarray', ('__array_prepare__',
+    """a.__array_prepare__(obj) -> Object of same type as ndarray object obj.
+
+    """))
+
+
 add_newdoc('numpy.core.multiarray', 'ndarray', ('__array_wrap__',
-    """a.__array_wrap__(obj) -> Object of same type as a from ndarray obj.
+    """a.__array_wrap__(obj) -> Object of same type as ndarray object a.
 
     """))
 
@@ -3943,7 +3949,7 @@
 
     """)
 
-   
+
 ##############################################################################
 #
 # Documentation for `generic` attributes and methods
@@ -3955,7 +3961,7 @@
     """)
 
 # Attributes
-    
+
 add_newdoc('numpy.core.numerictypes', 'generic', ('T',
     """
     """))
@@ -4246,8 +4252,8 @@
 add_newdoc('numpy.core.numerictypes', 'generic', ('view',
     """
     """))
-   
 
+
 ##############################################################################
 #
 # Documentation for other scalar classes
Index: numpy/doc/subclassing.py
===================================================================
--- numpy/doc/subclassing.py	(revision 7180)
+++ numpy/doc/subclassing.py	(working copy)
@@ -229,8 +229,8 @@
   <class '__main__.RealisticInfoArray'>
   information
 
-``__array_wrap__`` for ufuncs
------------------------------
+``__array_prepare__`` and ``__array_wrap__`` for ufuncs
+-------------------------------------------------------
 
 Let's say you have an instance ``obj`` of your new subclass,
 ``RealisticInfoArray``, and you pass it into a ufunc with another
@@ -239,17 +239,24 @@
   arr = np.arange(5)
   ret = np.multiply.outer(arr, obj)
 
-When a numpy ufunc is called on a subclass of ndarray, the
-__array_wrap__ method is called to transform the result into a new
-instance of the subclass. By default, __array_wrap__ will call
-__array_finalize__, and the attributes will be inherited.
+When a numpy ufunc is called on a subclass of ndarray, the output
+array ``ret`` is created and passed to the __array_prepare__ method,
+in order to give you a chance to transform the array's type and
+determine metadata before the ufunc actually performs its
+calculations. Once the ufunc has finished, it calls the
+__array_wrap__ method, again giving you a chance to update metadata
+and transform the array's type. By default, __array_prepare__ simply
+returns the array unmodified, and __array_wrap__ calls
+__array_finalize__ to transform the array into an instance of your
+subclass and inherit the attributes from ``obj``.
 
-By defining a specific __array_wrap__ method for our subclass, we can
-tweak the output. The __array_wrap__ method requires one argument, the
-object on which the ufunc is applied, and an optional parameter
-*context*. This parameter is returned by some ufuncs as a 3-element
-tuple: (name of the ufunc, argument of the ufunc, domain of the
-ufunc). See the masked array subclass for an implementation.
+By defining specific __array_prepare and/or __array_wrap__ methods for
+our subclass, we can tweak the output. These methods require one
+argument, the return object on which the ufunc is applied, and an
+optional parameter *context*. This parameter is returned by some ufuncs
+as a 3-element tuple: (name of the ufunc, arguments passed to the ufunc,
+and domain of the ufunc). See the masked array subclass for an
+implementation of __array_wrap__.
 
 Extra gotchas - custom __del__ methods and ndarray.base
 -------------------------------------------------------
