[Numpy-svn] r2991 - in trunk/numpy: linalg linalg/lapack_lite numarray

numpy-svn at scipy.org numpy-svn at scipy.org
Thu Aug 10 17:45:36 CDT 2006


Author: oliphant
Date: 2006-08-10 17:45:29 -0500 (Thu, 10 Aug 2006)
New Revision: 2991

Added:
   trunk/numpy/numarray/fft.py
   trunk/numpy/numarray/linear_algebra.py
   trunk/numpy/numarray/ma.py
   trunk/numpy/numarray/matrix.py
   trunk/numpy/numarray/mlab.py
   trunk/numpy/numarray/random_array.py
Modified:
   trunk/numpy/linalg/lapack_lite/wrapped_routines
   trunk/numpy/linalg/lapack_litemodule.c
   trunk/numpy/linalg/linalg.py
Log:
Add qr decomposition to linalg

Modified: trunk/numpy/linalg/lapack_lite/wrapped_routines
===================================================================
--- trunk/numpy/linalg/lapack_lite/wrapped_routines	2006-08-10 21:16:20 UTC (rev 2990)
+++ trunk/numpy/linalg/lapack_lite/wrapped_routines	2006-08-10 22:45:29 UTC (rev 2991)
@@ -12,6 +12,8 @@
 zpotrf
 dgesdd
 zgesdd
+dgeqrf
+zgeqrf
 # need this b/c it's not properly declared as external in the BLAS source
 dcabs1
 IGNORE: dlamch

Modified: trunk/numpy/linalg/lapack_litemodule.c
===================================================================
--- trunk/numpy/linalg/lapack_litemodule.c	2006-08-10 21:16:20 UTC (rev 2990)
+++ trunk/numpy/linalg/lapack_litemodule.c	2006-08-10 22:45:29 UTC (rev 2991)
@@ -16,11 +16,11 @@
 typedef struct { double r, i; } f2c_doublecomplex;
 /* typedef long int (*L_fp)(); */
 
