{"action":"created","comment":{"url":"https://api.github.com/repos/scipy/scipy/pulls/comments/26446344","id":26446344,"diff_hunk":"@@ -0,0 +1,2416 @@\n+\"\"\"\n+Routines for updating QR decompositions\n+\n+.. versionadded: 0.16.0\n+\n+\"\"\"\n+# \n+# Copyright (C) 2014 Eric Moore\n+#\n+# A few references for Updating QR factorizations:\n+# \n+# 1, 2, and 3 cover updates to full decompositons (q is square) and 4 and 5\n+# cover updates to thin (economic) decompositions (r is square). Reference 3\n+# additonally covers updating complete orthogonal factorizations and cholesky\n+# decompositions (i.e. updating R alone).\n+#\n+# 1. Golub, G. H. & Loan, C. F. van V. Matrix Computations, 3rd Ed.\n+#    (Johns Hopkins University Press, 1996).\n+#\n+# 2. Hammarling, S. & Lucas, C. Updating the QR factorization and the least\n+#    squares problem. 1-73 (The University of Manchester, 2008).\n+#    at <http://eprints.ma.man.ac.uk/1192/>\n+#\n+# 3. Gill, P. E., Golub, G. H., Murray, W. & Saunders, M. A. Methods for\n+#    modifying matrix factorizations. Math. Comp. 28, 505-535 (1974).\n+#\n+# 4. Daniel, J. W., Gragg, W. B., Kaufman, L. & Stewart, G. W.\n+#    Reorthogonalization and stable algorithms for updating the Gram-Schmidt QR\n+#    factorization. Math. Comput. 30, 772-795 (1976).\n+#\n+# 5. Reichel, L. & Gragg, W. B. Algorithm 686: FORTRAN Subroutines for\n+#    Updating the QR Decomposition. ACM Trans. Math. Softw. 16, 369–377 (1990).\n+# \n+\n+cimport cython\n+cimport libc.stdlib\n+cimport libc.limits as limits\n+from libc.math cimport sqrt, fabs, hypot\n+from libc.string cimport memset\n+cimport numpy as cnp\n+\n+cdef int NPY_ANYORDER = 0  \n+cdef int MEMORY_ERROR = limits.INT_MAX\n+\n+# These are commented out in the numpy support we cimported above.\n+# Here I have declared them as taking void* instead of PyArrayDescr\n+# and object. In this file, only NULL is passed to these parameters.\n+cdef extern from *:\n+    cnp.ndarray PyArray_CheckFromAny(object, void*, int, int, int, void*)\n+    cnp.ndarray PyArray_FromArray(cnp.ndarray, void*, int)\n+\n+from scipy.linalg cimport blas_pointers\n+from scipy.linalg cimport lapack_pointers\n+\n+import numpy as np\n+\n+#------------------------------------------------------------------------------\n+# These are a set of fused type wrappers around the BLAS and LAPACK calls used. \n+#------------------------------------------------------------------------------\n+\n+ctypedef float complex float_complex\n+ctypedef double complex double_complex\n+ctypedef fused blas_t:\n+    float\n+    double\n+    float_complex\n+    double_complex\n+\n+cdef inline blas_t* index2(blas_t* a, int* as, int i, int j) nogil:\n+    return a + i*as[0] + j*as[1]\n+\n+cdef inline blas_t* index1(blas_t* a, int* as, int i) nogil:\n+    return a + i*as[0]\n+\n+cdef inline blas_t* row(blas_t* a, int* as, int i) nogil:\n+    return a + i*as[0]\n+\n+cdef inline blas_t* col(blas_t* a, int* as, int j) nogil:\n+    return a + j*as[1]\n+    \n+cdef inline void copy(int n, blas_t* x, int incx, blas_t* y, int incy) nogil:\n+    if blas_t is float:\n+        blas_pointers.scopy(&n, x, &incx, y, &incy) \n+    elif blas_t is double:\n+        blas_pointers.dcopy(&n, x, &incx, y, &incy)\n+    elif blas_t is float_complex:\n+        blas_pointers.ccopy(&n, x, &incx, y, &incy)\n+    else:\n+        blas_pointers.zcopy(&n, x, &incx, y, &incy)\n+\n+cdef inline void swap(int n, blas_t* x, int incx, blas_t* y, int incy) nogil:\n+    if blas_t is float:\n+        blas_pointers.sswap(&n, x, &incx, y, &incy) \n+    elif blas_t is double:\n+        blas_pointers.dswap(&n, x, &incx, y, &incy)\n+    elif blas_t is float_complex:\n+        blas_pointers.cswap(&n, x, &incx, y, &incy)\n+    else:\n+        blas_pointers.zswap(&n, x, &incx, y, &incy)\n+\n+cdef inline void scal(int n, blas_t a, blas_t* x, int incx) nogil:\n+    if blas_t is float:\n+        blas_pointers.sscal(&n, &a, x, &incx)\n+    elif blas_t is double:\n+        blas_pointers.dscal(&n, &a, x, &incx)\n+    elif blas_t is float_complex:\n+        blas_pointers.cscal(&n, &a, x, &incx)\n+    else:\n+        blas_pointers.zscal(&n, &a, x, &incx)\n+\n+cdef inline void axpy(int n, blas_t a, blas_t* x, int incx,\n+                      blas_t* y, int incy) nogil:\n+    if blas_t is float:\n+        blas_pointers.saxpy(&n, &a, x, &incx, y, &incy)\n+    elif blas_t is double:\n+        blas_pointers.daxpy(&n, &a, x, &incx, y, &incy)\n+    elif blas_t is float_complex:\n+        blas_pointers.caxpy(&n, &a, x, &incx, y, &incy)\n+    else:\n+        blas_pointers.zaxpy(&n, &a, x, &incx, y, &incy)\n+\n+cdef inline blas_t nrm2(int n, blas_t* x, int incx) nogil:\n+    # since dnrm2 etc are functions, we can't at present call them...\n+    # this is mostly cribbed from netlib's dnrm2\n+    cdef double norm, scale, ssq, absxi\n+    cdef float fnorm, fscale, fssq, fabsxi\n+    if blas_t is float:\n+        if n < 1 or incx < 1:\n+            fnorm = 0\n+        elif n == 1:\n+            fnorm = fabs(x[0])\n+        else:\n+            fscale = 0.0\n+            fssq = 1.0\n+    \n+            for k in range(n):\n+                if index1(x, &incx, k)[0] != 0:\n+                    fabsxi = fabs(index1(x, &incx, k)[0])\n+                    if fscale < fabsxi:\n+                        fssq = 1 + fssq * (fscale/fabsxi)**2\n+                        fscale = fabsxi\n+                    else:\n+                        fssq = fssq + (fabsxi/fscale)**2\n+            fnorm = fscale*sqrt(fssq)\n+        return fnorm\n+    elif blas_t is double:\n+        if n < 1 or incx < 1:\n+            norm = 0\n+        elif n == 1:\n+            norm = fabs(x[0])\n+        else:\n+            scale = 0.0\n+            ssq = 1.0\n+    \n+            for k in range(n):\n+                if index1(x, &incx, k)[0] != 0:\n+                    absxi = fabs(index1(x, &incx, k)[0])\n+                    if scale < absxi:\n+                        ssq = 1 + ssq * (scale/absxi)**2\n+                        scale = absxi\n+                    else:\n+                        ssq = ssq + (absxi/scale)**2\n+            norm = scale*sqrt(ssq)\n+        return norm\n+    elif blas_t is float_complex:\n+        if n < 1 or incx < 1:\n+            fnorm = 0\n+        elif n == 1:\n+            fnorm = hypot(x[0].real, x[0].imag)\n+        else:\n+            fscale = 0.0\n+            fssq = 1.0\n+    \n+            for k in range(n):\n+                if index1(x, &incx, k)[0] != 0:\n+                    fabsxi = hypot(index1(x, &incx, k)[0].real, index1(x, &incx, k)[0].imag)\n+                    if fscale < fabsxi:\n+                        fssq = 1 + fssq * (fscale/fabsxi)**2\n+                        fscale = fabsxi\n+                    else:\n+                        fssq = fssq + (fabsxi/fscale)**2\n+            fnorm = fscale*sqrt(fssq)\n+        return fnorm\n+    else:\n+        if n < 1 or incx < 1:\n+            norm = 0\n+        elif n == 1:\n+            norm = hypot(x[0].real, x[0].imag)\n+        else:\n+            scale = 0.0\n+            ssq = 1.0\n+    \n+            for k in range(n):\n+                if index1(x, &incx, k)[0] != 0:\n+                    absxi = hypot(index1(x, &incx, k)[0].real, index1(x, &incx, k)[0].imag)\n+                    if scale < absxi:\n+                        ssq = 1 + ssq * (scale/absxi)**2\n+                        scale = absxi\n+                    else:\n+                        ssq = ssq + (absxi/scale)**2\n+            norm = scale*sqrt(ssq)\n+        return norm\n+\n+cdef inline void lartg(blas_t* a, blas_t* b, blas_t* c, blas_t* s) nogil:\n+    cdef blas_t g\n+    if blas_t is float:\n+        lapack_pointers.slartg(a, b, c, s, &g)\n+    elif blas_t is double:\n+        lapack_pointers.dlartg(a, b, c, s, &g)\n+    elif blas_t is float_complex:\n+        lapack_pointers.clartg(a, b, <float*>c, s, &g)\n+    else:\n+        lapack_pointers.zlartg(a, b, <double*>c, s, &g)\n+    # make this function more like the BLAS drotg\n+    a[0] = g\n+    b[0] = 0\n+\n+cdef inline void rot(int n, blas_t* x, int incx, blas_t* y, int incy,\n+                     blas_t c, blas_t s) nogil:\n+    if blas_t is float:\n+        blas_pointers.srot(&n, x, &incx, y, &incy, &c, &s) \n+    elif blas_t is double:\n+        blas_pointers.drot(&n, x, &incx, y, &incy, &c, &s)\n+    elif blas_t is float_complex:\n+        lapack_pointers.crot(&n, x, &incx, y, &incy, <float*>&c, &s)\n+    else:\n+        lapack_pointers.zrot(&n, x, &incx, y, &incy, <double*>&c, &s)\n+\n+cdef inline void larfg(int n, blas_t* alpha, blas_t* x, int incx,\n+                       blas_t* tau) nogil:\n+    if blas_t is float:\n+        lapack_pointers.slarfg(&n, alpha, x, &incx, tau)\n+    elif blas_t is double:\n+        lapack_pointers.dlarfg(&n, alpha, x, &incx, tau)\n+    elif blas_t is float_complex:\n+        lapack_pointers.clarfg(&n, alpha, x, &incx, tau)\n+    else:\n+        lapack_pointers.zlarfg(&n, alpha, x, &incx, tau)\n+\n+cdef inline void larf(char* side, int m, int n, blas_t* v, int incv, blas_t tau,\n+    blas_t* c, int ldc, blas_t* work) nogil:\n+    if blas_t is float:\n+        lapack_pointers.slarf(side, &m, &n, v, &incv, &tau, c, &ldc, work)\n+    elif blas_t is double:\n+        lapack_pointers.dlarf(side, &m, &n, v, &incv, &tau, c, &ldc, work)\n+    elif blas_t is float_complex:\n+        lapack_pointers.clarf(side, &m, &n, v, &incv, &tau, c, &ldc, work)\n+    else:\n+        lapack_pointers.zlarf(side, &m, &n, v, &incv, &tau, c, &ldc, work)\n+\n+cdef inline void ger(int m, int n, blas_t alpha, blas_t* x, int incx, blas_t* y,\n+        int incy, blas_t* a, int lda) nogil:\n+    if blas_t is float:\n+        blas_pointers.sger(&m, &n, &alpha, x, &incx, y, &incy, a, &lda)\n+    elif blas_t is double:\n+        blas_pointers.dger(&m, &n, &alpha, x, &incx, y, &incy, a, &lda)\n+    elif blas_t is float_complex:\n+        blas_pointers.cgeru(&m, &n, &alpha, x, &incx, y, &incy, a, &lda)\n+    else:\n+        blas_pointers.zgeru(&m, &n, &alpha, x, &incx, y, &incy, a, &lda)\n+\n+cdef inline void gemv(char* trans, int m, int n, blas_t alpha, blas_t* a,\n+        int lda, blas_t* x, int incx, blas_t beta, blas_t* y, int incy) nogil:\n+    if blas_t is float:\n+        blas_pointers.sgemv(trans, &m, &n, &alpha, a, &lda, x, &incx, &beta,\n+                y, &incy)\n+    elif blas_t is double:\n+        blas_pointers.dgemv(trans, &m, &n, &alpha, a, &lda, x, &incx, &beta,\n+                y, &incy)\n+    elif blas_t is float_complex:\n+        blas_pointers.cgemv(trans, &m, &n, &alpha, a, &lda, x, &incx, &beta,\n+                y, &incy)\n+    else:\n+        blas_pointers.zgemv(trans, &m, &n, &alpha, a, &lda, x, &incx, &beta,\n+                y, &incy)\n+\n+cdef inline void gemm(char* transa, char* transb, int m, int n, int k,\n+        blas_t alpha, blas_t* a, int lda, blas_t* b, int ldb, blas_t beta,\n+        blas_t* c, int ldc) nogil:\n+    if blas_t is float:\n+        blas_pointers.sgemm(transa, transb, &m, &n, &k, &alpha, a, &lda,\n+                b, &ldb, &beta, c, &ldc)\n+    elif blas_t is double:\n+        blas_pointers.dgemm(transa, transb, &m, &n, &k, &alpha, a, &lda,\n+                b, &ldb, &beta, c, &ldc)\n+    elif blas_t is float_complex:\n+        blas_pointers.cgemm(transa, transb, &m, &n, &k, &alpha, a, &lda,\n+                b, &ldb, &beta, c, &ldc)\n+    else:\n+        blas_pointers.zgemm(transa, transb, &m, &n, &k, &alpha, a, &lda,\n+                b, &ldb, &beta, c, &ldc)\n+\n+cdef inline void trmm(char* side, char* uplo, char* transa, char* diag, int m,\n+        int n, blas_t alpha, blas_t* a, int lda, blas_t* b, int ldb) nogil:\n+    if blas_t is float:\n+        blas_pointers.strmm(side, uplo, transa, diag, &m, &n, &alpha, a, &lda,\n+                b, &ldb)\n+    elif blas_t is double:\n+        blas_pointers.dtrmm(side, uplo, transa, diag, &m, &n, &alpha, a, &lda,\n+                b, &ldb)\n+    elif blas_t is float_complex:\n+        blas_pointers.ctrmm(side, uplo, transa, diag, &m, &n, &alpha, a, &lda,\n+                b, &ldb)\n+    else:\n+        blas_pointers.ztrmm(side, uplo, transa, diag, &m, &n, &alpha, a, &lda,\n+                b, &ldb)\n+\n+cdef inline int geqrf(int m, int n, blas_t* a, int lda, blas_t* tau,\n+                      blas_t* work, int lwork) nogil:\n+    cdef int info\n+    if blas_t is float:\n+        lapack_pointers.sgeqrf(&m, &n, a, &lda, tau, work, &lwork, &info)\n+    elif blas_t is double:\n+        lapack_pointers.dgeqrf(&m, &n, a, &lda, tau, work, &lwork, &info)\n+    elif blas_t is float_complex:\n+        lapack_pointers.cgeqrf(&m, &n, a, &lda, tau, work, &lwork, &info)\n+    else:\n+        lapack_pointers.zgeqrf(&m, &n, a, &lda, tau, work, &lwork, &info)\n+    return info\n+\n+cdef inline int ormqr(char* side, char* trans, int m, int n, int k, blas_t* a,\n+    int lda, blas_t* tau, blas_t* c, int ldc, blas_t* work, int lwork) nogil:\n+    cdef int info = 0 \n+    if blas_t is float:\n+        lapack_pointers.sormqr(side, trans, &m, &n, &k, a, &lda, tau, c, &ldc,\n+                work, &lwork, &info)\n+    elif blas_t is double:\n+        lapack_pointers.dormqr(side, trans, &m, &n, &k, a, &lda, tau, c, &ldc,\n+                work, &lwork, &info)\n+    elif blas_t is float_complex:\n+        lapack_pointers.cunmqr(side, trans, &m, &n, &k, a, &lda, tau, c, &ldc,\n+                work, &lwork, &info)\n+    else:\n+        lapack_pointers.zunmqr(side, trans, &m, &n, &k, a, &lda, tau, c, &ldc,\n+                work, &lwork, &info)\n+    return info\n+\n+#------------------------------------------------------------------------------\n+# Utility routines\n+#------------------------------------------------------------------------------\n+\n+cdef cnp.ndarray PyArray_FromArraySafe(cnp.ndarray arr, void* newtype, int flags):\n+    \"\"\"In Numpy 1.5.1, Use of the NPY_F_CONTIGUOUS flag is broken when used\n+    with a 1D input array.  The work around is to pass NPY_C_CONTIGUOUS since\n+    for 1D arrays these are equivalent.  This is Numpy's gh-2287, fixed in\n+    9b8ff38.\n+\n+    FIXME: Is it worth only applying this for numpy 1.5.1? \n+    \"\"\"\n+    if arr.ndim == 1 and (flags & cnp.NPY_F_CONTIGUOUS):\n+        flags = flags & ~cnp.NPY_F_CONTIGUOUS\n+        flags |= cnp.NPY_C_CONTIGUOUS\n+    return PyArray_FromArray(arr, newtype, flags)\n+\n+cdef void blas_t_conj(int n, blas_t* x, int* xs) nogil:\n+    cdef int j\n+    if blas_t is float_complex or blas_t is double_complex:\n+        for j in range(n):\n+            index1(x, xs, j)[0] = index1(x, xs, j)[0].conjugate()\n+\n+cdef void blas_t_2d_conj(int m, int n, blas_t* x, int* xs) nogil:\n+    cdef int i, j\n+    if blas_t is float_complex or blas_t is double_complex:\n+        for i in range(m):\n+            for j in range(n):\n+                index2(x, xs, i, j)[0] = index2(x, xs, i, j)[0].conjugate()\n+\n+cdef blas_t blas_t_sqrt(blas_t x) nogil:\n+    if blas_t is float:\n+        return sqrt(x)\n+    elif blas_t is double:\n+        return sqrt(x)\n+    elif blas_t is float_complex:\n+        return <float_complex>sqrt(<double>((<float*>&x)[0]))\n+    else:\n+        return sqrt((<double*>&x)[0])\n+\n+cdef bint blas_t_less_than(blas_t x, blas_t y) nogil:\n+    if blas_t is float or blas_t is double:\n+        return x < y\n+    else:\n+        return x.real < y.real\n+\n+cdef int to_lwork(blas_t a, blas_t b) nogil:\n+    cdef int ai, bi\n+    if blas_t is float or blas_t is double:\n+        ai = <int>a\n+        bi = <int>b\n+    elif blas_t is float_complex:\n+        ai = <int>((<float*>&a)[0])\n+        bi = <int>((<float*>&b)[0])\n+    elif blas_t is double_complex:\n+        ai = <int>((<double*>&a)[0])\n+        bi = <int>((<double*>&b)[0])\n+    return max(ai, bi)\n+\n+#------------------------------------------------------------------------------\n+# QR update routines start here.\n+#------------------------------------------------------------------------------\n+\n+cdef bint reorthx(int m, int n, blas_t* q, int* qs, bint qisF, int j, blas_t* u, blas_t* s) nogil:\n+    # U should be all zeros on entry.\n+    cdef blas_t unorm, snorm, wnorm, wpnorm, sigma_max, sigma_min, rc\n+    cdef char* T = 'T'\n+    cdef char* N = 'N'\n+    cdef char* C = 'C'\n+    cdef int ss = 1\n+    cdef blas_t inv_root2 = <blas_t>sqrt(2)\n+\n+    # u starts out as the jth basis vector.\n+    u[j] = 1\n+\n+    # s = Q.T.dot(u) = jth row of Q.\n+    copy(n, row(q, qs, j), qs[1], s, 1)\n+    blas_t_conj(n, s, &ss)\n+\n+    # make u be the part of u that is not in span(q)\n+    # i.e. u -= q.dot(s)\n+    if qisF:\n+        gemv(N, m, n, -1, q, qs[1], s, 1, 1, u, 1)\n+    else:\n+        gemv(T, n, m, -1, q, n, s, 1, 1, u, 1)\n+    wnorm = nrm2(m, u, 1)\n+\n+    if blas_t_less_than(inv_root2, wnorm):\n+        scal(m, 1/wnorm, u, 1)\n+        s[n] = wnorm\n+        return True\n+\n+    # if the above check failed, try one reorthogonalization\n+    if qisF:\n+        if blas_t is float or blas_t is double:\n+            gemv(T, m, n, 1, q, qs[1], u, 1, 0, s+n, 1)\n+        else:\n+            gemv(C, m, n, 1, q, qs[1], u, 1, 0, s+n, 1)\n+        gemv(N, m, n, -1, q, qs[1], s+n, 1, 1, u, 1)\n+    else:\n+        if blas_t is float or blas_t is double:\n+            gemv(N, n, m, 1, q, n, u, 1, 0, s+n, 1)\n+        else:\n+            blas_t_conj(m, u, &ss)\n+            gemv(N, n, m, 1, q, n, u, 1, 0, s+n, 1)\n+            blas_t_conj(m, u, &ss)\n+            blas_t_conj(n, s+n, &ss)\n+        gemv(T, n, m, -1, q, n, s+n, 1, 1, u, 1)\n+        \n+    wpnorm = nrm2(m, u, 1) \n+    \n+    if blas_t_less_than(wpnorm, wnorm/inv_root2): # u lies in span(q) \n+        scal(m, 0, u, 1)\n+        axpy(n, 1, s, 1, s+n, 1)\n+        s[n] = 0\n+        return False\n+    scal(m, 1/wpnorm, u, 1)\n+    axpy(n, 1, s, 1, s+n, 1)\n+    s[n] = wpnorm\n+    return True\n+\n+cdef int thin_qr_row_delete(int m, int n, blas_t* q, int* qs, bint qisF, blas_t* r, int* rs, int k, int p_eco, int p_full) nogil:\n+    cdef int i, j, argmin_row_norm \n+    cdef size_t usize = (m + 3*n + 1) * sizeof(blas_t)\n+    cdef blas_t* s\n+    cdef blas_t* u\n+    cdef blas_t* s1\n+    cdef int us[1]\n+    cdef int ss[1]\n+    cdef blas_t c, sn, min_row_norm, row_norm\n+\n+    u = <blas_t*>libc.stdlib.malloc(usize)\n+    if not u:\n+        return MEMORY_ERROR\n+    s = u + m\n+    ss[0] = 1\n+    ss[1] = 0\n+    us[0] = 1\n+    us[1] = 0\n+\n+    for i in range(p_eco):\n+        memset(u, 0, usize)\n+        # permute q such that row k is the last row.\n+        if k != m-1:\n+            for j in range(k, m-1):\n+                swap(n, row(q, qs, j), qs[1], row(q, qs, j+1), qs[1])\n+\n+        if not reorthx(m, n, q, qs, qisF, m-1, u, s):\n+            # if we get here it means that this basis vector lies in span(q).\n+            # we want to use s[:n+1] but we need a vector into null(q)\n+            # find the row of q with the smallest norm and try that. (Daniel, p785)\n+            min_row_norm = nrm2(n, row(q, qs, 0), qs[1])\n+            argmin_row_norm = 0\n+            for j in range(1, m):\n+                row_norm = nrm2(n, row(q, qs, j), qs[1])\n+                if blas_t_less_than(row_norm, min_row_norm):\n+                    min_row_norm = row_norm\n+                    argmin_row_norm = j\n+            memset(u, 0, m*sizeof(blas_t))\n+            if not reorthx(m, n, q, qs, qisF, argmin_row_norm, u, s+n+1):\n+                # failed, quit.\n+                libc.stdlib.free(u)\n+                return 0\n+\n+        memset(s+2*n, 0, n*sizeof(blas_t))\n+\n+        # what happens here...\n+        for j in range(n-1, -1, -1):\n+            lartg(index1(s, ss, n), index1(s, ss, j), &c, &sn)\n+            rot(n-j, index1(s+2*n, ss, j), ss[0], index2(r, rs,j, j), rs[1], c, sn)\n+            rot(m-1, u, us[0], col(q, qs, j), qs[0], c, sn.conjugate())\n+        m -= 1\n+\n+    libc.stdlib.free(u)\n+\n+    if p_full:\n+        qr_block_row_delete(m, n, q, qs, r, rs, k, p_full)\n+    return 1\n+\n+cdef void qr_block_row_delete(int m, int n, blas_t* q, int* qs,\n+                              blas_t* r, int* rs, int k, int p) nogil:\n+    cdef int i, j\n+    cdef blas_t c,s\n+    cdef blas_t* W\n+    cdef int* ws\n+\n+    if k != 0:\n+        for j in range(k, 0, -1):\n+            swap(m, row(q, qs, j+p-1), qs[1], row(q, qs, j-1), qs[1])\n+    \n+    # W is the block of rows to be removed from q, has shape, (p,m)\n+    W = q\n+    ws = qs\n+\n+    for j in range(p):\n+        blas_t_conj(m, row(W, ws, j), &ws[1])\n+\n+    for i in range(p):\n+        for j in range(m-2, i-1, -1):\n+            lartg(index2(W, ws, i, j), index2(W, ws, i, j+1), &c, &s)\n+            \n+            # update W\n+            if i+1 < p:\n+                rot(p-i-1, index2(W, ws, i+1, j), ws[0],\n+                    index2(W, ws, i+1, j+1), ws[0], c, s)\n+\n+            # update r if there is a nonzero row.\n+            if j-i < n:\n+                rot(n-j+i, index2(r, rs, j, j-i), rs[1],\n+                    index2(r, rs, j+1, j-i), rs[1], c, s)\n+\n+            # update q\n+            rot(m-p, index2(q, qs, p, j), qs[0], index2(q, qs, p, j+1), qs[0],\n+                c, s.conjugate())\n+\n+cdef void qr_col_delete(int m, int o, int n, blas_t* q, int* qs, blas_t* r,\n+                        int* rs, int k) nogil:\n+    \"\"\"\n+        Here we support both full and economic decomposition, q is (m,o), and r\n+        is (o, n).\n+    \"\"\"\n+    cdef int j\n+    cdef int limit = min(o, n)\n+\n+    for j in range(k, n-1):\n+        copy(limit, col(r, rs, j+1), rs[0], col(r, rs, j), rs[0])\n+\n+    hessenberg_qr(m, n-1, q, qs, r, rs, k)\n+\n+cdef int qr_block_col_delete(int m, int o, int n, blas_t* q, int* qs,\n+                              blas_t* r, int* rs, int k, int p) nogil:\n+    \"\"\"\n+        Here we support both full and economic decomposition, q is (m,o), and r\n+        is (o, n).\n+    \"\"\"\n+    cdef int j\n+    cdef int limit = min(o, n)\n+    cdef blas_t* work\n+    cdef int worksize = max(m, n)\n+\n+    work = <blas_t*>libc.stdlib.malloc(worksize*sizeof(blas_t))\n+    if not work:\n+        return MEMORY_ERROR\n+\n+    # move the columns to removed to the end\n+    for j in range(k, n-p):\n+        copy(limit, col(r, rs, j+p), rs[0], col(r, rs, j), rs[0])\n+\n+    p_subdiag_qr(m, o, n-p, q, qs, r, rs, k, p, work)\n+\n+    libc.stdlib.free(work)\n+    return 0\n+\n+cdef void thin_qr_row_insert(int m, int n, blas_t* q, int* qs, blas_t* r,\n+        int* rs, blas_t* u, int* us, int k) nogil:\n+    cdef int j\n+    cdef blas_t c, s\n+\n+    for j in range(n):\n+        lartg(index2(r, rs, j, j), index1(u, us, j), &c, &s)\n+        rot(n-j-1, index2(r, rs, j, j+1), rs[1], index1(u, us, j+1), us[0],\n+                c, s)\n+        rot(m, col(q, qs, j), qs[0], col(q, qs, n), qs[0], c, s.conjugate())\n+\n+    # permute q\n+    for j in range(m-1, k, -1):\n+        swap(n, row(q, qs, j), qs[1], row(q, qs, j-1), qs[1])\n+\n+cdef void qr_row_insert(int m, int n, blas_t* q, int* qs, blas_t* r, int* rs,\n+                        int k) nogil:\n+    cdef int j\n+    cdef blas_t c, s\n+    cdef int limit = min(m-1, n)\n+\n+    for j in range(limit):\n+        lartg(index2(r, rs, j, j), index2(r, rs, m-1, j), &c, &s)\n+        rot(n-j-1, index2(r, rs, j, j+1), rs[1], index2(r, rs, m-1, j+1), rs[1],\n+                c, s)\n+        rot(m, col(q, qs, j), qs[0], col(q, qs, m-1), qs[0], c, s.conjugate())\n+\n+    # permute q \n+    for j in range(m-1, k, -1):\n+        swap(m, row(q, qs, j), qs[1], row(q, qs, j-1), qs[1])\n+\n+cdef int thin_qr_block_row_insert(int m, int n, blas_t* q, int* qs, blas_t* r,\n+        int* rs, blas_t* u, int* us, int k, int p) nogil:\n+    # as below this should someday call lapack's xtpqrt.\n+    cdef int j\n+    cdef blas_t rjj, tau\n+    cdef blas_t* work\n+    cdef char* T = 'T'\n+    cdef char* N = 'N'\n+    cdef size_t worksize = m * sizeof(blas_t)\n+\n+    work = <blas_t*>libc.stdlib.malloc(worksize)\n+    if not work:\n+        return MEMORY_ERROR\n+   \n+    # possible FIX\n+    # as this is written it requires F order q, r, and u.  But thats not \n+    # strictly necessary. C order should also work too with a little fiddling.\n+    for j in range(n):\n+        rjj = index2(r, rs, j, j)[0]\n+        larfg(p+1, &rjj, col(u, us, j), us[0], &tau)\n+\n+        # here we apply the reflector by hand instead of calling larf\n+        # since we need to apply it to a stack of r atop u, and these\n+        # are separate.  This also permits the reflector to always be\n+        # p+1 long, rather than having a max of n+p.\n+        copy(n-j, index2(r, rs, j, j+1), rs[1], work, 1)\n+        blas_t_conj(p, col(u, us, j), &us[0])\n+        gemv(T, p, n-j-1, 1, index2(u, us, 0, j+1), us[1], col(u, us, j), us[0],\n+             1, work, 1)\n+        blas_t_conj(p, col(u, us, j), &us[0])\n+        ger(p, n-j-1, -tau.conjugate(), col(u, us, j), us[0], work, 1,\n+            index2(u, us, 0, j+1), us[1])\n+        axpy(n-j-1, -tau.conjugate(), work, 1, index2(r, rs, j, j+1), rs[1])\n+        index2(r, rs, j, j)[0] = rjj\n+        \n+        # now apply this reflector to q \n+        copy(m, col(q, qs, j), qs[0], work, 1)\n+        gemv(N, m, p, 1, index2(q, qs, 0, n), qs[1], col(u, us, j), us[0],\n+             1, work, 1)\n+        blas_t_conj(p, col(u, us, j), &us[0])\n+        ger(m, p, -tau, work, 1, col(u, us, j), us[0],\n+            index2(q, qs, 0, n), qs[1])\n+        axpy(m, -tau, work, 1, col(q, qs, j), qs[0])\n+\n+    # permute the rows of q, work columnwise, since q is fortran order\n+    if k != m-p:\n+        for j in range(n):\n+            copy(m-k-p, index2(q, qs, k, j), qs[0], work, 1)\n+            copy(p, index2(q, qs, m-p, j), qs[0], index2(q, qs, k, j), qs[0])\n+            copy(m-k-p, work, 1, index2(q, qs, k+p, j), qs[0])\n+\n+    libc.stdlib.free(work)\n+\n+cdef int qr_block_row_insert(int m, int n, blas_t* q, int* qs,\n+                              blas_t* r, int* rs, int k, int p) nogil:\n+    # this should someday call lapack's xtpqrt (requires lapack >= 3.4\n+    # released nov 11). RHEL6's atlas doesn't seem to have it.\n+    # On input this looks something like this:\n+    # q = x x x x 0 0 0  r = x x x\n+    #     x x x x 0 0 0      0 x x\n+    #     x x x x 0 0 0      0 0 x\n+    #     x x x x 0 0 0      0 0 0\n+    #     0 0 0 0 1 0 0      * * *\n+    #     0 0 0 0 0 1 0      * * *\n+    #     0 0 0 0 0 0 1      * * *\n+    #\n+    # The method will be to apply a series of reflectors to re triangularize r.\n+    # followed by permuting the rows of q to put the new rows in the requested\n+    # position.\n+    cdef int j, hlen\n+    cdef blas_t rjj, tau\n+    cdef blas_t* work\n+    cdef char* sideL = 'L'\n+    cdef char* sideR = 'R'\n+    # for tall or sqr + rows should be n. for fat + rows should be new m\n+    cdef int limit = min(m, n)\n+\n+    work = <blas_t*>libc.stdlib.malloc(max(m,n)*sizeof(blas_t))\n+    if not work:\n+        return MEMORY_ERROR\n+\n+    for j in range(limit):\n+        rjj = index2(r, rs, j, j)[0]\n+        hlen = m-j\n+        larfg(hlen, &rjj, index2(r, rs, j+1, j), rs[0], &tau)\n+        index2(r, rs, j, j)[0] = 1\n+        if j+1 < n:\n+            larf(sideL, hlen, n-j-1, index2(r, rs, j, j), rs[0],\n+                    tau.conjugate(), index2(r, rs, j, j+1), rs[1], work)\n+        larf(sideR, m, hlen, index2(r, rs, j, j), rs[0], tau,\n+                index2(q, qs, 0, j), qs[1], work)\n+        memset(index2(r, rs, j, j), 0, hlen*sizeof(blas_t))\n+        index2(r, rs, j, j)[0] = rjj\n+\n+    # permute the rows., work columnwise, since q is fortran order\n+    if k != m-p:\n+        for j in range(m):\n+            copy(m-k-p, index2(q, qs, k, j), qs[0], work, 1)\n+            copy(p, index2(q, qs, m-p, j), qs[0], index2(q, qs, k, j), qs[0])\n+            copy(m-k-p, work, 1, index2(q, qs, k+p, j), qs[0])\n+\n+    libc.stdlib.free(work)\n+    return 0\n+\n+cdef void thin_qr_col_insert(int m, int n, blas_t* q, int* qs, blas_t* r,\n+        int* rs, blas_t* u, int* us, int k, int p_eco, int p_full,\n+        blas_t* rcond) nogil:\n+    # here q and r will always be fortran ordered since we have to allocate them\n+    cdef int i, j, info\n+    cdef blas_t c, sn \n+    cdef blas_t rc\n+    cdef blas_t* s\n+    cdef char* N = 'N'\n+    cdef char* T = 'T'\n+    cdef char* C = 'C'\n+    cdef char* TC\n+    \n+    if blas_t is float or blas_t is double:\n+        TC = T\n+    else:\n+        TC = C\n+    \n+    # on entry, Q and R have both been increased in size, Q via the appending\n+    # columns of zeros, and R by the addition of both columns and rows of\n+    # zeros.  In R, the new columns are located from column k to k + p. and \n+    # the new rows are at the bottom.  m and n refer to the size of the\n+    # original system, not the new system.\n+    \n+    s = <blas_t*>libc.stdlib.malloc(2*(n+p_eco)*sizeof(blas_t))\n+\n+    for j in range(p_eco):\n+        rc = rcond[0]\n+        info = reorth(m, n+j, q, qs, True, col(u, us, j), us, s,  &rc)\n+        #libc.stdio.printf('info = %d, rcond = %g\\n', info, rc)\n+        copy(m, col(u, us, j), us[0], col(q, qs, n+j), qs[0])\n+        copy(n+j+1, s, 1, col(r, rs, k+j), rs[0]) \n+\n+        for i in range(n-2+1, k-1, -1):\n+            lartg(index2(r, rs, i+j, k+j), index2(r, rs, i+j+1, k+j), &c, &sn)\n+            rot(n-i, index2(r, rs, i+j, i+p_eco+p_full), rs[1],\n+                index2(r, rs, i+j+1, i+p_eco+p_full), rs[1], c, sn)\n+            rot(m, col(q, qs, i+j), qs[0], col(q, qs, i+j+1), qs[0],\n+                c, sn.conjugate())\n+    libc.stdlib.free(s)\n+    \n+    if p_full > 0:\n+        # if this is true, we have ensured the u is also F contiguous.\n+        gemm(TC, N, m, p_full, m, 1, q, qs[1], col(u, us, p_eco), us[1], 0,\n+             col(r, rs, k+p_eco), rs[1]) \n+        qr_block_col_insert(m, n+p_eco+p_full, q, qs, r, rs, k+p_eco, p_full)\n+\n+cdef void qr_col_insert(int m, int n, blas_t* q, int* qs, blas_t* r, int* rs,\n+                        int k) nogil:\n+    cdef int j\n+    cdef blas_t c, s, temp, tau\n+    cdef blas_t* work\n+    \n+    for j in range(m-2, k-1, -1):\n+        lartg(index2(r, rs, j, k), index2(r, rs, j+1, k), &c, &s)\n+\n+        # update r if j is a nonzero row\n+        if j+1 < n:\n+            rot(n-j-1, index2(r, rs, j, j+1), rs[1],\n+                    index2(r, rs, j+1, j+1), rs[1], c, s)\n+\n+        # update the columns of q\n+        rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0], c, s.conjugate())\n+\n+cdef int qr_block_col_insert(int m, int n, blas_t* q, int* qs,\n+                              blas_t* r, int* rs, int k, int p) nogil:\n+    cdef int i, j\n+    cdef blas_t c, s\n+    cdef blas_t* tau = NULL\n+    cdef blas_t* work = NULL\n+    cdef int info, lwork\n+    cdef char* side = 'R'\n+    cdef char* trans = 'N'\n+\n+    if m >= n:\n+        # if m > n, r looks like this.\n+        # x x x x x x x x x x\n+        #   x x x x x x x x x\n+        #     x x x x x x x x\n+        #       x x x x x x x\n+        #       x x x   x x x\n+        #       x x x     x x\n+        #       x x x       x\n+        #       x x x \n+        #       x x x \n+        #       x x x \n+        #       x x x \n+        #       x x x \n+        #\n+        # First zero the lower part of the new columns using a qr.  \n+\n+        # query the workspace, \n+        info = geqrf(m-n+p, p, index2(r, rs, n-p, k), rs[1], tau, &c, -1)\n+        info = ormqr(side, trans, m, m-(n-p), p, index2(r, rs, n-p, k), rs[1],\n+                tau, index2(q, qs, 0, n-p), qs[1], &s, -1)\n+\n+        # we're only doing one allocation, so use the larger\n+        lwork = to_lwork(c, s)\n+\n+        # allocate the workspace + tau\n+        work = <blas_t*>libc.stdlib.malloc((lwork+min(m-n+p, p))*sizeof(blas_t))\n+        if not work:\n+            return MEMORY_ERROR\n+        tau = work + lwork\n+\n+        # qr\n+        info = geqrf(m-n+p, p, index2(r, rs, n-p, k), rs[1], tau, work, lwork)\n+        if info < 0: \n+            return libc.stdlib.abs(info)\n+\n+        # apply the Q from this small qr to the last (m-(n-p)) columns of q.\n+        info = ormqr(side, trans, m, m-(n-p), p, index2(r, rs, n-p, k), rs[1],\n+                tau, index2(q, qs, 0, n-p), qs[1], work, lwork)\n+        if info < 0:\n+            return info\n+\n+        libc.stdlib.free(work)\n+\n+        # zero the reflectors since we're done with them\n+        # memset can be used here, since r is always fortan order\n+        for j in range(p):\n+            memset(index2(r, rs, n-p+1+j, k+j), 0, (m-(n-p+1+j))*sizeof(blas_t))\n+\n+        # now we have something that looks like \n+        # x x x x x x x x x x\n+        #   x x x x x x x x x\n+        #     x x x x x x x x\n+        #       x x x x x x x\n+        #       x x x   x x x\n+        #       x x x     x x\n+        #       x x x       x\n+        #       x x x         \n+        #       0 x x        \n+        #       0 0 x         \n+        #       0 0 0 \n+        #       0 0 0 \n+        #\n+        # and the rest of the columns need to be eliminated using rotations.\n+\n+        for i in range(p):\n+            for j in range(n-p+i-1, k+i-1, -1):\n+                lartg(index2(r, rs, j, k+i), index2(r, rs, j+1, k+i), &c, &s)\n+                if j+1 < n:\n+                    rot(n-k-i-1, index2(r, rs, j, k+i+1), rs[1],\n+                            index2(r, rs, j+1, k+i+1), rs[1], c, s)\n+                rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0],\n+                        c, s.conjugate())\n+    else: \n+        # this case we can only uses givens rotations.\n+        for i in range(p):\n+            for j in range(m-2, k+i-1, -1):\n+                lartg(index2(r, rs, j, k+i), index2(r, rs, j+1, k+i), &c, &s)\n+                if j+1 < n:\n+                    rot(n-k-i-1, index2(r, rs, j, k+i+1), rs[1],\n+                            index2(r, rs, j+1, k+i+1), rs[1], c, s)\n+                rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0],\n+                        c, s.conjugate())\n+    return 0\n+\n+cdef void thin_qr_rank_1_update(int m, int n, blas_t* q, int* qs, bint qisF, \n+    blas_t* r, int* rs, blas_t* u, int* us, blas_t* v, int* vs, blas_t* s,\n+    int* ss) nogil:\n+    \"\"\"Assume that q is (M,N) and either C or F contiguous, r is (N,N), u is M,\n+       and V is N.  s is a 2*n work array.\n+    \"\"\"\n+    cdef int j, info\n+    cdef blas_t c, sn, rlast, t, rcond = 0.0\n+    \n+    info = reorth(m, n, q, qs, qisF, u, us, s, &rcond)\n+\n+    # reduce s with givens, using u as the n+1 column of q\n+    # do the first one since the rots will be different.\n+    lartg(index1(s, ss, n-1), index1(s, ss, n), &c, &sn)\n+    t = index2(r, rs, n-1, n-1)[0] \n+    rlast = -t * sn.conjugate()\n+    index2(r, rs, n-1, n-1)[0] = t * c\n+    rot(m, col(q, qs, n-1), qs[0], u, us[0], c, sn.conjugate())\n+\n+    for j in range(n-2, -1, -1):\n+        lartg(index1(s, ss, j), index1(s, ss, j+1), &c, &sn)\n+        rot(n-j, index2(r, rs, j, j), rs[1],\n+                index2(r, rs, j+1, j), rs[1], c, sn)\n+        rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0], c, sn.conjugate())\n+\n+    # add v to the first row of r\n+    blas_t_conj(n, v, vs)\n+    axpy(n, s[0],  v, vs[0], row(r, rs, 0), rs[1])\n+\n+    # now r is upper hessenberg with the only value in the last row stored in \n+    # rlast (This is very similar to hessenberg_qr below, but this loop ends \n+    # at n-1 instead of n)\n+    for j in range(n-1):\n+        lartg(index2(r, rs, j, j), index2(r, rs, j+1, j), &c, &sn)\n+        rot(n-j-1, index2(r, rs, j, j+1), rs[1],\n+                index2(r, rs, j+1, j+1), rs[1], c, sn)\n+        rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0], c, sn.conjugate())\n+\n+    # handle the extra value in rlast\n+    lartg(index2(r, rs, n-1, n-1), &rlast, &c, &sn)\n+    rot(m, col(q, qs, n-1), qs[0], u, us[0], c, sn.conjugate())\n+\n+cdef void thin_qr_rank_p_update(int m, int n, int p, blas_t* q, int* qs,\n+    bint qisF, blas_t* r, int* rs, blas_t* u, int* us, blas_t* v, int* vs,\n+    blas_t* s, int* ss) nogil:\n+    \"\"\"Assume that q is (M,N) and either C or F contiguous, r is (N,N), u is\n+       (M,p) and V is (N,p).  s is a 2*n work array.\n+    \"\"\"\n+    cdef int j\n+\n+    for j in range(p):\n+        thin_qr_rank_1_update(m, n, q, qs, qisF, r, rs, col(u, us, j), us,\n+                              col(v, vs, j), vs, s, ss)\n+  \n+cdef void qr_rank_1_update(int m, int n, blas_t* q, int* qs, blas_t* r, int* rs,\n+                           blas_t* u, int* us, blas_t* v, int* vs) nogil:\n+    \"\"\" here we will assume that the u = Q.T.dot(u) and not the bare u.\n+        if A is MxN then q is MxM, r is MxN, u is M and v is N.\n+        e.g. currently assuming full matrices.\n+    \"\"\"\n+    cdef int j\n+    cdef blas_t c, s\n+\n+    # The technique here is to reduce u to a series of givens rotations followed\n+    # by a scalar e.g. [u1,u2,u3] --> [u,0,0].  Applying these rotations to r as\n+    # we go.  Then we will have the update be adding v scaled by the remainder\n+    # of u to the first row of r, which will be upper hessenberg due to the\n+    # givens applied to reduce u. We then reduce the upper hessenberg r to upper\n+    # triangular.\n+\n+    for j in range(m-2, -1, -1):\n+        lartg(index1(u, us, j), index1(u, us, j+1), &c, &s)\n+\n+        # update jth and (j+1)th rows of r.\n+        if n-j > 0:\n+            rot(n-j, index2(r, rs, j, j), rs[1], index2(r, rs, j+1, j), rs[1], c, s)\n+\n+        # update jth and (j+1)th cols of q.\n+        rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0], c, s.conjugate())\n+\n+    # add v to the first row\n+    blas_t_conj(n, v, vs)\n+    axpy(n, u[0],  v, vs[0], row(r, rs, 0), rs[1])\n+    \n+    # return to q, r form\n+    hessenberg_qr(m, n, q, qs, r, rs, 0)\n+    # no return, return q, r from python driver.\n+\n+cdef int qr_rank_p_update(int m, int n, int p, blas_t* q, int* qs, blas_t* r,\n+                        int* rs, blas_t* u, int* us, blas_t* v, int* vs) nogil:\n+    cdef int i, j\n+    cdef blas_t c, s\n+    cdef blas_t* tau = NULL\n+    cdef blas_t* work = NULL\n+    cdef int info, lwork\n+    cdef char* sideR = 'R'\n+    cdef char* sideL = 'L'\n+    cdef char* uplo = 'U'\n+    cdef char* trans = 'N'\n+    cdef char* diag = 'N'\n+\n+    if m > n:\n+        # query the workspace\n+        # below p_subdiag_qr will need workspace of size m, which is the\n+        # minimum, ormqr will also require.\n+        info = geqrf(m-n, p, index2(u, us, n, 0), us[1], tau, &c, -1)\n+        if info < 0:\n+            return libc.stdlib.abs(info)\n+        info = ormqr(sideR, trans, m, m-n, p, index2(u, us, n, 0), us[1], tau,\n+                index2(q, qs, 0, n), qs[1], &s, -1)\n+        if info < 0:\n+            return info\n+ \n+        # we're only doing one allocation, so use the larger\n+        lwork = to_lwork(c, s)\n+\n+        # allocate the workspace + tau\n+        work = <blas_t*>libc.stdlib.malloc((lwork+min(m-n, p))*sizeof(blas_t))\n+        if not work:\n+            return MEMORY_ERROR\n+        tau = work + lwork\n+\n+        # qr\n+        info = geqrf(m-n, p, index2(u, us, n, 0), us[1], tau, work, lwork)\n+        if info < 0:\n+            libc.stdlib.free(work)\n+            return libc.stdlib.abs(info)\n+\n+        # apply the Q from this small qr to the last (m-n) columns of q.\n+        info = ormqr(sideR, trans, m, m-n, p, index2(u, us, n, 0), us[1], tau,\n+                index2(q, qs, 0, n), qs[1], work, lwork)\n+        if info < 0:\n+            libc.stdlib.free(work)\n+            return info\n+\n+        # reduce u the rest of the way to upper triangular using givens.\n+        for i in range(p):\n+            for j in range(n+i-1, i-1, -1):\n+                lartg(index2(u, us, j, i), index2(u, us, j+1, i), &c, &s)\n+                if p-i-1:\n+                    rot(p-i-1, index2(u, us, j, i+1), us[1],\n+                            index2(u, us, j+1, i+1), us[1], c, s)\n+                rot(n, row(r, rs, j), rs[1], row(r, rs, j+1), rs[1], c, s)\n+                rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0],\n+                        c, s.conjugate())\n+\n+    else: # m == n or m < n\n+        # reduce u to upper triangular using givens.\n+        for i in range(p):\n+            for j in range(m-2, i-1, -1):\n+                lartg(index2(u, us, j, i), index2(u, us, j+1, i), &c, &s)\n+                if p-i-1:\n+                    rot(p-i-1, index2(u, us, j, i+1), us[1],\n+                            index2(u, us, j+1, i+1), us[1], c, s)\n+                rot(n, row(r, rs, j), rs[1], row(r, rs, j+1), rs[1], c, s)\n+                rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0],\n+                        c, s.conjugate())\n+\n+        # allocate workspace\n+        work = <blas_t*>libc.stdlib.malloc(n*sizeof(blas_t))\n+        if not work:\n+            return MEMORY_ERROR\n+\n+    # now form UV**H and add it to R.\n+    # This won't fill in any more of R than we have already.\n+    blas_t_2d_conj(p, n, v, vs)\n+    trmm(sideL, uplo, trans, diag, p, n, 1, u, us[1], v, vs[1])\n+\n+    # (should this be n, p length adds instead since these are fortan contig?)\n+    for j in range(p):\n+        axpy(n, 1, row(v, vs, j), vs[1], row(r, rs, j), rs[1])\n+\n+    # now r has p subdiagonals, eliminate them with reflectors.\n+    p_subdiag_qr(m, m, n, q, qs, r, rs, 0, p, work)\n+\n+    libc.stdlib.free(work)\n+    return 0\n+\n+cdef void hessenberg_qr(int m, int n, blas_t* q, int* qs, blas_t* r, int* rs,\n+                        int k) nogil:\n+    \"\"\"Reduce an upper hessenberg matrix r, to upper triangluar, starting in\n+       row j.  Apply these transformation to q as well. Both full and economic\n+       decompositions are supported here. \n+    \"\"\"\n+    cdef int j\n+    cdef blas_t c, s\n+    cdef int limit = min(m-1, n)\n+\n+    for j in range(k, limit):\n+        lartg(index2(r, rs, j, j), index2(r, rs, j+1, j), &c, &s)\n+\n+        # update the rest of r\n+        if j+1 < m:\n+            rot(n-j-1, index2(r, rs, j, j+1), rs[1],\n+                    index2(r, rs, j+1, j+1), rs[1], c, s) \n+\n+        # update q\n+        rot(m, col(q, qs, j), qs[0], col(q, qs, j+1), qs[0], c, s.conjugate())\n+\n+cdef void p_subdiag_qr(int m, int o, int n, blas_t* q, int* qs, blas_t* r, int* rs,\n+                       int k, int p, blas_t* work) nogil:\n+    \"\"\" Reduce a matrix r to upper triangular form by eliminating the lower p \n+        subdiagionals using reflectors. Both full and economic decompositions\n+        are supported here.  In either case, q is (m,o) and r is (o,n)\n+        \n+        q and r must be fortran order here, with work at least max(m,n) long.\n+    \"\"\"\n+    cdef int j\n+    cdef int last\n+    cdef blas_t tau\n+    cdef blas_t rjj \n+    cdef int limit = min(m-1, n)\n+    cdef char* sideR = 'R'\n+    cdef char* sideL = 'L'\n+\n+    # R now has p subdiagonal values to be removed starting from col k.\n+    for j in range(k, limit):\n+        # length of the reflector\n+        last = min(p+1, o-j)\n+        rjj = index2(r, rs, j, j)[0]\n+        larfg(last, &rjj, index2(r, rs, j+1, j), rs[0], &tau)\n+        index2(r, rs, j, j)[0] = 1\n+\n+        # apply the reflector to r if necessary\n+        if j+1 < n:\n+            larf(sideL, last, n-j-1, index2(r, rs, j, j), rs[0],\n+                    tau.conjugate(), index2(r, rs, j, j+1), rs[1], work)\n+\n+        # apply the reflector to q\n+        larf(sideR, m, last, index2(r, rs, j, j), rs[0], tau,\n+                index2(q, qs, 0, j), qs[1], work)\n+\n+        # rezero the householder vector we no longer need.\n+        memset(index2(r, rs, j+1, j), 0, (last-1)*sizeof(blas_t))\n+\n+        # restore the rjj element\n+        index2(r, rs, j, j)[0] = rjj\n+\n+def _reorth(cnp.ndarray q, cnp.ndarray u, rcond):\n+    cdef cnp.ndarray s\n+    cdef double rc_d\n+    cdef float rc_f\n+    cdef int m, n\n+    cdef void* qp\n+    cdef void* up\n+    cdef void* sp\n+    cdef int qs[2]\n+    cdef int us[2]\n+    cdef int ss[2]\n+    cdef cnp.npy_intp size\n+    cdef bint qisF\n+\n+    if q.ndim != 2: raise ValueError('q must be 2d')\n+    m = q.shape[0]\n+    n = q.shape[1]\n+    if cnp.PyArray_CHKFLAGS(q, cnp.NPY_F_CONTIGUOUS):\n+        qisF = True\n+    elif cnp.PyArray_CHKFLAGS(q, cnp.NPY_C_CONTIGUOUS):\n+        qisF = False\n+    else:\n+        raise ValueError('q must be one segment.')\n+    if u.ndim != 1: raise ValueError('u must be 1d')\n+    if u.shape[0] != m: raise ValueError('u.shape[0] must be q.shape[0]')\n+    typecode = cnp.PyArray_TYPE(q)\n+    if cnp.PyArray_TYPE(u) != typecode:\n+        raise ValueError('q and u must have the same type.')\n+\n+    if not (typecode == cnp.NPY_FLOAT or typecode == cnp.NPY_DOUBLE \\\n+            or typecode == cnp.NPY_CFLOAT or typecode == cnp.NPY_CDOUBLE):\n+        raise ValueError('q and u must be a blas compatible type: f d F or D')\n+\n+    q = validate_array(q, True)\n+    u = validate_array(u, True)\n+\n+    qp = extract(q, qs)\n+    up = extract(u, us)\n+\n+    size = 2*n\n+    s = cnp.PyArray_ZEROS(1, &size, typecode, 1)\n+    sp = extract(s, ss)\n+\n+    if typecode == cnp.NPY_FLOAT:\n+        rc_f = rcond\n+        info = reorth(m, n, <float*>qp, qs, qisF, <float*>up, us, <float*>sp, &rc_f)\n+        return u, s[:n+1], rc_f, info\n+    if typecode == cnp.NPY_DOUBLE:\n+        rc_d = rcond\n+        info = reorth(m, n, <double*>qp, qs, qisF, <double*>up, us, <double*>sp, &rc_d)\n+        return u, s[:n+1], rc_d, info\n+    if typecode == cnp.NPY_CFLOAT:\n+        rc_f = rcond\n+        info = reorth(m, n, <float_complex*>qp, qs, qisF, <float_complex*>up, us, <float_complex*>sp, <float_complex*>&rc_f)\n+        return u, s[:n+1], rc_f, info\n+    if typecode == cnp.NPY_CDOUBLE:\n+        rc_d = rcond\n+        info = reorth(m, n, <double_complex*>qp, qs, qisF, <double_complex*>up, us, <double_complex*>sp, <double_complex*>&rc_d)\n+        return u, s[:n+1], rc_d, info\n+\n+cdef int reorth(int m, int n, blas_t* q, int* qs, bint qisF, blas_t* u,\n+                int* us, blas_t* s, blas_t* RCOND) nogil:\n+    \"\"\"Given a (m,n) matrix q with orthonormal columns and a (m,) vector u,\n+       find vectors s, w and scalar p such that u = Qs + pw where w is of unit\n+       length and orthogonal to the columns of q.\n+\n+       FIX comment on return values, and RCOND \n+\n+       The method used for orthogonalizing u against q is described in [5]\n+       listed in the file header.\n+    \"\"\"\n+    cdef blas_t unorm, snorm, wnorm, wpnorm, sigma_max, sigma_min, rc\n+    cdef char* T = 'T'\n+    cdef char* N = 'N'\n+    cdef char* C = 'C'\n+    cdef int ss = 1\n+    cdef blas_t inv_root2 = <blas_t>sqrt(2)\n+\n+    # normalize u\n+    unorm = nrm2(m, u, us[0])\n+    scal(m, 1/unorm, u, us[0])\n+\n+    # decompose u into q's columns.\n+    if qisF:\n+        if blas_t is float or blas_t is double:\n+            gemv(T, m, n, 1, q, qs[1], u, us[0], 0, s, 1)\n+        else:\n+            gemv(C, m, n, 1, q, qs[1], u, us[0], 0, s, 1)\n+    else:\n+        if blas_t is float or blas_t is double:\n+            gemv(N, n, m, 1, q, n, u, us[0], 0, s, 1)\n+        else:\n+            blas_t_conj(m, u, us)\n+            gemv(N, n, m, 1, q, n, u, us[0], 0, s, 1)\n+            blas_t_conj(m, u, us)\n+            blas_t_conj(n, s, &ss)\n+\n+    # sigma_max is the largest singular value of q augmented with u/unorm\n+    snorm = nrm2(n, s, 1)\n+    sigma_max = blas_t_sqrt(1 + snorm)\n+\n+    # make u be the part of u that is not in span(q)\n+    # i.e. u -= q.dot(s)\n+    if qisF:\n+        gemv(N, m, n, -1, q, qs[1], s, 1, 1, u, us[0])\n+    else:\n+        gemv(T, n, m, -1, q, n, s, 1, 1, u, us[0])\n+    wnorm = nrm2(m, u, us[0])\n+\n+    # sigma_min is the smallest singular value of q qugmented with u/unorm\n+    # the others are == 1, since q is orthonormal.\n+    sigma_min = wnorm / sigma_max\n+    rc = sigma_min / sigma_max\n+\n+    # check the conditioning of the problem.\n+    if blas_t_less_than(rc, RCOND[0]):\n+        RCOND[0] = rc\n+        return 2\n+    RCOND[0] = rc\n+\n+    if blas_t_less_than(inv_root2, wnorm):\n+        scal(m, 1/wnorm, u, us[0])\n+        scal(n, unorm, s, 1)\n+        s[n] = unorm*wnorm\n+        return 0\n+\n+    # if the above check failed, try one reorthogonalization\n+    if qisF:\n+        if blas_t is float or blas_t is double:\n+            gemv(T, m, n, 1, q, qs[1], u, us[0], 0, s+n, 1)\n+        else:\n+            gemv(C, m, n, 1, q, qs[1], u, us[0], 0, s+n, 1)\n+        gemv(N, m, n, -1, q, qs[1], s+n, 1, 1, u, us[0])\n+    else:\n+        if blas_t is float or blas_t is double:\n+            gemv(N, n, m, 1, q, n, u, us[0], 0, s+n, 1)\n+        else:\n+            blas_t_conj(m, u, us)\n+            gemv(N, n, m, 1, q, n, u, us[0], 0, s+n, 1)\n+            blas_t_conj(m, u, us)\n+            blas_t_conj(n, s+n, &ss)\n+        gemv(T, n, m, -1, q, n, s+n, 1, 1, u, us[0])\n+        \n+    wpnorm = nrm2(m, u, us[0]) \n+    \n+    if blas_t_less_than(wpnorm, wnorm/inv_root2): # u lies in span(q) \n+        scal(m, 0, u, us[0])\n+        axpy(n, 1, s, 1, s+n, 1)\n+        scal(n, unorm, s, 1)\n+        s[n] = 0\n+        return 1\n+    scal(m, 1/wpnorm, u, us[0])\n+    axpy(n, 1, s, 1, s+n, 1)\n+    scal(n, unorm, s, 1)\n+    s[n] = wpnorm*unorm\n+    return 0\n+\n+def _form_qTu(object a, object b):\n+    \"\"\" this function only exists to expose the cdef version below for testing\n+        purposes. Here we perform minimal input validation to ensure that the\n+        inputs meet the requirements below.\n+    \"\"\"\n+    cdef cnp.ndarray q, u, qTu\n+    cdef int typecode\n+    cdef void* qTuvoid\n+    cdef int qTus[2]\n+\n+    if not cnp.PyArray_Check(a) or not cnp.PyArray_Check(b):\n+        raise ValueError('Inputs must be arrays')\n+\n+    q = a\n+    u = b\n+    \n+    typecode = cnp.PyArray_TYPE(q)\n+    if cnp.PyArray_TYPE(u) != typecode:\n+        raise ValueError('q and u must have the same type.')\n+\n+    if not (typecode == cnp.NPY_FLOAT or typecode == cnp.NPY_DOUBLE \\\n+            or typecode == cnp.NPY_CFLOAT or typecode == cnp.NPY_CDOUBLE):\n+        raise ValueError('q and u must be a blas compatible type: f d F or D')\n+\n+    q = validate_array(q, True)\n+    u = validate_array(u, True)\n+\n+    qTu = cnp.PyArray_ZEROS(u.ndim, u.shape, typecode, 1)\n+    qTuvoid = extract(qTu, qTus)\n+    form_qTu(q, u, qTuvoid, qTus, 0)\n+    return qTu\n+\n+cdef form_qTu(cnp.ndarray q, cnp.ndarray u, void* qTuvoid, int* qTus,\n+              int k):\n+    \"\"\" assuming here that q and u have compatible shapes, and are the same type\n+        + Q is contiguous.  This function is preferable over simply calling\n+        np.dot for two reasons: 1) this output is always in F order, 2) no\n+        copies need be made if Q is complex.  Point 2 in particular makes this\n+        a good bit faster for complex inputs.\n+    \"\"\"\n+    cdef int m = q.shape[0]\n+    cdef int n = q.shape[1]\n+    cdef int typecode = cnp.PyArray_TYPE(q)\n+    cdef cnp.ndarray qTu\n+    cdef char* T = 'T'\n+    cdef char* C = 'C'\n+    cdef char* N = 'N'\n+    cdef void* qvoid\n+    cdef void* uvoid\n+    cdef int qs[2]\n+    cdef int us[2]\n+    cdef int ldu\n+\n+    if cnp.PyArray_CHKFLAGS(q, cnp.NPY_F_CONTIGUOUS):\n+        qvoid = extract(q, qs)\n+        if u.ndim == 1:\n+            uvoid = extract(u, us)\n+            if typecode == cnp.NPY_FLOAT:\n+                gemv(T, m, n, 1, <float*>qvoid, qs[1],\n+                        <float*>uvoid, us[0], 0, col(<float*>qTuvoid, qTus, k), qTus[0])\n+            if typecode == cnp.NPY_DOUBLE:\n+                gemv(T, m, n, 1, <double*>qvoid, qs[1],\n+                        <double*>uvoid, us[0], 0, col(<double*>qTuvoid, qTus, k), qTus[0])\n+            if typecode == cnp.NPY_CFLOAT:\n+                gemv(C, m, n, 1, <float_complex*>qvoid, qs[1],\n+                        <float_complex*>uvoid, us[0], 0,\n+                        col(<float_complex*>qTuvoid, qTus, k), qTus[0])\n+            if typecode == cnp.NPY_CDOUBLE:\n+                gemv(C, m, n, 1, <double_complex*>qvoid, qs[1],\n+                        <double_complex*>uvoid, us[0], 0,\n+                        col(<double_complex*>qTuvoid, qTus, k), qTus[0])\n+        elif u.ndim == 2:\n+            p = u.shape[1]\n+            if cnp.PyArray_CHKFLAGS(u, cnp.NPY_F_CONTIGUOUS):\n+                utrans = N\n+                uvoid = extract(u, us)\n+                ldu = us[1]\n+            elif cnp.PyArray_CHKFLAGS(u, cnp.NPY_C_CONTIGUOUS):\n+                utrans = T\n+                uvoid = extract(u, us)\n+                ldu = us[0]\n+            else:\n+                u = PyArray_FromArraySafe(u, NULL, cnp.NPY_F_CONTIGUOUS)\n+                utrans = N\n+                uvoid = extract(u, us)\n+                ldu = us[1]\n+            if typecode == cnp.NPY_FLOAT:\n+                gemm(T, utrans, m, p, m, 1, <float*>qvoid, qs[1],\n+                        <float*>uvoid, ldu, 0, col(<float*>qTuvoid, qTus, k), qTus[1])\n+            if typecode == cnp.NPY_DOUBLE:\n+                gemm(T, utrans, m, p, m, 1, <double*>qvoid, qs[1],\n+                        <double*>uvoid, ldu, 0, col(<double*>qTuvoid, qTus, k), qTus[1])\n+            if typecode == cnp.NPY_CFLOAT:\n+                gemm(C, utrans, m, p, m, 1, <float_complex*>qvoid, qs[1],\n+                        <float_complex*>uvoid, ldu, 0,\n+                        col(<float_complex*>qTuvoid, qTus, k), qTus[1])\n+            if typecode == cnp.NPY_CDOUBLE:\n+                gemm(C, utrans, m, p, m, 1, <double_complex*>qvoid, qs[1],\n+                        <double_complex*>uvoid, ldu, 0,\n+                        col(<double_complex*>qTuvoid, qTus, k), qTus[1])\n+    \n+    elif cnp.PyArray_CHKFLAGS(q, cnp.NPY_C_CONTIGUOUS):\n+        qvoid = extract(q, qs)\n+        if u.ndim == 1:\n+            uvoid = extract(u, us)\n+            if typecode == cnp.NPY_FLOAT:\n+                gemv(N, m, n, 1, <float*>qvoid, qs[0],\n+                        <float*>uvoid, us[0], 0, col(<float*>qTuvoid, qTus, k), qTus[0])\n+            if typecode == cnp.NPY_DOUBLE:\n+                gemv(N, m, n, 1, <double*>qvoid, qs[0],\n+                        <double*>uvoid, us[0], 0, col(<double*>qTuvoid, qTus, k), qTus[0])\n+            if typecode == cnp.NPY_CFLOAT:\n+                blas_t_conj(m, <float_complex*>uvoid, us)\n+                gemv(N, m, n, 1, <float_complex*>qvoid, qs[0],\n+                        <float_complex*>uvoid, us[0], 0,\n+                        col(<float_complex*>qTuvoid, qTus, k), qTus[0])\n+                blas_t_conj(m, col(<float_complex*>qTuvoid, qTus, k), qTus)\n+            if typecode == cnp.NPY_CDOUBLE:\n+                blas_t_conj(m, <double_complex*>uvoid, us)\n+                gemv(N, m, n, 1, <double_complex*>qvoid, qs[0],\n+                        <double_complex*>uvoid, us[0], 0,\n+                        col(<double_complex*>qTuvoid, qTus, k), qTus[0])\n+                blas_t_conj(m, col(<double_complex*>qTuvoid, qTus, k), qTus)\n+        elif u.ndim == 2:\n+            p = u.shape[1]\n+            if cnp.PyArray_CHKFLAGS(u, cnp.NPY_F_CONTIGUOUS):\n+                utrans = N\n+                uvoid = extract(u, us)\n+                ldu = us[1]\n+            elif cnp.PyArray_CHKFLAGS(u, cnp.NPY_C_CONTIGUOUS):\n+                utrans = T\n+                uvoid = extract(u, us)\n+                ldu = us[0]\n+            else:\n+                u = PyArray_FromArraySafe(u, NULL, cnp.NPY_F_CONTIGUOUS)\n+                utrans = N\n+                uvoid = extract(u, us)\n+                ldu = us[1]\n+            if typecode == cnp.NPY_FLOAT:\n+                gemm(N, utrans, m, p, m, 1, <float*>qvoid, qs[0],\n+                        <float*>uvoid, ldu, 0, col(<float*>qTuvoid, qTus, k), qTus[1])\n+            elif typecode == cnp.NPY_DOUBLE:\n+                gemm(N, utrans, m, p, m, 1, <double*>qvoid, qs[0],\n+                        <double*>uvoid, ldu, 0, col(<double*>qTuvoid, qTus, k), qTus[1])\n+            elif typecode == cnp.NPY_CFLOAT:\n+                blas_t_2d_conj(m, p, <float_complex*>uvoid, us)\n+                gemm(N, utrans, m, p, m, 1, <float_complex*>qvoid, qs[0],\n+                        <float_complex*>uvoid, ldu, 0,\n+                        col(<float_complex*>qTuvoid, qTus, k), qTus[1])\n+                blas_t_2d_conj(m, p, col(<float_complex*>qTuvoid, qTus, k), qTus)\n+            elif typecode == cnp.NPY_CDOUBLE:\n+                blas_t_2d_conj(m, p, <double_complex*>uvoid, us)\n+                gemm(N, utrans, m, p, m, 1, <double_complex*>qvoid, qs[0],\n+                        <double_complex*>uvoid, ldu, 0,\n+                        col(<double_complex*>qTuvoid, qTus, k), qTus[1])\n+                blas_t_2d_conj(m, p, col(<double_complex*>qTuvoid, qTus, k), qTus)\n+        else:\n+            raise ValueError('1 <= u.ndim <= 2')\n+    else:\n+        raise ValueError('q must be either F or C contig')\n+\n+cdef validate_array(cnp.ndarray a, bint chkfinite):\n+    # here we check that a has positive strides and that its size is small\n+    # enough to fit in into an int, as BLAS/LAPACK require\n+    cdef bint copy = False\n+    cdef bint too_large = False\n+    cdef int j\n+\n+    for j in range(a.ndim):\n+        if a.strides[j] <= 0 or \\\n+                (a.strides[j] / a.descr.itemsize) >= limits.INT_MAX:\n+            copy = True\n+        if a.shape[j] >= limits.INT_MAX:\n+            raise ValueError('Input array to large for use with BLAS')\n+\n+    if chkfinite:\n+        if not np.isfinite(a).all():\n+            raise ValueError('array must not contain infs or NaNs')\n+\n+    if copy:\n+            return PyArray_FromArraySafe(a, NULL, cnp.NPY_F_CONTIGUOUS)\n+    return a\n+\n+cdef validate_qr(object q0, object r0, bint overwrite_q, int q_order,\n+                 bint overwrite_r, int r_order, bint chkfinite):\n+    cdef cnp.ndarray Q\n+    cdef cnp.ndarray R\n+    cdef int typecode\n+    cdef bint economic = False\n+\n+    q_order |= cnp.NPY_BEHAVED_NS | cnp.NPY_ELEMENTSTRIDES  \n+    r_order |= cnp.NPY_BEHAVED_NS | cnp.NPY_ELEMENTSTRIDES\n+\n+    if not overwrite_q:\n+        q_order |= cnp.NPY_ENSURECOPY\n+\n+    if not overwrite_r:\n+        r_order |= cnp.NPY_ENSURECOPY\n+\n+    # in the interests of giving better error messages take any number of\n+    # dimensions here.\n+    Q = PyArray_CheckFromAny(q0, NULL, 0, 0, q_order, NULL)\n+    R = PyArray_CheckFromAny(r0, NULL, 0, 0, r_order, NULL)\n+\n+    if Q.ndim != 2 or R.ndim != 2:\n+        raise ValueError('Q and R must be 2d')","path":"scipy/linalg/_decomp_update.pyx","position":1485,"original_position":1485,"commit_id":"ff7d83a579753edf98a1b0d440ff314650049bb0","original_commit_id":"ff7d83a579753edf98a1b0d440ff314650049bb0","user":{"login":"rgommers","id":98330,"avatar_url":"https://avatars.githubusercontent.com/u/98330?v=3","gravatar_id":"","url":"https://api.github.com/users/rgommers","html_url":"https://github.com/rgommers","followers_url":"https://api.github.com/users/rgommers/followers","following_url":"https://api.github.com/users/rgommers/following{/other_user}","gists_url":"https://api.github.com/users/rgommers/gists{/gist_id}","starred_url":"https://api.github.com/users/rgommers/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/rgommers/subscriptions","organizations_url":"https://api.github.com/users/rgommers/orgs","repos_url":"https://api.github.com/users/rgommers/repos","events_url":"https://api.github.com/users/rgommers/events{/privacy}","received_events_url":"https://api.github.com/users/rgommers/received_events","type":"User","site_admin":false},"body":"nitpick: ``2d`` -->  ``2-D``","created_at":"2015-03-15T10:18:25Z","updated_at":"2015-03-15T10:18:25Z","html_url":"https://github.com/scipy/scipy/pull/4249#discussion_r26446344","pull_request_url":"https://api.github.com/repos/scipy/scipy/pulls/4249","_links":{"self":{"href":"https://api.github.com/repos/scipy/scipy/pulls/comments/26446344"},"html":{"href":"https://github.com/scipy/scipy/pull/4249#discussion_r26446344"},"pull_request":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4249"}}},"pull_request":{"url":"https://api.github.com/repos/scipy/scipy/pulls/4249","id":25745696,"html_url":"https://github.com/scipy/scipy/pull/4249","diff_url":"https://github.com/scipy/scipy/pull/4249.diff","patch_url":"https://github.com/scipy/scipy/pull/4249.patch","issue_url":"https://api.github.com/repos/scipy/scipy/issues/4249","number":4249,"state":"open","locked":false,"title":"Add routines for updating QR decompositions","user":{"login":"ewmoore","id":2314641,"avatar_url":"https://avatars.githubusercontent.com/u/2314641?v=3","gravatar_id":"","url":"https://api.github.com/users/ewmoore","html_url":"https://github.com/ewmoore","followers_url":"https://api.github.com/users/ewmoore/followers","following_url":"https://api.github.com/users/ewmoore/following{/other_user}","gists_url":"https://api.github.com/users/ewmoore/gists{/gist_id}","starred_url":"https://api.github.com/users/ewmoore/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/ewmoore/subscriptions","organizations_url":"https://api.github.com/users/ewmoore/orgs","repos_url":"https://api.github.com/users/ewmoore/repos","events_url":"https://api.github.com/users/ewmoore/events{/privacy}","received_events_url":"https://api.github.com/users/ewmoore/received_events","type":"User","site_admin":false},"body":"This pull request goes on top of #4021, accordingly only the last few commits are relevant here.\r\n\r\nHere three new functions are added to perform various modifications of qr decompositions: `qr_update`, `qr_insert` and `qr_delete` for performing rank-k updates, adding row(s) or column(s) and removing row(s) or column(s).  \r\n\r\nThe methods used are essentially those in [1], with some input from [2].\r\n\r\nI've endeavored to write this such that updates can be performed in-place as much as possible.  This is, for instance, why the rank-1 update and the rank-k update algorithms are slightly different.  This allows the rank-1 update to support updating `q` and `r` when they are stored in C order.  (In this spirit, it might be worth changing `linalg.qr` to return fortran ordered 'r', since at present it always returns a C ordered `r`,  and a F ordered `q`.  Numpy's qr returns both as C order.)\r\n\r\nSome outstanding things:\r\n\r\n1.  I suppose I'll need to add a 'check_finite'  flag to these.  \r\n2.  `overwrite_qr` should probably default to `False` instead of true.\r\n\r\nReferences\r\n\r\n1.  Golub, G. H. & Loan, C. F. van V. Matrix Computations, 3rd Ed. (Johns Hopkins University Press, 1996).\r\n2.  Hammarling, S. & Lucas, C. Updating the QR factorization and the least squares problem. 1-73 (The University of Manchester, 2008). at <http://eprints.ma.man.ac.uk/1192/>\r\n","created_at":"2014-12-09T14:18:50Z","updated_at":"2015-03-15T10:18:25Z","closed_at":null,"merged_at":null,"merge_commit_sha":"17a841665ec903797737bc471a2d90921363c200","assignee":null,"milestone":null,"commits_url":"https://api.github.com/repos/scipy/scipy/pulls/4249/commits","review_comments_url":"https://api.github.com/repos/scipy/scipy/pulls/4249/comments","review_comment_url":"https://api.github.com/repos/scipy/scipy/pulls/comments{/number}","comments_url":"https://api.github.com/repos/scipy/scipy/issues/4249/comments","statuses_url":"https://api.github.com/repos/scipy/scipy/statuses/ff7d83a579753edf98a1b0d440ff314650049bb0","head":{"label":"ewmoore:updateqr2","ref":"updateqr2","sha":"ff7d83a579753edf98a1b0d440ff314650049bb0","user":{"login":"ewmoore","id":2314641,"avatar_url":"https://avatars.githubusercontent.com/u/2314641?v=3","gravatar_id":"","url":"https://api.github.com/users/ewmoore","html_url":"https://github.com/ewmoore","followers_url":"https://api.github.com/users/ewmoore/followers","following_url":"https://api.github.com/users/ewmoore/following{/other_user}","gists_url":"https://api.github.com/users/ewmoore/gists{/gist_id}","starred_url":"https://api.github.com/users/ewmoore/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/ewmoore/subscriptions","organizations_url":"https://api.github.com/users/ewmoore/orgs","repos_url":"https://api.github.com/users/ewmoore/repos","events_url":"https://api.github.com/users/ewmoore/events{/privacy}","received_events_url":"https://api.github.com/users/ewmoore/received_events","type":"User","site_admin":false},"repo":{"id":5864721,"name":"scipy","full_name":"ewmoore/scipy","owner":{"login":"ewmoore","id":2314641,"avatar_url":"https://avatars.githubusercontent.com/u/2314641?v=3","gravatar_id":"","url":"https://api.github.com/users/ewmoore","html_url":"https://github.com/ewmoore","followers_url":"https://api.github.com/users/ewmoore/followers","following_url":"https://api.github.com/users/ewmoore/following{/other_user}","gists_url":"https://api.github.com/users/ewmoore/gists{/gist_id}","starred_url":"https://api.github.com/users/ewmoore/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/ewmoore/subscriptions","organizations_url":"https://api.github.com/users/ewmoore/orgs","repos_url":"https://api.github.com/users/ewmoore/repos","events_url":"https://api.github.com/users/ewmoore/events{/privacy}","received_events_url":"https://api.github.com/users/ewmoore/received_events","type":"User","site_admin":false},"private":false,"html_url":"https://github.com/ewmoore/scipy","description":"Scipy main repository","fork":true,"url":"https://api.github.com/repos/ewmoore/scipy","forks_url":"https://api.github.com/repos/ewmoore/scipy/forks","keys_url":"https://api.github.com/repos/ewmoore/scipy/keys{/key_id}","collaborators_url":"https://api.github.com/repos/ewmoore/scipy/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/ewmoore/scipy/teams","hooks_url":"https://api.github.com/repos/ewmoore/scipy/hooks","issue_events_url":"https://api.github.com/repos/ewmoore/scipy/issues/events{/number}","events_url":"https://api.github.com/repos/ewmoore/scipy/events","assignees_url":"https://api.github.com/repos/ewmoore/scipy/assignees{/user}","branches_url":"https://api.github.com/repos/ewmoore/scipy/branches{/branch}","tags_url":"https://api.github.com/repos/ewmoore/scipy/tags","blobs_url":"https://api.github.com/repos/ewmoore/scipy/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/ewmoore/scipy/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/ewmoore/scipy/git/refs{/sha}","trees_url":"https://api.github.com/repos/ewmoore/scipy/git/trees{/sha}","statuses_url":"https://api.github.com/repos/ewmoore/scipy/statuses/{sha}","languages_url":"https://api.github.com/repos/ewmoore/scipy/languages","stargazers_url":"https://api.github.com/repos/ewmoore/scipy/stargazers","contributors_url":"https://api.github.com/repos/ewmoore/scipy/contributors","subscribers_url":"https://api.github.com/repos/ewmoore/scipy/subscribers","subscription_url":"https://api.github.com/repos/ewmoore/scipy/subscription","commits_url":"https://api.github.com/repos/ewmoore/scipy/commits{/sha}","git_commits_url":"https://api.github.com/repos/ewmoore/scipy/git/commits{/sha}","comments_url":"https://api.github.com/repos/ewmoore/scipy/comments{/number}","issue_comment_url":"https://api.github.com/repos/ewmoore/scipy/issues/comments{/number}","contents_url":"https://api.github.com/repos/ewmoore/scipy/contents/{+path}","compare_url":"https://api.github.com/repos/ewmoore/scipy/compare/{base}...{head}","merges_url":"https://api.github.com/repos/ewmoore/scipy/merges","archive_url":"https://api.github.com/repos/ewmoore/scipy/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/ewmoore/scipy/downloads","issues_url":"https://api.github.com/repos/ewmoore/scipy/issues{/number}","pulls_url":"https://api.github.com/repos/ewmoore/scipy/pulls{/number}","milestones_url":"https://api.github.com/repos/ewmoore/scipy/milestones{/number}","notifications_url":"https://api.github.com/repos/ewmoore/scipy/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/ewmoore/scipy/labels{/name}","releases_url":"https://api.github.com/repos/ewmoore/scipy/releases{/id}","created_at":"2012-09-19T00:02:55Z","updated_at":"2014-06-18T15:15:18Z","pushed_at":"2015-03-12T15:32:32Z","git_url":"git://github.com/ewmoore/scipy.git","ssh_url":"git@github.com:ewmoore/scipy.git","clone_url":"https://github.com/ewmoore/scipy.git","svn_url":"https://github.com/ewmoore/scipy","homepage":"http://scipy.org/","size":53867,"stargazers_count":0,"watchers_count":0,"language":"C","has_issues":false,"has_downloads":true,"has_wiki":false,"has_pages":false,"forks_count":0,"mirror_url":null,"open_issues_count":0,"forks":0,"open_issues":0,"watchers":0,"default_branch":"master"}},"base":{"label":"scipy:master","ref":"master","sha":"a17d7144511162516e000fa487992edd7379ab22","user":{"login":"scipy","id":288277,"avatar_url":"https://avatars.githubusercontent.com/u/288277?v=3","gravatar_id":"","url":"https://api.github.com/users/scipy","html_url":"https://github.com/scipy","followers_url":"https://api.github.com/users/scipy/followers","following_url":"https://api.github.com/users/scipy/following{/other_user}","gists_url":"https://api.github.com/users/scipy/gists{/gist_id}","starred_url":"https://api.github.com/users/scipy/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/scipy/subscriptions","organizations_url":"https://api.github.com/users/scipy/orgs","repos_url":"https://api.github.com/users/scipy/repos","events_url":"https://api.github.com/users/scipy/events{/privacy}","received_events_url":"https://api.github.com/users/scipy/received_events","type":"Organization","site_admin":false},"repo":{"id":1460385,"name":"scipy","full_name":"scipy/scipy","owner":{"login":"scipy","id":288277,"avatar_url":"https://avatars.githubusercontent.com/u/288277?v=3","gravatar_id":"","url":"https://api.github.com/users/scipy","html_url":"https://github.com/scipy","followers_url":"https://api.github.com/users/scipy/followers","following_url":"https://api.github.com/users/scipy/following{/other_user}","gists_url":"https://api.github.com/users/scipy/gists{/gist_id}","starred_url":"https://api.github.com/users/scipy/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/scipy/subscriptions","organizations_url":"https://api.github.com/users/scipy/orgs","repos_url":"https://api.github.com/users/scipy/repos","events_url":"https://api.github.com/users/scipy/events{/privacy}","received_events_url":"https://api.github.com/users/scipy/received_events","type":"Organization","site_admin":false},"private":false,"html_url":"https://github.com/scipy/scipy","description":"Scipy library main repository","fork":false,"url":"https://api.github.com/repos/scipy/scipy","forks_url":"https://api.github.com/repos/scipy/scipy/forks","keys_url":"https://api.github.com/repos/scipy/scipy/keys{/key_id}","collaborators_url":"https://api.github.com/repos/scipy/scipy/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/scipy/scipy/teams","hooks_url":"https://api.github.com/repos/scipy/scipy/hooks","issue_events_url":"https://api.github.com/repos/scipy/scipy/issues/events{/number}","events_url":"https://api.github.com/repos/scipy/scipy/events","assignees_url":"https://api.github.com/repos/scipy/scipy/assignees{/user}","branches_url":"https://api.github.com/repos/scipy/scipy/branches{/branch}","tags_url":"https://api.github.com/repos/scipy/scipy/tags","blobs_url":"https://api.github.com/repos/scipy/scipy/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/scipy/scipy/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/scipy/scipy/git/refs{/sha}","trees_url":"https://api.github.com/repos/scipy/scipy/git/trees{/sha}","statuses_url":"https://api.github.com/repos/scipy/scipy/statuses/{sha}","languages_url":"https://api.github.com/repos/scipy/scipy/languages","stargazers_url":"https://api.github.com/repos/scipy/scipy/stargazers","contributors_url":"https://api.github.com/repos/scipy/scipy/contributors","subscribers_url":"https://api.github.com/repos/scipy/scipy/subscribers","subscription_url":"https://api.github.com/repos/scipy/scipy/subscription","commits_url":"https://api.github.com/repos/scipy/scipy/commits{/sha}","git_commits_url":"https://api.github.com/repos/scipy/scipy/git/commits{/sha}","comments_url":"https://api.github.com/repos/scipy/scipy/comments{/number}","issue_comment_url":"https://api.github.com/repos/scipy/scipy/issues/comments{/number}","contents_url":"https://api.github.com/repos/scipy/scipy/contents/{+path}","compare_url":"https://api.github.com/repos/scipy/scipy/compare/{base}...{head}","merges_url":"https://api.github.com/repos/scipy/scipy/merges","archive_url":"https://api.github.com/repos/scipy/scipy/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/scipy/scipy/downloads","issues_url":"https://api.github.com/repos/scipy/scipy/issues{/number}","pulls_url":"https://api.github.com/repos/scipy/scipy/pulls{/number}","milestones_url":"https://api.github.com/repos/scipy/scipy/milestones{/number}","notifications_url":"https://api.github.com/repos/scipy/scipy/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/scipy/scipy/labels{/name}","releases_url":"https://api.github.com/repos/scipy/scipy/releases{/id}","created_at":"2011-03-09T18:52:03Z","updated_at":"2015-03-14T23:07:50Z","pushed_at":"2015-03-14T23:07:50Z","git_url":"git://github.com/scipy/scipy.git","ssh_url":"git@github.com:scipy/scipy.git","clone_url":"https://github.com/scipy/scipy.git","svn_url":"https://github.com/scipy/scipy","homepage":"http://scipy.org/scipylib/","size":84233,"stargazers_count":1536,"watchers_count":1536,"language":"Python","has_issues":true,"has_downloads":true,"has_wiki":true,"has_pages":false,"forks_count":812,"mirror_url":null,"open_issues_count":784,"forks":812,"open_issues":784,"watchers":1536,"default_branch":"master"}},"_links":{"self":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4249"},"html":{"href":"https://github.com/scipy/scipy/pull/4249"},"issue":{"href":"https://api.github.com/repos/scipy/scipy/issues/4249"},"comments":{"href":"https://api.github.com/repos/scipy/scipy/issues/4249/comments"},"review_comments":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4249/comments"},"review_comment":{"href":"https://api.github.com/repos/scipy/scipy/pulls/comments{/number}"},"commits":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4249/commits"},"statuses":{"href":"https://api.github.com/repos/scipy/scipy/statuses/ff7d83a579753edf98a1b0d440ff314650049bb0"}}}}