Add qr decomposition to linalg

1 files changed, 160 insertions, 41 deletions

diff --git a/numpy/linalg/lapack_litemodule.c b/numpy/linalg/lapack_litemodule.c
index 0f1211b07..12484480e 100644
--- a/numpy/linalg/lapack_litemodule.c
+++ b/numpy/linalg/lapack_litemodule.c
@@ -16,11 +16,11 @@ typedef struct { float r, i; } f2c_complex;
 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 @@ extern void FNAME(zgeev)(char *jobvl, char *jobvr, int *n,
                          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 @@ lapack_lite_dgeev(PyObject *self, PyObject *args)
     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 @@ lapack_lite_dsyevd(PyObject *self, PyObject *args)
     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 @@ lapack_lite_zheevd(PyObject *self, PyObject *args)
     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 @@ lapack_lite_dgelsd(PyObject *self, PyObject *args)
     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 @@ lapack_lite_dgesv(PyObject *self, PyObject *args)
     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 @@ lapack_lite_dgesdd(PyObject *self, PyObject *args)
     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 @@ lapack_lite_dgetrf(PyObject *self, PyObject *args)
     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 @@ lapack_lite_dpotrf(PyObject *self, PyObject *args)
     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 @@ lapack_lite_zgeev(PyObject *self, PyObject *args)
     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 @@ lapack_lite_zgelsd(PyObject *self, PyObject *args)
     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 @@ lapack_lite_zgesv(PyObject *self, PyObject *args)
     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 @@ lapack_lite_zgesdd(PyObject *self, PyObject *args)
     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 @@ lapack_lite_zgetrf(PyObject *self, PyObject *args)
     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 @@ lapack_lite_zpotrf(PyObject *self, PyObject *args)
 
     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 @@ static struct PyMethodDef lapack_lite_module_methods[] = {
     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}
 };