-extern void FNAME(dgeev)(char *jobvl, char *jobvr, int *n,
+extern int FNAME(dgeev)(char *jobvl, char *jobvr, int *n,
                          double a[], int *lda, double wr[], double wi[],
                          double vl[], int *ldvl, double vr[], int *ldvr,
                          double work[], int lwork[], int *info);
-extern void FNAME(zgeev)(char *jobvl, char *jobvr, int *n,
+extern int FNAME(zgeev)(char *jobvl, char *jobvr, int *n,
                          f2c_doublecomplex a[], int *lda,
                          f2c_doublecomplex w[],
                          f2c_doublecomplex vl[], int *ldvl,
@@ -28,54 +28,70 @@
                          f2c_doublecomplex work[], int *lwork,
                          double rwork[], int *info);
 
-extern void FNAME(dsyevd)(char *jobz, char *uplo, int *n,
+extern int FNAME(dsyevd)(char *jobz, char *uplo, int *n,
                           double a[], int *lda, double w[], double work[],
                           int *lwork, int iwork[], int *liwork, int *info);
-extern void FNAME(zheevd)(char *jobz, char *uplo, int *n,
+extern int FNAME(zheevd)(char *jobz, char *uplo, int *n,
                           f2c_doublecomplex a[], int *lda,
                           double w[], f2c_doublecomplex work[],
                           int *lwork, double rwork[], int *lrwork, int iwork[],
                           int *liwork, int *info);
 
-extern void FNAME(dgelsd)(int *m, int *n, int *nrhs,
+extern int FNAME(dgelsd)(int *m, int *n, int *nrhs,
                           double a[], int *lda, double b[], int *ldb,
                           double s[], double *rcond, int *rank,
                           double work[], int *lwork, int iwork[], int *info);
-extern void FNAME(zgelsd)(int *m, int *n, int *nrhs,
+extern int FNAME(zgelsd)(int *m, int *n, int *nrhs,
                           f2c_doublecomplex a[], int *lda,
                           f2c_doublecomplex b[], int *ldb,
                           double s[], double *rcond, int *rank,
                           f2c_doublecomplex work[], int *lwork,
                           double rwork[], int iwork[], int *info);
 
-extern void FNAME(dgesv)(int *n, int *nrhs,
+extern int FNAME(dgesv)(int *n, int *nrhs,
                          double a[], int *lda, int ipiv[],
                          double b[], int *ldb, int *info);
-extern void FNAME(zgesv)(int *n, int *nrhs,
+extern int FNAME(zgesv)(int *n, int *nrhs,
                          f2c_doublecomplex a[], int *lda, int ipiv[],
                          f2c_doublecomplex b[], int *ldb, int *info);
 
-extern void FNAME(dgetrf)(int *m, int *n,
+extern int FNAME(dgetrf)(int *m, int *n,
                           double a[], int *lda, int ipiv[], int *info);
-extern void FNAME(zgetrf)(int *m, int *n,
+extern int FNAME(zgetrf)(int *m, int *n,
                           f2c_doublecomplex a[], int *lda, int ipiv[],
                           int *info);
 
-extern void FNAME(dpotrf)(char *uplo, int *n, double a[], int *lda, int *info);
-extern void FNAME(zpotrf)(char *uplo, int *n,
+extern int FNAME(dpotrf)(char *uplo, int *n, double a[], int *lda, int *info);
+extern int FNAME(zpotrf)(char *uplo, int *n,
                           f2c_doublecomplex a[], int *lda, int *info);
 
-extern void FNAME(dgesdd)(char *jobz, int *m, int *n,
+extern int FNAME(dgesdd)(char *jobz, int *m, int *n,
                           double a[], int *lda, double s[], double u[],
                           int *ldu, double vt[], int *ldvt, double work[],
                           int *lwork, int iwork[], int *info);
-extern void FNAME(zgesdd)(char *jobz, int *m, int *n,
+extern int FNAME(zgesdd)(char *jobz, int *m, int *n,
                           f2c_doublecomplex a[], int *lda,
                           double s[], f2c_doublecomplex u[], int *ldu,
                           f2c_doublecomplex vt[], int *ldvt,
                           f2c_doublecomplex work[], int *lwork,
                           double rwork[], int iwork[], int *info);
 
+extern int FNAME(dgeqrf)(int *m, int *n, double a[], int *lda,
+                          double tau[], double work[],
+                          int *lwork, int *info);
+
+extern int FNAME(zgeqrf)(int *m, int *n, f2c_doublecomplex a[], int *lda,
+                          f2c_doublecomplex tau[], f2c_doublecomplex work[],
+                          int *lwork, int *info);
+
+extern int FNAME(dorgqr)(int *m, int *n, int *k, double a[], int *lda,
+                          double tau[], double work[], 
+                          int *lwork, int *info);
+
+extern int FNAME(zungqr)(int *m, int *n, int *k, f2c_doublecomplex a[],
+                          int *lda, f2c_doublecomplex tau[], 
+                          f2c_doublecomplex work[], int *lwork, int *info);
+
 static PyObject *LapackError;
 
 #define TRY(E) if (!(E)) return NULL
@@ -141,9 +157,10 @@
     TRY(check_object(vr,PyArray_DOUBLE,"vr","PyArray_DOUBLE","dgeev"));
     TRY(check_object(work,PyArray_DOUBLE,"work","PyArray_DOUBLE","dgeev"));
 
-    lapack_lite_status__ = 0;
-    FNAME(dgeev)(&jobvl,&jobvr,&n,DDATA(a),&lda,DDATA(wr),DDATA(wi),
-                 DDATA(vl),&ldvl,DDATA(vr),&ldvr,DDATA(work),&lwork,&info);
+    lapack_lite_status__ = \
+            FNAME(dgeev)(&jobvl,&jobvr,&n,DDATA(a),&lda,DDATA(wr),DDATA(wi),
+                         DDATA(vl),&ldvl,DDATA(vr),&ldvr,DDATA(work),&lwork,
+                         &info);
 
     return Py_BuildValue("{s:i,s:c,s:c,s:i,s:i,s:i,s:i,s:i,s:i}","dgeev_",
                          lapack_lite_status__,"jobvl",jobvl,"jobvr",jobvr,
@@ -224,9 +241,9 @@
     TRY(check_object(work,PyArray_DOUBLE,"work","PyArray_DOUBLE","dsyevd"));
     TRY(check_object(iwork,PyArray_INT,"iwork","PyArray_INT","dsyevd"));
 
-    lapack_lite_status__ = 0;
-    FNAME(dsyevd)(&jobz,&uplo,&n,DDATA(a),&lda,DDATA(w),DDATA(work),
-                  &lwork,IDATA(iwork),&liwork,&info);
+    lapack_lite_status__ = \
+            FNAME(dsyevd)(&jobz,&uplo,&n,DDATA(a),&lda,DDATA(w),DDATA(work),
+                          &lwork,IDATA(iwork),&liwork,&info);
 
     return Py_BuildValue("{s:i,s:c,s:c,s:i,s:i,s:i,s:i,s:i}","dsyevd_",
                          lapack_lite_status__,"jobz",jobz,"uplo",uplo,
@@ -311,7 +328,7 @@
     TRY(check_object(w,PyArray_DOUBLE,"rwork","PyArray_DOUBLE","zheevd"));
     TRY(check_object(iwork,PyArray_INT,"iwork","PyArray_INT","zheevd"));
 
-    lapack_lite_status__ = 0;
+    lapack_lite_status__ = \
     FNAME(zheevd)(&jobz,&uplo,&n,ZDATA(a),&lda,DDATA(w),ZDATA(work),
                   &lwork,DDATA(rwork),&lrwork,IDATA(iwork),&liwork,&info);
 
@@ -349,10 +366,11 @@
     TRY(check_object(work,PyArray_DOUBLE,"work","PyArray_DOUBLE","dgelsd"));
     TRY(check_object(iwork,PyArray_INT,"iwork","PyArray_INT","dgelsd"));
 
-    lapack_lite_status__ = 0;
-    FNAME(dgelsd)(&m,&n,&nrhs,DDATA(a),&lda,DDATA(b),&ldb,
-                  DDATA(s),&rcond,&rank,DDATA(work),&lwork,IDATA(iwork),&info);
-
+    lapack_lite_status__ = \
+            FNAME(dgelsd)(&m,&n,&nrhs,DDATA(a),&lda,DDATA(b),&ldb,
+                          DDATA(s),&rcond,&rank,DDATA(work),&lwork,
+                          IDATA(iwork),&info);
+    
     return Py_BuildValue("{s:i,s:i,s:i,s:i,s:i,s:i,s:d,s:i,s:i,s:i}","dgelsd_",
                          lapack_lite_status__,"m",m,"n",n,"nrhs",nrhs,
                          "lda",lda,"ldb",ldb,"rcond",rcond,"rank",rank,
@@ -377,7 +395,7 @@
     TRY(check_object(ipiv,PyArray_INT,"ipiv","PyArray_INT","dgesv"));
     TRY(check_object(b,PyArray_DOUBLE,"b","PyArray_DOUBLE","dgesv"));
 
-    lapack_lite_status__ = 0;
+    lapack_lite_status__ = \
     FNAME(dgesv)(&n,&nrhs,DDATA(a),&lda,IDATA(ipiv),DDATA(b),&ldb,&info);
 
     return Py_BuildValue("{s:i,s:i,s:i,s:i,s:i,s:i}","dgesv_",
@@ -414,9 +432,10 @@
     TRY(check_object(work,PyArray_DOUBLE,"work","PyArray_DOUBLE","dgesdd"));
     TRY(check_object(iwork,PyArray_INT,"iwork","PyArray_INT","dgesdd"));
 
-    lapack_lite_status__ = 0;
-    FNAME(dgesdd)(&jobz,&m,&n,DDATA(a),&lda,DDATA(s),DDATA(u),&ldu,
-                  DDATA(vt),&ldvt,DDATA(work),&lwork,IDATA(iwork),&info);
+    lapack_lite_status__ = \
+            FNAME(dgesdd)(&jobz,&m,&n,DDATA(a),&lda,DDATA(s),DDATA(u),&ldu,
+                          DDATA(vt),&ldvt,DDATA(work),&lwork,IDATA(iwork),
+                          &info);
 
     return Py_BuildValue("{s:i,s:c,s:i,s:i,s:i,s:i,s:i,s:i,s:i}","dgesdd_",
                          lapack_lite_status__,"jobz",jobz,"m",m,"n",n,
@@ -439,8 +458,8 @@
     TRY(check_object(a,PyArray_DOUBLE,"a","PyArray_DOUBLE","dgetrf"));
     TRY(check_object(ipiv,PyArray_INT,"ipiv","PyArray_INT","dgetrf"));
 
-    lapack_lite_status__ = 0;
-    FNAME(dgetrf)(&m,&n,DDATA(a),&lda,IDATA(ipiv),&info);
+    lapack_lite_status__ = \
+            FNAME(dgetrf)(&m,&n,DDATA(a),&lda,IDATA(ipiv),&info);
 
     return Py_BuildValue("{s:i,s:i,s:i,s:i,s:i}","dgetrf_",lapack_lite_status__,
                          "m",m,"n",n,"lda",lda,"info",info);
@@ -459,14 +478,61 @@
     TRY(PyArg_ParseTuple(args,"ciOii",&uplo,&n,&a,&lda,&info));
     TRY(check_object(a,PyArray_DOUBLE,"a","PyArray_DOUBLE","dpotrf"));
 
-    lapack_lite_status__ = 0;
-    FNAME(dpotrf)(&uplo,&n,DDATA(a),&lda,&info);
+    lapack_lite_status__ = \
+            FNAME(dpotrf)(&uplo,&n,DDATA(a),&lda,&info);
 
     return Py_BuildValue("{s:i,s:i,s:i,s:i}","dpotrf_",lapack_lite_status__,
                          "n",n,"lda",lda,"info",info);
 }
 
 static PyObject *
+lapack_lite_dgeqrf(PyObject *self, PyObject *args)
+{
+        int  lapack_lite_status__;
+        int m, n, lwork;
+        PyObject *a, *tau, *work;
+        int lda;
+        int info;
+
+        TRY(PyArg_ParseTuple(args,"llOlOOll",&m,&n,&a,&lda,&tau,&work,&lwork,&info));
+
+        /* check objects and convert to right storage order */
+        TRY(check_object(a,PyArray_DOUBLE,"a","PyArray_DOUBLE","dgeqrf"));
+        TRY(check_object(tau,PyArray_DOUBLE,"tau","PyArray_DOUBLE","dgeqrf"));
+        TRY(check_object(work,PyArray_DOUBLE,"work","PyArray_DOUBLE","dgeqrf"));
+
+        lapack_lite_status__ = \
+                FNAME(dgeqrf)(&m, &n, DDATA(a), &lda, DDATA(tau), 
+                              DDATA(work), &lwork, &info);
+
+        return Py_BuildValue("{s:l,s:l,s:l,s:l,s:l,s:l}","dgeqrf_",
+                             lapack_lite_status__,"m",m,"n",n,"lda",lda,
+                             "lwork",lwork,"info",info);
+}
+
+
+static PyObject *
+lapack_lite_dorgqr(PyObject *self, PyObject *args)
+{
+        int  lapack_lite_status__;
+        int m, n, k, lwork;
+        PyObject *a, *tau, *work;
+        int lda;
+        int info;
+        
+        TRY(PyArg_ParseTuple(args,"lllOlOOll",  &m, &n, &k, &a, &lda, &tau, &work, &lwork, &info));
+        TRY(check_object(a,PyArray_DOUBLE,"a","PyArray_DOUBLE","dorgqr"));
+        TRY(check_object(tau,PyArray_DOUBLE,"tau","PyArray_DOUBLE","dorgqr"));
+        TRY(check_object(work,PyArray_DOUBLE,"work","PyArray_DOUBLE","dorgqr"));
+        lapack_lite_status__ = \
+        FNAME(dorgqr)(&m, &n, &k, DDATA(a), &lda, DDATA(tau), DDATA(work), &lwork, &info);
+
+        return Py_BuildValue("{s:l,s:l}","dorgqr_",lapack_lite_status__,
+                             "info",info);
+}
+
+
+static PyObject *
 lapack_lite_zgeev(PyObject *self, PyObject *args)
 {
     int  lapack_lite_status__;
@@ -495,9 +561,10 @@
     TRY(check_object(work,PyArray_CDOUBLE,"work","PyArray_CDOUBLE","zgeev"));
     TRY(check_object(rwork,PyArray_DOUBLE,"rwork","PyArray_DOUBLE","zgeev"));
 
-    lapack_lite_status__ = 0;
-    FNAME(zgeev)(&jobvl,&jobvr,&n,ZDATA(a),&lda,ZDATA(w),ZDATA(vl),&ldvl,
-                 ZDATA(vr),&ldvr,ZDATA(work),&lwork,DDATA(rwork),&info);
+    lapack_lite_status__ = \
+            FNAME(zgeev)(&jobvl,&jobvr,&n,ZDATA(a),&lda,ZDATA(w),ZDATA(vl),
+                         &ldvl,ZDATA(vr),&ldvr,ZDATA(work),&lwork,
+                         DDATA(rwork),&info);
 
     return Py_BuildValue("{s:i,s:c,s:c,s:i,s:i,s:i,s:i,s:i,s:i}","zgeev_",
                          lapack_lite_status__,"jobvl",jobvl,"jobvr",jobvr,
@@ -535,7 +602,7 @@
     TRY(check_object(rwork,PyArray_DOUBLE,"rwork","PyArray_DOUBLE","zgelsd"));
     TRY(check_object(iwork,PyArray_INT,"iwork","PyArray_INT","zgelsd"));
 
-    lapack_lite_status__ = 0;
+    lapack_lite_status__ = \
     FNAME(zgelsd)(&m,&n,&nrhs,ZDATA(a),&lda,ZDATA(b),&ldb,DDATA(s),&rcond,
                   &rank,ZDATA(work),&lwork,DDATA(rwork),IDATA(iwork),&info);
 
@@ -562,7 +629,7 @@
     TRY(check_object(ipiv,PyArray_INT,"ipiv","PyArray_INT","zgesv"));
     TRY(check_object(b,PyArray_CDOUBLE,"b","PyArray_CDOUBLE","zgesv"));
 
-    lapack_lite_status__ = 0;
+    lapack_lite_status__ = \
     FNAME(zgesv)(&n,&nrhs,ZDATA(a),&lda,IDATA(ipiv),ZDATA(b),&ldb,&info);
 
     return Py_BuildValue("{s:i,s:i,s:i,s:i,s:i,s:i}","zgesv_",
@@ -601,7 +668,7 @@
     TRY(check_object(rwork,PyArray_DOUBLE,"rwork","PyArray_DOUBLE","zgesdd"));
     TRY(check_object(iwork,PyArray_INT,"iwork","PyArray_INT","zgesdd"));
 
-    lapack_lite_status__ = 0;
+    lapack_lite_status__ = \
     FNAME(zgesdd)(&jobz,&m,&n,ZDATA(a),&lda,DDATA(s),ZDATA(u),&ldu,
                   ZDATA(vt),&ldvt,ZDATA(work),&lwork,DDATA(rwork),
                   IDATA(iwork),&info);
@@ -627,7 +694,7 @@
     TRY(check_object(a,PyArray_CDOUBLE,"a","PyArray_CDOUBLE","zgetrf"));
     TRY(check_object(ipiv,PyArray_INT,"ipiv","PyArray_INT","zgetrf"));
 
-    lapack_lite_status__ = 0;
+    lapack_lite_status__ = \
     FNAME(zgetrf)(&m,&n,ZDATA(a),&lda,IDATA(ipiv),&info);
 
     return Py_BuildValue("{s:i,s:i,s:i,s:i,s:i}","zgetrf_",
@@ -646,13 +713,61 @@
 
     TRY(PyArg_ParseTuple(args,"ciOii",&uplo,&n,&a,&lda,&info));
     TRY(check_object(a,PyArray_CDOUBLE,"a","PyArray_CDOUBLE","zpotrf"));
-    lapack_lite_status__ = 0;
+    lapack_lite_status__ = \
     FNAME(zpotrf)(&uplo,&n,ZDATA(a),&lda,&info);
 
     return Py_BuildValue("{s:i,s:i,s:i,s:i}","zpotrf_",
                          lapack_lite_status__,"n",n,"lda",lda,"info",info);
 }
 
+static PyObject *
+lapack_lite_zgeqrf(PyObject *self, PyObject *args)
+{
+        int  lapack_lite_status__;
+        int m, n, lwork;
+        PyObject *a, *tau, *work;
+        int lda;
+        int info;
+        
+        TRY(PyArg_ParseTuple(args,"llOlOOll",&m,&n,&a,&lda,&tau,&work,&lwork,&info));
+
+/* check objects and convert to right storage order */
+        TRY(check_object(a,PyArray_CDOUBLE,"a","PyArray_CDOUBLE","zgeqrf"));
+        TRY(check_object(tau,PyArray_CDOUBLE,"tau","PyArray_CDOUBLE","zgeqrf"));
+        TRY(check_object(work,PyArray_CDOUBLE,"work","PyArray_CDOUBLE","zgeqrf"));
+
+        lapack_lite_status__ = \
+        FNAME(zgeqrf)(&m, &n, ZDATA(a), &lda, ZDATA(tau), ZDATA(work), &lwork, &info);
+        
+        return Py_BuildValue("{s:l,s:l,s:l,s:l,s:l,s:l}","zgeqrf_",lapack_lite_status__,"m",m,"n",n,"lda",lda,"lwork",lwork,"info",info);
+}
+
+
+static PyObject *
+lapack_lite_zungqr(PyObject *self, PyObject *args)
+{
+        int  lapack_lite_status__;
+        int m, n, k, lwork;
+        PyObject *a, *tau, *work;
+        int lda;
+        int info;
+
+        TRY(PyArg_ParseTuple(args,"lllOlOOll",  &m, &n, &k, &a, &lda, &tau, &work, &lwork, &info));
+        TRY(check_object(a,PyArray_CDOUBLE,"a","PyArray_CDOUBLE","zungqr"));
+        TRY(check_object(tau,PyArray_CDOUBLE,"tau","PyArray_CDOUBLE","zungqr"));
+        TRY(check_object(work,PyArray_CDOUBLE,"work","PyArray_CDOUBLE","zungqr"));
+
+
+        lapack_lite_status__ = \
+        FNAME(zungqr)(&m, &n, &k, ZDATA(a), &lda, ZDATA(tau), ZDATA(work), 
+                      &lwork, &info);
+        
+        return Py_BuildValue("{s:l,s:l}","zungqr_",lapack_lite_status__,
+                             "info",info);
+}
+
+
+
 #define STR(x) #x
 #define lameth(name) {STR(name), lapack_lite_##name, METH_VARARGS, NULL}
 static struct PyMethodDef lapack_lite_module_methods[] = {
@@ -664,12 +779,16 @@
     lameth(dgesdd),
     lameth(dgetrf),
     lameth(dpotrf),
+    lameth(dgeqrf),
+    lameth(dorgqr),
     lameth(zgeev),
     lameth(zgelsd),
     lameth(zgesv),
     lameth(zgesdd),
     lameth(zgetrf),
     lameth(zpotrf),
+    lameth(zgeqrf),
+    lameth(zungqr),
     { NULL,NULL,0}
 };
 

Modified: trunk/numpy/linalg/linalg.py
===================================================================
--- trunk/numpy/linalg/linalg.py	2006-08-10 21:16:20 UTC (rev 2990)
+++ trunk/numpy/linalg/linalg.py	2006-08-10 22:45:29 UTC (rev 2991)
@@ -12,6 +12,7 @@
            'eigvalsh', 'pinv',
            'det', 'svd',
            'eig', 'eigh','lstsq', 'norm',
+           'qr',
            'LinAlgError'
            ]
 
@@ -177,6 +178,95 @@
         return s.astype(result_t)
     return s
 
+# QR decompostion
+
+def qr(a, mode='full'):
+    """cacluates A=QR, Q orthonormal, R upper triangular
+    
+    mode:  'full' --> (Q,R)
+           'r'    --> R
+           'economic' --> A2 where the diagonal + upper triangle
+                   part of A2 is R. This is faster if you only need R
+    """
+    _assertRank2(a)
+    t, result_t = _commonType(a)
+    a = _fastCopyAndTranspose(t, a)
+    m,n = a.shape
+    mn = min(m,n)
+    tau = zeros((mn,), t)
+    if isComplexType(t):
+        lapack_routine = lapack_lite.zgeqrf
+        routine_name = 'zgeqrf'
+    else:
+        lapack_routine = lapack_lite.dgeqrf
+        routine_name = 'dgeqrf'        
+        
+    # calculate optimal size of work data 'work'
+    lwork = 1
+    work = zeros((lwork,), t)
+    results=lapack_routine(m, n, a, m, tau, work, -1, 0)
+    if results['info'] > 0:
+        raise LinAlgError, '%s returns %d' % (routine_name, results['info'])
+
+    # do qr decomposition
+    lwork = int(abs(work[0]))
+    work = zeros((lwork,),t)
+    results=lapack_routine(m, n, a, m, tau, work, lwork, 0)
+
+    if results['info'] > 0:
+        raise LinAlgError, '%s returns %d' % (routine_name, results['info'])
+
+    #  atemp: convert fortrag storing order to num storing order
+    atemp = a.transpose()
+
+    if atemp.dtype != result_t:
+        atemp = atemp.astype(result_t)
+
+    #  economic mode
+    if mode[0]=='e':
+        return atemp
+
+    #  generate r
+    r = zeros((mn,n), result_t)
+    for i in range(mn):
+            r[i, i:] = atemp[i, i:]
+
+    #  'r'-mode, that is, calculate only r
+    if mode[0]=='r':
+        return r
+
+    #  from here on: build orthonormal matrix q from a
+
+    if isComplexType(t):
+        lapack_routine = lapack_lite.zungqr
+        routine_name = 'zungqr'                
+    else:
+        lapack_routine = lapack_lite.dorgqr
+        routine_name = 'dorgqr'
+
+    # determine optimal lwork
+    lwork = 1
+    work=zeros((lwork,), t)
+    results=lapack_routine(m,mn,mn, a, m, tau, work, -1, 0)
+    if results['info'] > 0:
+        raise LinAlgError, '%s returns %d' % (routine_name, results['info'])
+
+    # compute q
+    lwork = int(abs(work[0]))
+    work=zeros((lwork,), t)
+    results=lapack_routine(m,mn,mn, a, m, tau, work, lwork, 0)
+
+    if results['info'] > 0:
+        raise LinAlgError, '%s returns %d' % (routine_name, results['info'])
+    
+    q = a[:mn,:].transpose()
+
+    if (q.dtype != result_t):
+        q = q.astype(result_t)
+
+    return q,r
+
+
 # Eigenvalues
 def eigvals(a):
     _assertRank2(a)

Added: trunk/numpy/numarray/fft.py
===================================================================

Added: trunk/numpy/numarray/linear_algebra.py
===================================================================
--- trunk/numpy/numarray/linear_algebra.py	2006-08-10 21:16:20 UTC (rev 2990)
+++ trunk/numpy/numarray/linear_algebra.py	2006-08-10 22:45:29 UTC (rev 2991)
@@ -0,0 +1,15 @@
+
+from numpy.oldnumeric.linear_algebra import *
+
+import numpy.oldnumeric.linear_algebra as nol
+
+__all__ = nol.__all__
+__all__ += ['qr_decomposition']
+
+from numpy.linalg import qr as _qr
+
+def qr_decomposition(a, mode='full'):
+    res = _qr(a, mode)
+    if mode == 'full':
+        return res
+    return (None, res)

Added: trunk/numpy/numarray/ma.py
===================================================================

Added: trunk/numpy/numarray/matrix.py
===================================================================

Added: trunk/numpy/numarray/mlab.py
===================================================================

Added: trunk/numpy/numarray/random_array.py
===================================================================



More information about the Numpy-svn mailing list