Merge trunk (r2124:2142) to power optimization branch

23 files changed, 566 insertions, 101 deletions

diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py
index f3e25d1d9..ac26cc6cc 100644
--- a/numpy/core/code_generators/generate_umath.py
+++ b/numpy/core/code_generators/generate_umath.py
@@ -9,9 +9,21 @@ class TypeDescription(object):
     def __init__(self, type, f=None, in_=None, out=None):
         self.type = type
         self.func_data = f
+        if in_ is not None:
+            in_ = in_.replace('.', type)
         self.in_ = in_
+        if out is not None:
+            out = out.replace('.', type)
         self.out = out
 
+    def finish_signature(self, nin, nout):
+        if self.in_ is None:
+            self.in_ = self.type * nin
+        assert len(self.in_) == nin
+        if self.out is None:
+            self.out = self.type * nout
+        assert len(self.out) == nout
+
 _fdata_map = dict(f='%sf', d='%s', g='%sl',
                   F='nc_%sf', D='nc_%s', G='nc_%sl')
 def build_func_data(types, f):
@@ -55,6 +67,8 @@ class Ufunc(object):
         self.type_descriptions = []
         for td in type_descriptions:
             self.type_descriptions.extend(td)
+        for td in self.type_descriptions:
+            td.finish_signature(self.nin, self.nout)
 
 #each entry in defdict is
 
@@ -153,26 +167,19 @@ defdict = {
           'compute x**2.',
           TD(nobool_or_obj),
           TD(O, f='Py_square'),
-          ),

-'reciprocal' :     

+          ),
+'reciprocal' :
     Ufunc(1, 1, None,
           'compute 1/x',
           TD(nobool_or_obj),
           TD(O, f='Py_reciprocal'),
-          ),

-'ones_like' :     

+          ),
+'ones_like' :
     Ufunc(1, 1, None,
           'return 1',
           TD(nobool_or_obj),
           TD(O, f='Py_get_one'),
           ),
