[Numpy-svn] r2991 - in trunk/numpy: linalg linalg/lapack_lite numarray
numpy-svn at scipy.org
numpy-svn at scipy.org
Thu Aug 10 18:45:36 EDT 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