-#'ipower' :
-#    Ufunc(2, 1, One,
-#          'computes x1**n elementwise, where n is an integer or integer array.',
-#          TD(intflt),
-#          TD(cmplx, f='ipow'),
-#          TD(O, f='PyNumber_Power'),
-#          ),
 'power' :
     Ufunc(2, 1, One,
           'computes x1**x2 elementwise.',
@@ -526,7 +533,10 @@ def make_arrays(funcdict):
     #
     code1list = []
     code2list = []
-    for name, uf in funcdict.iteritems():
+    names = funcdict.keys()
+    names.sort()
+    for name in names:
+        uf = funcdict[name]
         funclist = []
         datalist = []
         siglist = []
@@ -560,17 +570,7 @@ def make_arrays(funcdict):
                 tname = chartoname[t.type].upper()
                 funclist.append('%s_%s' % (tname, name))
 
-            if t.in_ is None:
-                xlist = t.type * uf.nin
-            else:
-                xlist = t.in_
-            for x in xlist:
-                siglist.append('PyArray_%s' % (chartoname[x].upper(),))
-            if t.out is None:
-                xlist = t.type * uf.nout
-            else:
-                xlist = t.out
-            for x in xlist:
+            for x in t.in_ + t.out:
                 siglist.append('PyArray_%s' % (chartoname[x].upper(),))
 
             k += 1
@@ -588,7 +588,10 @@ def make_arrays(funcdict):
 
 def make_ufuncs(funcdict):
     code3list = []
-    for name, uf in funcdict.items():
+    names = funcdict.keys()
+    names.sort()
+    for name in names:
+        uf = funcdict[name]
         mlist = []
         mlist.append(\
 r"""f = PyUFunc_FromFuncAndData(%s_functions, %s_data, %s_signatures, %d,
diff --git a/numpy/core/defmatrix.py b/numpy/core/defmatrix.py
index d506220a2..9234d2f53 100644
--- a/numpy/core/defmatrix.py
+++ b/numpy/core/defmatrix.py
@@ -188,6 +188,17 @@ class matrix(N.ndarray):
     def __str__(self):
         return str(self.__array__())
 
+    def sum(self, axis=None, dtype=None):
+        """Sum the matrix over the given axis.  If the axis is None, sum
+        over all dimensions.  This preserves the orientation of the
+        result as a row or column.
+        """
+        s = N.ndarray.sum(self, axis, dtype)
+        if axis==1:
+            return s.transpose()
+        else:
+            return s
+    
     # Needed becase tolist method expects a[i] 
     #  to have dimension a.ndim-1
     def tolist(self):
diff --git a/numpy/core/include/numpy/fenv/fenv.c b/numpy/core/include/numpy/fenv/fenv.c
new file mode 100644
index 000000000..b15ecdc2f
--- /dev/null
+++ b/numpy/core/include/numpy/fenv/fenv.c
@@ -0,0 +1,38 @@
+/*-
+ * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD$
+ */
+
+#include <sys/types.h>
+#include "fenv.h"
+
+const fenv_t __fe_dfl_env = {
+	0xffff0000,
+	0xffff0000,
+	0xffffffff,
+	{ 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff, 0xff }
+};
diff --git a/numpy/core/include/numpy/fenv/fenv.h b/numpy/core/include/numpy/fenv/fenv.h
new file mode 100644
index 000000000..792995b78
--- /dev/null
+++ b/numpy/core/include/numpy/fenv/fenv.h
@@ -0,0 +1,224 @@
+/*-
+ * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD$
+ */
+
+#ifndef	_FENV_H_
+#define	_FENV_H_
+
+#include <sys/cdefs.h>
+#include <sys/types.h>
+
+typedef struct {
+	__uint32_t	__control;
+	__uint32_t	__status;
+	__uint32_t	__tag;
+	char		__other[16];
+} fenv_t;
+
+typedef	__uint16_t	fexcept_t;
+
+/* Exception flags */
+#define	FE_INVALID	0x01
+#define	FE_DENORMAL	0x02
+#define	FE_DIVBYZERO	0x04
+#define	FE_OVERFLOW	0x08
+#define	FE_UNDERFLOW	0x10
+#define	FE_INEXACT	0x20
+#define	FE_ALL_EXCEPT	(FE_DIVBYZERO | FE_DENORMAL | FE_INEXACT | \
+			 FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW)
+
+/* Rounding modes */
+#define	FE_TONEAREST	0x0000
+#define	FE_DOWNWARD	0x0400
+#define	FE_UPWARD	0x0800
+#define	FE_TOWARDZERO	0x0c00
+#define	_ROUND_MASK	(FE_TONEAREST | FE_DOWNWARD | \
+			 FE_UPWARD | FE_TOWARDZERO)
+
+__BEGIN_DECLS
+
+/* Default floating-point environment */
+extern const fenv_t	__fe_dfl_env;
+#define	FE_DFL_ENV	(&__fe_dfl_env)
+
+#define	__fldcw(__cw)		__asm __volatile("fldcw %0" : : "m" (__cw))
+#define	__fldenv(__env)		__asm __volatile("fldenv %0" : : "m" (__env))
+#define	__fnclex()		__asm __volatile("fnclex")
+#define	__fnstenv(__env)	__asm("fnstenv %0" : "=m" (*(__env)))
+#define	__fnstcw(__cw)		__asm("fnstcw %0" : "=m" (*(__cw)))
+#define	__fnstsw(__sw)		__asm("fnstsw %0" : "=am" (*(__sw)))
+#define	__fwait()		__asm __volatile("fwait")
+
+static __inline int
+feclearexcept(int __excepts)
+{
+	fenv_t __env;
+
+	if (__excepts == FE_ALL_EXCEPT) {
+		__fnclex();
+	} else {
+		__fnstenv(&__env);
+		__env.__status &= ~__excepts;
+		__fldenv(__env);
+	}
+	return (0);
+}
+
+static __inline int
+fegetexceptflag(fexcept_t *__flagp, int __excepts)
+{
+	int __status;
+
+	__fnstsw(&__status);
+	*__flagp = __status & __excepts;
+	return (0);
+}
+
+static __inline int
+fesetexceptflag(const fexcept_t *__flagp, int __excepts)
+{
+	fenv_t __env;
+
+	__fnstenv(&__env);
+	__env.__status &= ~__excepts;
+	__env.__status |= *__flagp & __excepts;
+	__fldenv(__env);
+	return (0);
+}
+
+static __inline int
+feraiseexcept(int __excepts)
+{
+	fexcept_t __ex = __excepts;
+
+	fesetexceptflag(&__ex, __excepts);
+	__fwait();
+	return (0);
+}
+
+static __inline int
+fetestexcept(int __excepts)
+{
+	int __status;
+
+	__fnstsw(&__status);
+	return (__status & __excepts);
+}
+
+static __inline int
+fegetround(void)
+{
+	int __control;
+
+	__fnstcw(&__control);
+	return (__control & _ROUND_MASK);
+}
+
+static __inline int
+fesetround(int __round)
+{
+	int __control;
+
+	if (__round & ~_ROUND_MASK)
+		return (-1);
+	__fnstcw(&__control);
+	__control &= ~_ROUND_MASK;
+	__control |= __round;
+	__fldcw(__control);
+	return (0);
+}
+
+static __inline int
+fegetenv(fenv_t *__envp)
+{
+	int __control;
+
+	/*
+	 * fnstenv masks all exceptions, so we need to save and
+	 * restore the control word to avoid this side effect.
+	 */
+	__fnstcw(&__control);
+	__fnstenv(__envp);
+	__fldcw(__control);
+	return (0);
+}
+
+static __inline int
+feholdexcept(fenv_t *__envp)
+{
+
+	__fnstenv(__envp);
+	__fnclex();
+	return (0);
+}
+
+static __inline int
+fesetenv(const fenv_t *__envp)
+{
+
+	__fldenv(*__envp);
+	return (0);
+}
+
+static __inline int
+feupdateenv(const fenv_t *__envp)
+{
+	int __status;
+
+	__fnstsw(&__status);
+	__fldenv(*__envp);
+	feraiseexcept(__status & FE_ALL_EXCEPT);
+	return (0);
+}
+
+#if __BSD_VISIBLE
+
+static __inline int
+fesetmask(int __mask)
+{
+	int __control;
+
+	__fnstcw(&__control);
+	__mask = (__control | FE_ALL_EXCEPT) & ~__mask;
+	__fldcw(__mask);
+	return (~__control & FE_ALL_EXCEPT);
+}
+
+static __inline int
+fegetmask(void)
+{
+	int __control;
+
+	__fnstcw(&__control);
+	return (~__control & FE_ALL_EXCEPT);
+}
+
+#endif /* __BSD_VISIBLE */
+
+__END_DECLS
+
+#endif	/* !_FENV_H_ */
diff --git a/numpy/core/include/numpy/ufuncobject.h b/numpy/core/include/numpy/ufuncobject.h
index 028ec5338..eb195a4ab 100644
--- a/numpy/core/include/numpy/ufuncobject.h
+++ b/numpy/core/include/numpy/ufuncobject.h
@@ -267,7 +267,7 @@ typedef struct {
 #if defined(__GLIBC__) || defined(__APPLE__) || defined(__MINGW32__)
 #include <fenv.h>
 #elif defined(__CYGWIN__)
-#include <mingw/fenv.h>
+#include "fenv/fenv.c"
 #endif
 
 #define UFUNC_CHECK_STATUS(ret) {                                       \
diff --git a/numpy/core/ma.py b/numpy/core/ma.py
index 92ce83ae8..ae10ca681 100644
--- a/numpy/core/ma.py
+++ b/numpy/core/ma.py
@@ -2142,9 +2142,21 @@ array.fill = _m(not_implemented)
 array.flags = property(_m(not_implemented))
 array.flatten = _m(ravel)
 array.getfield = _m(not_implemented)
-array.max = _m(maximum)
+def _max(a, axis=None):
+    if axis is None:
+        return maximum(a)
+    else:
+        return maximum.reduce(a, axis)
+array.max = _m(_max)
+del _max
+def _min(a, axis=None):
+    if axis is None:
+        return minimum(a)
+    else:
+        return minimum.reduce(a, axis)
+array.min = _m(_min)
+del _min
 array.mean = _m(average)
-array.min = _m(minimum)
 array.nbytes = property(_m(not_implemented))
 array.ndim = _m(not_implemented)
 array.newbyteorder = _m(not_implemented)
diff --git a/numpy/core/setup.py b/numpy/core/setup.py
index 27c06ecbf..dce7477d1 100644
--- a/numpy/core/setup.py
+++ b/numpy/core/setup.py
@@ -151,10 +151,16 @@ def configuration(parent_package='',top_path=None):
             join('src','_isnan.c'),
             join('src','ucsnarrow.c'),
             join('include','numpy','*object.h'),
+            'include/numpy/fenv/fenv.c',
+            'include/numpy/fenv/fenv.h',
 	    join(codegen_dir,'genapi.py'),
 	    join(codegen_dir,'*.txt')
             ]
 
+    # Don't install fenv unless we need them.
+    if sys.platform == 'cygwin':
+        config.add_data_dir('include/numpy/fenv')
+        
     config.add_extension('multiarray',
                          sources = [join('src','multiarraymodule.c'),
                                     generate_config_h,
diff --git a/numpy/core/src/arraymethods.c b/numpy/core/src/arraymethods.c
index 4c077bbcd..fe87009e0 100644
--- a/numpy/core/src/arraymethods.c
+++ b/numpy/core/src/arraymethods.c
@@ -1198,7 +1198,7 @@ static char doc_sum[] = "a.sum(axis=None, dtype=None)\n\n"\
   "2.0\n"\
   ">>> array([0.5, 1.5]).sum(dtype=int32)\n"\
   "1\n"\
-  ">>> array([[0, 1], [0, 5]]).sum()\n"\
+  ">>> array([[0, 1], [0, 5]]).sum(axis=0)\n"\
   "array([0, 6])\n"\
   ">>> array([[0, 1], [0, 5]]).sum(axis=1)\n"\
   "array([1, 5])";
diff --git a/numpy/core/src/arrayobject.c b/numpy/core/src/arrayobject.c
index ef9fdd6e9..ba47ffd19 100644
--- a/numpy/core/src/arrayobject.c
+++ b/numpy/core/src/arrayobject.c
@@ -5771,7 +5771,7 @@ PyArray_CastTo(PyArrayObject *out, PyArrayObject *mp)
 	if (simple) {
 		char *inptr;
 		char *optr = out->data;
-		intp obytes = out->descr->elsize * outsize;
+		intp obytes = out->descr->elsize * mpsize;
 		intp ncopies = outsize / mpsize;
 
 		while(ncopies--) {
diff --git a/numpy/core/src/arraytypes.inc.src b/numpy/core/src/arraytypes.inc.src
index 468728154..0fd105ab1 100644
--- a/numpy/core/src/arraytypes.inc.src
+++ b/numpy/core/src/arraytypes.inc.src
@@ -238,7 +238,7 @@ UNICODE_getitem(char *ip, PyArrayObject *ap)
 	}
 	else buffer = ip;
 #ifdef Py_UNICODE_WIDE
-	obj = PyUnicode_FromUnicode((const PyArray_UCS4 *)buffer, mysize);
+	obj = PyUnicode_FromUnicode((const Py_UNICODE *)buffer, mysize);
 #else
 	/* create new empty unicode object of length mysize*2 */
 	obj = MyPyUnicode_New(mysize*2);
diff --git a/numpy/core/src/multiarraymodule.c b/numpy/core/src/multiarraymodule.c
index b5bc6af6a..481e5709a 100644
--- a/numpy/core/src/multiarraymodule.c
+++ b/numpy/core/src/multiarraymodule.c
@@ -3385,16 +3385,16 @@ _use_inherit(PyArray_Descr *type, PyObject *newobj, int *errflag)
 		PyErr_SetString(PyExc_ValueError,
 				"cannot base a new descriptor on an"\
 				" OBJECT descriptor.");
-		return NULL;
+		goto fail;
 	}
 	new = PyArray_DescrNew(type);
-	if (new == NULL) return NULL;
+	if (new == NULL) goto fail;
 
 	if (new->elsize && new->elsize != conv->elsize) {
 		PyErr_SetString(PyExc_ValueError, 
 				"mismatch in size of old"\
 				"and new data-descriptor");
-		return NULL;
+		goto fail;
 	}
 	new->elsize = conv->elsize;
 	if (conv->fields != Py_None) {
@@ -3404,6 +3404,11 @@ _use_inherit(PyArray_Descr *type, PyObject *newobj, int *errflag)
 	Py_DECREF(conv);
 	*errflag = 0;
 	return new;
+
+ fail:
+	Py_DECREF(conv);
+	return NULL;
+
 }
 
 static PyArray_Descr *
@@ -3448,7 +3453,7 @@ _convert_from_tuple(PyObject *obj)
 	else {
 		/* interpret next item as shape (if it's a tuple)
 		   and reset the type to PyArray_VOID with 
-		   anew fields attribute. 
+		   a new fields attribute. 
 	       */
 		PyArray_Dims shape={NULL,-1};
 		PyArray_Descr *newdescr;
@@ -3460,7 +3465,8 @@ _convert_from_tuple(PyObject *obj)
 			goto fail;
 		}
 		/* If (type, 1) was given, it is equivalent to type... */
-		if (shape.len == 1 && shape.ptr[0] == 1 && PyNumber_Check(val)) {
+		if (shape.len == 1 && shape.ptr[0] == 1 && \
+		    PyNumber_Check(val)) {
 			PyDimMem_FREE(shape.ptr);
 			return type;
 		}
@@ -4203,7 +4209,7 @@ _prepend_ones(PyArrayObject *arr, int nd, int ndmin)
 static char doc_fromobject[] = "array(object, dtype=None, copy=1, fortran=0, "\
         "subok=0,ndmin=0)\n"\
         "will return a new array formed from the given object type given.\n"\
-        "Object can anything with an __array__ method, or any object\n"\
+        "Object can be anything with an __array__ method, or any object\n"\
         "exposing the array interface, or any (nested) sequence.\n"\
         "If no type is given, then the type will be determined as the\n"\
         "minimum type required to hold the objects in the sequence.\n"\
diff --git a/numpy/core/src/scalartypes.inc.src b/numpy/core/src/scalartypes.inc.src
index 2f48ff041..9926cd83c 100644
--- a/numpy/core/src/scalartypes.inc.src
+++ b/numpy/core/src/scalartypes.inc.src
@@ -205,8 +205,9 @@ PyArray_FromScalar(PyObject *scalar, PyArray_Descr *outcode)
 	if (PyArray_ISOBJECT(r)) {
 		Py_INCREF(*((PyObject **)memptr));
 	}
-
+#ifndef Py_UNICODE_WIDE
  finish:	
+#endif
 	if (outcode == NULL) return r;
 	
 	if (outcode->type_num == typecode->type_num) {
diff --git a/numpy/core/src/ufuncobject.c b/numpy/core/src/ufuncobject.c
index 4b6c152a5..e09b1ffd5 100644
--- a/numpy/core/src/ufuncobject.c
+++ b/numpy/core/src/ufuncobject.c
@@ -859,7 +859,7 @@ construct_matrices(PyUFuncLoopObject *loop, PyObject *args, PyArrayObject **mps)
 
         /* Check number of arguments */
         nargs = PyTuple_Size(args);
-        if ((nargs != self->nin) && (nargs != self->nargs)) {
+        if ((nargs < self->nin) || (nargs > self->nargs)) {
                 PyErr_SetString(PyExc_ValueError, 
 				"invalid number of arguments");
                 return -1;
@@ -2507,8 +2507,24 @@ PyUFunc_GenericReduction(PyUFuncObject *self, PyObject *args,
 	
 }
 
-static PyObject *
-_find_array_wrap(PyObject *args)
+/* This function analyzes the input arguments
+   and determines an appropriate __array_wrap__ function to call
+   for the outputs.   
+
+   If an output argument is provided, then it is wrapped 
+   with its own __array_wrap__ not with the one determined by
+   the input arguments. 
+
+   if the provided output argument is already an array,
+   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_wrap(PyObject *args, PyObject **output_wrap, int nin, int nout)  
 {
 	int nargs, i;
 	int np = 0;
@@ -2517,7 +2533,7 @@ _find_array_wrap(PyObject *args)
 	PyObject *obj, *wrap = NULL;
 
 	nargs = PyTuple_GET_SIZE(args);
-	for (i=0; i<nargs; i++) {
+	for (i=0; i<nin; i++) {
 		obj = PyTuple_GET_ITEM(args, i);
 		if (PyArray_CheckExact(obj) || PyBigArray_CheckExact(obj) || \
 		    PyArray_IsAnyScalar(obj))
@@ -2537,21 +2553,71 @@ _find_array_wrap(PyObject *args)
 			PyErr_Clear();
 		}
 	}
-	if (np < 2)
-		return wrap;
-	wrap = wraps[0];
-	maxpriority = PyArray_GetPriority(with_wrap[0], PyArray_SUBTYPE_PRIORITY);
-	for (i = 1; i < np; ++i) {
-		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]);
-		}
-	}
-	return wrap;
+	if (np >= 2) {
+                wrap = wraps[0];
+                maxpriority = PyArray_GetPriority(with_wrap[0], 
+                                                  PyArray_SUBTYPE_PRIORITY);
+                for (i = 1; i < np; ++i) {
+                        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_wrap__ method of the output object
+              passed in. -- 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) || 
+                            PyBigArray_CheckExact(obj)) {
+                                output_wrap[i] = Py_None;
+                        }
+                        else {
+                                PyObject *owrap;
+                                owrap = PyObject_GetAttrString  \
+                                        (obj,"__array_wrap__");
+                                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]);
+                }
+        }
+           
+	return;
 }
 
 static PyObject *
@@ -2561,8 +2627,8 @@ ufunc_generic_call(PyUFuncObject *self, PyObject *args)
 	PyTupleObject *ret;
 	PyArrayObject *mps[MAX_ARGS];
 	PyObject *retobj[MAX_ARGS];
+        PyObject *wraparr[MAX_ARGS];
 	PyObject *res;
-	PyObject *wrap;
         int errval;
         
 	/* Initialize all array objects to NULL to make cleanup easier 
@@ -2582,17 +2648,30 @@ ufunc_generic_call(PyUFuncObject *self, PyObject *args)
 	
 	for(i=0; i<self->nin; i++) Py_DECREF(mps[i]);
 
+
 	/*  Use __array_wrap__ on all outputs 
 	        if present on one of the input arguments.
 	    If present for multiple inputs:
 	        use __array_wrap__ of input object with largest 
 		__array_priority__ (default = 0.0)
 	 */
-	wrap = _find_array_wrap(args);
+
+        /* 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_wrap__
+           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_wrap__ method to call. 
+        */
+        _find_array_wrap(args, wraparr, self->nin, self->nout);
 	
 	/* wrap outputs */
 	for (i=0; i<self->nout; i++) {
 		int j=self->nin+i;
+                PyObject *wrap;
 		/* check to see if any UPDATEIFCOPY flags are set 
 		   which meant that a temporary output was generated 
 		*/
@@ -2603,7 +2682,13 @@ ufunc_generic_call(PyUFuncObject *self, PyObject *args)
 					      back into old */
 			mps[j] = (PyArrayObject *)old;
 		}
+                wrap = wraparr[i];
 		if (wrap != NULL) {
+                        if (wrap == Py_None) {
+                                Py_DECREF(wrap);
+                                retobj[i] = (PyObject *)mps[j];
+                                continue;
+                        }
 			res = PyObject_CallFunction(wrap, "O(OOi)",
 						    mps[j], self, args, i);
 			if (res == NULL && PyErr_ExceptionMatches(PyExc_TypeError)) {
@@ -2619,6 +2704,7 @@ ufunc_generic_call(PyUFuncObject *self, PyObject *args)
 				continue;
 			}
 		}
+                /* default behavior */
 		retobj[i] = PyArray_Return(mps[j]);
 	}
 	
diff --git a/numpy/core/src/umathmodule.c.src b/numpy/core/src/umathmodule.c.src
index fd6f0ab44..896371e65 100644
--- a/numpy/core/src/umathmodule.c.src
+++ b/numpy/core/src/umathmodule.c.src
@@ -7,7 +7,6 @@
 #include "abstract.h"
 #include <math.h>
 
-
 /* A whole slew of basic math functions are provided originally
    by Konrad Hinsen. */
 
diff --git a/numpy/core/tests/test_defmatrix.py b/numpy/core/tests/test_defmatrix.py
index 5b94d1e3c..48fd6079e 100644
--- a/numpy/core/tests/test_defmatrix.py
+++ b/numpy/core/tests/test_defmatrix.py
@@ -25,6 +25,21 @@ class test_ctor(ScipyTestCase):
         assert mvec.shape == (1,5)
         
 class test_properties(ScipyTestCase):
+    def check_sum(self):
+        """Test whether matrix.sum(axis=1) preserves orientation.
+        Fails in NumPy <= 0.9.6.2127.
+        """
+        M = matrix([[1,2,0,0],
+                   [3,4,0,0],
+                   [1,2,1,2],
+                   [3,4,3,4]])
+        sum0 = matrix([8,12,4,6])
+        sum1 = matrix([3,7,6,14]).T
+        sumall = 30
+        assert_array_equal(sum0, M.sum(axis=0))
+        assert_array_equal(sum1, M.sum(axis=1))
+        assert sumall == M.sum()
+
     def check_basic(self):
         import numpy.linalg as linalg
         
@@ -77,7 +92,7 @@ class test_properties(ScipyTestCase):
 	assert A.sum() == matrix(2)
 	assert A.mean() == matrix(0.5)
 
-class test_autocasting(ScipyTestCase):
+class test_casting(ScipyTestCase):
     def check_basic(self):
         A = arange(100).reshape(10,10)
         mA = matrix(A)
diff --git a/numpy/core/tests/test_ma.py b/numpy/core/tests/test_ma.py
index 385ff9866..03274237e 100644
--- a/numpy/core/tests/test_ma.py
+++ b/numpy/core/tests/test_ma.py
@@ -676,8 +676,16 @@ class test_ufuncs(ScipyTestCase):
         self.failUnlessEqual(sum(a[:3]), 0)
         self.failUnlessEqual(product(a), 0)
 
+    def test_minmax(self):
+        a = arange(1,13).reshape(3,4)
+        amask = masked_where(a < 5,a)
+        self.failUnlessEqual(amask.max(), a.max())
+        self.failUnlessEqual(amask.min(), 5)
+        self.failUnless((amask.max(0) == a.max(0)).all())
+        self.failUnless((amask.min(0) == [5,6,7,8]).all())
+        self.failUnless(amask.max(1)[0].mask)
+        self.failUnless(amask.min(1)[0].mask)
         
-
 def eqmask(m1, m2):
     if m1 is nomask:
         return m2 is nomask
diff --git a/numpy/distutils/ccompiler.py b/numpy/distutils/ccompiler.py
index f865416b3..efd458897 100644
--- a/numpy/distutils/ccompiler.py
+++ b/numpy/distutils/ccompiler.py
@@ -10,7 +10,7 @@ from distutils.version import LooseVersion
 
 import log
 from exec_command import exec_command
-from misc_util import cyg2win32, is_sequence
+from misc_util import cyg2win32, is_sequence, mingw32
 from distutils.spawn import _nt_quote_args
 
 # Using customized CCompiler.spawn.
@@ -235,8 +235,7 @@ if sys.platform == 'win32':
    compiler_class['mingw32'] = ('mingw32ccompiler', 'Mingw32CCompiler',
                                 "Mingw32 port of GNU C Compiler for Win32"\
                                 "(for MSC built Python)")
-   if os.environ.get('OSTYPE','')=='msys' or \
-          os.environ.get('MSYSTEM','')=='MINGW32':
+   if mingw32():
        # On windows platforms, we want to default to mingw32 (gcc)
        # because msvc can't build blitz stuff.
        log.info('Setting mingw32 as default compiler for nt.')
diff --git a/numpy/distutils/fcompiler/gnu.py b/numpy/distutils/fcompiler/gnu.py
index ef0426b5e..10bafdf8e 100644
--- a/numpy/distutils/fcompiler/gnu.py
+++ b/numpy/distutils/fcompiler/gnu.py
@@ -7,6 +7,7 @@ import warnings
 from numpy.distutils.cpuinfo import cpu
 from numpy.distutils.fcompiler import FCompiler
 from numpy.distutils.exec_command import exec_command, find_executable
+from numpy.distutils.misc_util import mingw32
 
 class GnuFCompiler(FCompiler):
 
@@ -104,7 +105,10 @@ class GnuFCompiler(FCompiler):
             g2c = self.g2c
         
         if sys.platform=='win32':
-            opt.append('gcc')
+            # To avoid undefined reference __EH_FRAME_BEGIN__ linker error,
+            # don't use -lgcc option for mingw32 g77 linker.
+            if not mingw32():
+                opt.append('gcc')
         if g2c is not None:
             opt.append(g2c)
         if sys.platform == 'darwin':
diff --git a/numpy/distutils/misc_util.py b/numpy/distutils/misc_util.py
index 7fe329385..0de24cc94 100644
--- a/numpy/distutils/misc_util.py
+++ b/numpy/distutils/misc_util.py
@@ -9,7 +9,7 @@ __all__ = ['Configuration', 'get_numpy_include_dirs', 'default_config_dict',
            'dict_append', 'appendpath', 'generate_config_py',
            'get_cmd', 'allpath', 'get_mathlibs',
            'terminal_has_colors', 'red_text', 'green_text', 'yellow_text',
-           'blue_text', 'cyan_text', 'cyg2win32', 'all_strings',
+           'blue_text', 'cyan_text', 'cyg2win32','mingw32','all_strings',
            'has_f_sources', 'has_cxx_sources', 'filter_sources',
            'get_dependencies', 'is_local_src_dir', 'get_ext_source_files',
            'get_script_files', 'get_lib_source_files', 'get_data_files',
@@ -158,6 +158,16 @@ def cyg2win32(path):
         path = path[10] + ':' + os.path.normcase(path[11:])
     return path
 
+def mingw32():
+    """ Return true when using mingw32 environment.
+    """
+    if sys.platform=='win32':
+        if os.environ.get('OSTYPE','')=='msys':
+            return True
+        if os.environ.get('MSYSTEM','')=='MINGW32':
+            return True
+    return False
+
 #########################
 
 #XXX need support for .C that is also C++
diff --git a/numpy/lib/polynomial.py b/numpy/lib/polynomial.py
index 0978c6274..5bd393cd5 100644
--- a/numpy/lib/polynomial.py
+++ b/numpy/lib/polynomial.py
@@ -10,9 +10,9 @@ import re
 import numpy.core.numeric as NX
 
 from numpy.core import isscalar
-from twodim_base import diag, vander
-from shape_base import hstack, atleast_1d
-from function_base import trim_zeros, sort_complex
+from numpy.lib.twodim_base import diag, vander
+from numpy.lib.shape_base import hstack, atleast_1d
+from numpy.lib.function_base import trim_zeros, sort_complex
 eigvals = None
 lstsq = None
 
@@ -355,6 +355,9 @@ class poly1d(object):
     p = poly1d([1,2,3], variable='lambda') will use lambda in the
     string representation of p.
     """
+    coeffs = None
+    order = None
+    variable = None
     def __init__(self, c_or_r, r=0, variable=None):
         if isinstance(c_or_r, poly1d):
             for key in c_or_r.__dict__.keys():
@@ -464,7 +467,7 @@ class poly1d(object):
         if not isscalar(val) or int(val) != val or val < 0:
             raise ValueError, "Power to non-negative integers only."
         res = [1]
-        for k in range(val):
+        for _ in range(val):
             res = polymul(self.coeffs, res)
         return poly1d(res)
 
diff --git a/numpy/lib/scimath.py b/numpy/lib/scimath.py
index 1e06946b6..1875647aa 100644
--- a/numpy/lib/scimath.py
+++ b/numpy/lib/scimath.py
@@ -27,7 +27,7 @@ def _fix_real_lt_zero(x):
     x = asarray(x)
     if any(isreal(x) & (x<0)):
         x = _tocomplex(x)
-    return asscalar(x)
+    return x
 
 def _fix_real_abs_gt_1(x):
     x = asarray(x)
diff --git a/numpy/testing/numpytest.py b/numpy/testing/numpytest.py
index 9ddfbdd8c..127a5d3e2 100644
--- a/numpy/testing/numpytest.py
+++ b/numpy/testing/numpytest.py
@@ -1,5 +1,6 @@
 
 import os
+import re
 import sys
 import imp
 import types
@@ -207,18 +208,27 @@ class NumpyTest:
     Package is supposed to contain a directory tests/
     with test_*.py files where * refers to the names of submodules.
     See .rename() method to redefine name mapping between test_*.py files
-    and names of submodules.
+    and names of submodules. Pattern test_*.py can be overwritten by
+    redefining .get_testfile() method.
 
     test_*.py files are supposed to define a classes, derived
     from NumpyTestCase or unittest.TestCase, with methods having
-    names starting with test or bench or check.
+    names starting with test or bench or check. The names of TestCase
+    classes must have a prefix test. This can be overwritten by
+    redefining .check_testcase_name() method.
 
     And that is it! No need to implement test or test_suite functions
     in each .py file.
 
-    Also old styled test_suite(level=1) hooks are supported but
-    soon to be removed.
+    Also old styled test_suite(level=1) hooks are supported.
     """
+    _check_testcase_name = re.compile(r'test.*').match
+    def check_testcase_name(self, name):
+        return self._check_testcase_name(name) is not None
+
+    def get_testfile(self, short_module_name):
+        return 'test_' + short_module_name + '.py'
+
     def __init__(self, package=None):
         if package is None:
             from numpy.distutils.misc_util import get_frame
@@ -277,13 +287,14 @@ class NumpyTest:
         short_module_name = self._rename_map.get(short_module_name,short_module_name)
         if short_module_name is None:
             return []
+        full_module_name = module.__name__+'.'+short_module_name
 
         test_dir = os.path.join(d,'tests')
-        test_file = os.path.join(test_dir,'test_'+short_module_name+'.py')
+        fn = self.get_testfile(short_module_name)
+        test_file = os.path.join(test_dir,fn)
 
         local_test_dir = os.path.join(os.getcwd(),'tests')
-        local_test_file = os.path.join(local_test_dir,
-                                       'test_'+short_module_name+'.py')
+        local_test_file = os.path.join(local_test_dir, fn)
         if os.path.basename(os.path.dirname(local_test_dir)) \
                == os.path.basename(os.path.dirname(test_dir)) \
            and os.path.isfile(local_test_file):
@@ -308,19 +319,10 @@ class NumpyTest:
             return []
     
         try:
-            if sys.version[:3]=='2.1':
-                # Workaround for Python 2.1 .pyc file generator bug
-                import random
-                pref = '-nopyc'+`random.randint(1,100)`
-            else:
-                pref = ''
             f = open(test_file,'r')
-            test_module = imp.load_module(\
-                module.__name__+'.test_'+short_module_name+pref,
-                f, test_file+pref,('.py', 'r', 1))
+            test_module = imp.load_module(full_module_name,
+                                          f, test_file, ('.py', 'r', 1))
             f.close()
-            if sys.version[:3]=='2.1' and os.path.isfile(test_file+pref+'c'):
-                os.remove(test_file+pref+'c')
         except:
             self.warn('   !! FAILURE importing tests for %s' % mstr(module))
             output_exception(sys.stderr)
@@ -332,21 +334,14 @@ class NumpyTest:
 
     def _get_suite_list(self, test_module, level, module_name='__main__'):
         mstr = self._module_str
-        if hasattr(test_module,'test_suite'):
-            # Using old styled test suite
-            try:
-                total_suite = test_module.test_suite(level)
-                return total_suite._tests
-            except:
-                self.warn('   !! FAILURE building tests for %s' % mstr(test_module))
-                output_exception(sys.stderr)
-                return []
         suite_list = []
+        if hasattr(test_module,'test_suite'):
+            suite_list.extend(test_module.test_suite(level)._tests)
         for name in dir(test_module):
             obj = getattr(test_module, name)
             if type(obj) is not type(unittest.TestCase) \
                or not issubclass(obj, unittest.TestCase) \
-               or obj.__name__[:4] != 'test':
+               or not self.check_testcase_name(obj.__name__):
                 continue
             for mthname in self._get_method_names(obj,level):
                 suite = obj(mthname)
@@ -370,6 +365,8 @@ class NumpyTest:
             if package_name != name[:len(package_name)] \
                    or module is None:
                 continue
+            if not hasattr(module,'__file__'):
+                continue
             if os.path.basename(os.path.dirname(module.__file__))=='tests':
                 continue
             modules.append(module)
diff --git a/numpy/testing/utils.py b/numpy/testing/utils.py
index 19acd2654..4ea19c270 100644
--- a/numpy/testing/utils.py
+++ b/numpy/testing/utils.py
@@ -56,15 +56,58 @@ else:
         """ Return number of jiffies (1/100ths of a second) that this
     process has been scheduled in user mode. [Emulation with time.time]. """
         return int(100*(time.time()-_load_time))
-
     def memusage():
         """ Return memory usage of running python. [Not implemented]"""
         return
 
+if os.name=='nt':
+    # Code stolen from enthought/debug/memusage.py
+    import win32pdh    
+    # from win32pdhutil, part of the win32all package
+    def GetPerformanceAttributes(object, counter, instance = None,
+                                 inum=-1, format = win32pdh.PDH_FMT_LONG, machine=None):
+        # NOTE: Many counters require 2 samples to give accurate results,
+        # including "% Processor Time" (as by definition, at any instant, a
+        # thread's CPU usage is either 0 or 100).  To read counters like this,
+        # you should copy this function, but keep the counter open, and call
+        # CollectQueryData() each time you need to know.
+        # See http://msdn.microsoft.com/library/en-us/dnperfmo/html/perfmonpt2.asp
+        # My older explanation for this was that the "AddCounter" process forced
+        # the CPU to 100%, but the above makes more sense :)
+        path = win32pdh.MakeCounterPath( (machine,object,instance, None, inum,counter) )
+        hq = win32pdh.OpenQuery()
+        try:
+            hc = win32pdh.AddCounter(hq, path)
+            try:
+                win32pdh.CollectQueryData(hq)
+                type, val = win32pdh.GetFormattedCounterValue(hc, format)
+                return val
+            finally:
+                win32pdh.RemoveCounter(hc)
+        finally:
+            win32pdh.CloseQuery(hq)
+    
+    def memusage(processName="python", instance=0):
+        return GetPerformanceAttributes("Process", "Virtual Bytes",
+                                        processName, instance,
+                                        win32pdh.PDH_FMT_LONG, None)
+
 def assert_equal(actual,desired,err_msg='',verbose=1):
     """ Raise an assertion if two items are not
         equal.  I think this should be part of unittest.py
     """
+    if isinstance(desired, dict):
+        assert isinstance(actual, dict),`type(actual)`
+        assert_equal(len(actual),len(desired),err_msg,verbose)
+        for k,i in desired.items():
+            assert actual.has_key(k),`k`
+            assert_equal(actual[k], desired[k], 'key=%r\n%s' % (k,err_msg), verbose)
+        return
+    if isinstance(desired, list) and isinstance(actual, list):
+        assert_equal(len(actual),len(desired),err_msg,verbose)
+        for k in range(len(desired)):
+            assert_equal(actual[k], desired[k], 'item=%r\n%s' % (k,err_msg), verbose)
+        return
     from numpy.core import ArrayType
     if isinstance(actual, ArrayType) or isinstance(desired, ArrayType):
         return assert_array_equal(actual, desired, err_msg)
@@ -79,7 +122,7 @@ def assert_equal(actual,desired,err_msg='',verbose=1):
              + 'DESIRED: ' + repr(desired) \
              + '\nACTUAL: ' + repr(actual)
     assert desired == actual, msg
-
+    return
 
 def assert_almost_equal(actual,desired,decimal=7,err_msg='',verbose=1):
     """ Raise an assertion if two items are not