--- a/Makeconf.in	Tue Jan 20 21:15:17 2009 +0100
+++ b/Makeconf.in	Tue Jan 20 21:16:42 2009 +0100
@@ -239,6 +239,7 @@
 CHOLMOD_LIBS = @CHOLMOD_LIBS@
 CXSPARSE_LIBS = @CXSPARSE_LIBS@
 OPENGL_LIBS = @OPENGL_LIBS@
+QRUPDATE_LIBS = @QRUPDATE_LIBS@
 ARPACK_LIBS = @ARPACK_LIBS@
 LIBS = @LIBS@
 

--- a/configure.in	Tue Jan 20 21:15:17 2009 +0100
+++ b/configure.in	Tue Jan 20 21:16:42 2009 +0100
@@ -865,6 +865,34 @@
   AC_MSG_WARN($warn_blas_f77_incompatible)
 fi
 
+QRUPDATE_LIBS=
+AC_SUBST(QRUPDATE_LIBS)
+
+# Check for the qrupdate library
+AC_ARG_WITH(qrupdate,
+  [AS_HELP_STRING([--without-qrupdate],
+     [don't use qrupdate, disable QR & Cholesky updating functions])],
+  with_qrupdate=$withval, with_qrupdate=yes)
+
+warn_qrupdate="qrupdate not found. The QR & Cholesky updating function will not be available."
+if test "$with_qrupdate" = yes; then
+  with_qrupdate=no
+  if $have_fortran_compiler; then 
+    AC_F77_FUNC(sqr1up)
+  elif $have_f2c; then
+    sqr1up=sqr1up_
+  fi
+  AC_CHECK_LIB(qrupdate, $sqr1up, 
+    [QRUPDATE_LIBS="-lqrupdate"; with_qrupdate=yes], [], [$BLAS_LIBS $FLIBS])
+  if test "$with_qrupdate" = yes; then
+    AC_DEFINE(HAVE_QRUPDATE, 1, [Define if the qrupdate library is used.])
+    warn_qrupdate=
+  fi
+fi
+if test -n "$warn_qrupdate"; then
+  AC_MSG_WARN($warn_qrupdate)
+fi
+
 # Check for AMD library
 AMD_LIBS=
 AC_SUBST(AMD_LIBS)
@@ -1997,7 +2025,6 @@
   libcruft/lapack/Makefile 
   libcruft/misc/Makefile libcruft/odepack/Makefile
   libcruft/ordered-qz/Makefile libcruft/quadpack/Makefile
-  libcruft/qrupdate/Makefile
   libcruft/ranlib/Makefile libcruft/slatec-fn/Makefile
   libcruft/slatec-err/Makefile libcruft/villad/Makefile
   libcruft/blas-xtra/Makefile libcruft/lapack-xtra/Makefile])
@@ -2034,6 +2061,7 @@
   CHOLMOD libraries:    $CHOLMOD_LIBS
   CXSPARSE libraries:   $CXSPARSE_LIBS
   ARPACK libraries:     $ARPACK_LIBS
+  QRUPDATE libraries:	$QRUPDATE_LIBS
   HDF5 libraries:       $HDF5_LIBS
   CURL libraries:       $CURL_LIBS
   REGEX libraries:      $REGEX_LIBS
@@ -2120,6 +2148,11 @@
   warn_msg_printed=true
 fi
 
+if test -n "$warn_qrupdate"; then
+  AC_MSG_WARN($warn_qrupdate)
+  warn_msg_printed=true
+fi
+
 if test -n "$warn_amd"; then
   AC_MSG_WARN($warn_amd)
   warn_msg_printed=true

--- a/libcruft/Makefile.in	Tue Jan 20 21:15:17 2009 +0100
+++ b/libcruft/Makefile.in	Tue Jan 20 21:16:42 2009 +0100
@@ -43,7 +43,7 @@
 
 CRUFT_DIRS = amos @BLAS_DIR@ blas-xtra daspk dasrt dassl \
 	@FFT_DIR@ @LAPACK_DIR@ lapack-xtra \
-	misc odepack ordered-qz quadpack qrupdate ranlib \
+	misc odepack ordered-qz quadpack ranlib \
 	slatec-err slatec-fn villad
 
 SUBDIRS = $(CRUFT_DIRS)

--- a/libcruft/qrupdate/Makefile.in	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,58 +0,0 @@
-# Makefile for octave's libcruft/factupd directory
-#
-# Copyright (C) 2008 VZLU, a.s.
-#
-# This file is part of Octave.
-# 
-# Octave is free software; you can redistribute it and/or modify it
-# under the terms of the GNU General Public License as published by the
-# Free Software Foundation; either version 3 of the License, or (at
-# your option) any later version.
-# 
-# Octave is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-# for more details.
-# 
-# You should have received a copy of the GNU General Public License
-# along with Octave; see the file COPYING.  If not, see
-# <http://www.gnu.org/licenses/>.
-
-TOPDIR = ../..
-
-srcdir = @srcdir@
-top_srcdir = @top_srcdir@
-VPATH = @srcdir@
-
-EXTERNAL_DISTFILES = $(DISTFILES)
-
-FSRC = dqr1up.f zqr1up.f \
-       dqrinc.f zqrinc.f \
-       dqrdec.f zqrdec.f \
-       dqrinr.f zqrinr.f \
-       dqrder.f zqrder.f \
-       dqrshc.f zqrshc.f \
-       dch1up.f zch1up.f \
-       dch1dn.f zch1dn.f \
-       dchinx.f zchinx.f \
-       dchdex.f zchdex.f \
-       dqrqhu.f zqrqhu.f \
-       dqrqhv.f zqrqhv.f \
-       dqhqr.f zqhqr.f \
-       sch1up.f cch1up.f \
-       sqrinc.f cqrinc.f \
-       sqrdec.f cqrdec.f \
-       sqrinr.f cqrinr.f \
-       sqrder.f cqrder.f \
-       sqrshc.f cqrshc.f \
-       sqr1up.f cqr1up.f \
-       sch1dn.f cch1dn.f \
-       schinx.f cchinx.f \
-       schdex.f cchdex.f \
-       sqrqhu.f cqrqhu.f \
-       sqrqhv.f cqrqhv.f \
-       sqhqr.f cqhqr.f
-
-include $(TOPDIR)/Makeconf
-
-include ../Makerules

--- a/libcruft/qrupdate/cch1dn.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cch1dn(n,R,u,w,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a hermitian positive definite matrix A, i.e.
-c               A = R'*R, this subroutine downdates R -> R1 so that
-c               R1'*R1 = A - u*u' 
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n,info
-      complex R(n,n),u(n)
-      real w(n)
-      external ctrsv,clartg,scnrm2
-      real rho,scnrm2
-      complex crho,rr,ui,t
-      integer i,j
-
-c quick return if possible
-      if (n <= 0) return
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0e0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call ctrsv('U','C','N',n,R,n,u,1)
-      rho = scnrm2(n,u,1)
-c check positive definiteness      
-      rho = 1 - rho**2
-      if (rho <= 0e0) then
-        info = 1
-        return
-      end if
-      crho = sqrt(rho)
-c eliminate R' \ u
-      do i = n,1,-1
-        ui = u(i)
-c generate next rotation        
-        call clartg(crho,ui,w(i),u(i),rr)
-        crho = rr
-      end do
-c apply rotations
-      do i = n,1,-1
-        ui = 0e0
-        do j = i,1,-1
-          t = w(j)*ui + u(j)*R(j,i)
-          R(j,i) = w(j)*R(j,i) - conjg(u(j))*ui
-          ui = t
-        end do
-      end do
-
-      info = 0
-      end

--- a/libcruft/qrupdate/cch1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,56 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cch1up(n,R,u,w)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a hermitian positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A + u*u' or A - u*u'
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a real workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n
-      complex R(n,n),u(n)
-      real w(n)
-      external clartg
-      complex rr,ui,t
-      integer i,j
-      
-      do i = 1,n
-c apply stored rotations, column-wise
-        ui = conjg(u(i))
-        do j = 1,i-1
-          t = w(j)*R(j,i) + u(j)*ui
-          ui = w(j)*ui - conjg(u(j))*R(j,i)
-          R(j,i) = t
-        end do
-c generate next rotation        
-        call clartg(R(i,i),ui,w(i),u(i),rr)
-        R(i,i) = rr
-      end do
-      end
-

--- a/libcruft/qrupdate/cchdex.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,62 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine cchdex(n,R,R1,j)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the deleted row/column
-      integer n,j,info
-      complex R(n,n),R1(n-1,n-1)
-      real c
-      complex Qdum,s,rr
-      external xerbla,clacpy,cqhqr,clartg
-
-c quick return if possible
-      if (n == 1) return
-
-c check arguments      
-      info = 0
-      if (n <= 0) then
-        info = 1
-      else if (j < 1 .or. j > n) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('CCHDEX',info)
-      end if
-
-c setup the new matrix R1
-      if (j > 1) then
-        call clacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1)
-      end if
-      if (j < n) then
-        call clacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1)
-        call cqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1)
-c eliminate R(n,n)      
-        call clartg(R1(n-1,n-1),R(n,n),c,s,rr)
-        R1(n-1,n-1) = rr
-      endif
-      end

--- a/libcruft/qrupdate/cchinx.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,109 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine cchinx(n,R,R1,j,u,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
-c               jj = [1:j-1,j+1:n+1].
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the inserted row/column
-c u (in)        the vector (n+1) determining the rank-1 update
-c info (out)    on exit, if info = 1, the 
-c               definiteness.
-
-      integer n,j,info
-      complex R(n,n),R1(n+1,n+1),u(n+1)
-      real rho,scnrm2
-      complex Qdum,w
-      external xerbla,ccopy,clacpy,ctrsv,scnrm2,cqrqhu
-      integer jj
-
-c quick return if possible
-      if (n == 0) then
-        if (real(u(1)) <= 0) then
-          info = 1
-          return
-        else
-          R(1,1) = sqrt(real(u(1)))
-        end if
-      end if
-
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 1
-      else if (j < 1 .or. j > n+1) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('CCHINX',info)
-      end if
-
-c copy shifted vector
-      if (j > 1) then
-        call ccopy(j-1,u,1,R1(1,j),1)
-      end if
-      w = u(j)
-      if (j < n+1) then
-        call ccopy(n-j+1,u(j+1),1,R1(j,j),1)
-      end if
-      
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0e0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call ctrsv('U','T','N',n,R,n,R1(1,j),1)
-      rho = scnrm2(n,R1(1,j),1)
-c check positive definiteness      
-      rho = u(j) - rho**2
-      if (rho <= 0e0) then
-        info = 1
-        return
-      end if
-      R1(n+1,n+1) = sqrt(rho)
-
-c setup the new matrix R1
-      do i = 1,n+1
-        R1(n+1,i) = 0e0
-      end do
-      if (j > 1) then
-        call clacpy('0',j-1,n,R(1,1),n,R1(1,1),n+1)
-      end if
-      if (j <= n) then
-        call clacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1)
-c retriangularize
-        jj = min(j+1,n)
-        call cqrqhu(0,n+1-j,n-j,Qdum,1,R1(j,jj),n+1,R1(j,j),w)
-        R1(j,j) = w
-        do jj = j+1,n
-          R1(jj,j) = 0e0
-        end do
-      end if
-
-      end

--- a/libcruft/qrupdate/cqhqr.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqhqr(m,n,k,Q,ldq,R,ldr)
-c purpose:      given an k-by-n upper Hessenberg matrix R and
-c               an m-by-k matrix Q, this subroutine updates
-c               R -> R1 and Q -> Q1 so that R1 is upper 
-c               trapezoidal, R1 = G*R and Q1 = Q*G', where
-c               G is an unitary matrix, giving Q1*R1 = Q*R.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q
-c n (in)        number of columns of the matrix R
-c k (in)        number of columns of Q and rows of R.
-c Q (io)        on entry, the unitary matrix Q
-c               on exit, the updated matrix Q1
-c ldq (in)      leading dimension of Q
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated upper Hessenberg matrix R1
-c ldr (in)      leading dimension of R
-c
-      integer m,n,k,ldq,ldr
-      complex Q(ldq,*),R(ldr,*)
-      real c
-      complex s,rr
-      external xerbla,clartg,crot
-      integer info,i
-c quick return if possible.
-      if (n <= 0 .or. k <= 1) return
-c check arguments.
-      info = 0
-      if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('CQHQR',info)
-      end if
-c triangularize      
-      do i = 1,min(k-1,n)
-        call clartg(R(i,i),R(i+1,i),c,s,rr)
-        R(i,i) = rr
-        R(i+1,i) = 0e0
-        if (i < n) then
-          call crot(n-i,R(i,i+1),ldr,R(i+1,i+1),ldr,c,s)
-        end if
-c apply rotation to Q        
-        if (m > 0) then
-          call crot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s))
-        end if
-      end do
-      end 
-

--- a/libcruft/qrupdate/cqr1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,53 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqr1up(m,n,k,Q,R,u,v)
-c purpose:      updates a QR factorization after rank-1 modification
-c               i.e., given a m-by-k unitary Q and m-by-n upper 
-c               trapezoidal R, an m-vector u and n-vector v, 
-c               this subroutine updates Q -> Q1 and R -> R1 so that
-c               Q1*R1 = Q*R + Q*Q'u*v', and Q1 is again unitary
-c               and R1 upper trapezoidal.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R. k <= m.
-c Q (io)        on entry, the unitary m-by-k matrix Q.
-c               on exit, the updated matrix Q1.
-c R (io)        on entry, the upper trapezoidal m-by-n matrix R.
-c               on exit, the updated matrix R1.
-c u (in)        the left m-vector.
-c v (in)        the right n-vector.
-c
-      integer m,n,k
-      complex Q(m,k),R(k,n),u(m),v(n)
-      complex w
-      external cqrqhv,cqhqr
-      integer i
-c quick return if possible      
-      if (m <= 0 .or. n <= 0) return
-c eliminate tail of Q'*u
-      call cqrqhv(m,n,k,Q,m,R,m,u,w)
-c update R      
-      do i = 1,n
-        R(1,i) = R(1,i) + w*conjg(v(i))
-      end do
-c retriangularize R
-      call cqhqr(m,n,k,Q,m,R,k)
-      end 

--- a/libcruft/qrupdate/cqrdec.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,66 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqrdec(m,n,k,Q,R,R1,j)
-c purpose:      updates a QR factorization after deleting
-c               a column.      
-c               i.e., given an m-by-k unitary matrix Q, an k-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               forms an m-by-(n-1) matrix R1 so that Q1 remains
-c               unitary, R1 is upper trapezoidal, and 
-c               Q1*R1 = [A(:,1:j-1) A(:,j+1:n)], where A = Q*R.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R.
-c Q (io)        on entry, the unitary m-by-k matrix Q.
-c               on exit, the updated matrix Q1.
-c R (in)        the original upper trapezoidal matrix R.
-c R1 (out)      the updated matrix R1.
-c j (in)        the position of the deleted column in R.
-c               1 <= j <= n.
-c
-      integer m,n,k,j
-      complex Q(m,k),R(k,n),R1(k,n-1)
-      external xerbla,ccopy,cqhqr
-      integer info
-c quick return if possible      
-      if (m <= 0 .or. k <= 0 .or. n == 1) return
-c check arguments      
-      info = 0
-      if (n < 1) then
-        info = 2
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('CQRDEC',info)
-      end if
-c copy leading portion
-      call ccopy(k*(j-1),R,1,R1,1)
-      if (j < n) then
-c copy trailing portion of R        
-        call ccopy(k*(n-j),R(1,j+1),1,R1(1,j),1)
-c if necessary, retriangularize R1(j:k,j:n-1) and update Q(:,j:k)
-        if (j < k) then
-          call cqhqr(m,n-j,k-j+1,Q(1,j),m,R1(j,j),k)
-        end if
-      end if
-      end 

--- a/libcruft/qrupdate/cqrder.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,93 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqrder(m,n,Q,Q1,R,R1,j)
-c purpose:      updates a QR factorization after deleting a row.      
-c               i.e., given an m-by-m unitary matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:m, this subroutine forms the (m-1)-by-(m-1) matrix 
-c               Q1 and an (m-1)-by-n matrix R1 so that Q1 is again 
-c               unitary, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R.
-c               (complex version)
-c               
-c arguments:
-c m (in)        number of rows of the matrix R. 
-c n (in)        number of columns of the matrix R
-c Q (in)        the unitary matrix Q
-c Q1 (out)      the updated matrix Q1
-c R (in)        the upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new row in R1
-c
-      integer m,n,j
-      complex Q(m,m),Q1(m-1,m-1),R(m,n),R1(m-1,n)
-      real c
-      complex s,rr,w
-      external xerbla,clacpy,clartg,crot,csscal,caxpy
-      integer i
-c quick return if possible      
-      if (m == 1) return
-c check arguments      
-      info = 0
-      if (m < 1) then
-        info = 1
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('CQRDER',info)
-      end if
-c setup the new matrix Q1
-c permute the columns of Q and rows of R so that the deleted row ends 
-c up being the topmost row.      
-      if (j > 1) then
-        call clacpy('0',j-1,m-1,Q(1,2),m,Q1(1,1),m-1)
-      end if
-      if (j < m) then
-        call clacpy('0',m-j,m-1,Q(j+1,2),m,Q1(j,1),m-1)
-      end if
-c setup the new matrix R1
-      call clacpy('0',m-1,n,R(2,1),m,R1(1,1),m-1)
-c eliminate Q(j,2:m)
-      w = Q(j,m)
-      do i = m-1,2,-1
-        call clartg(Q(j,i),w,c,s,rr)
-        w = rr
-c apply rotation to rows of R1
-        if (i <= n) then
-          call crot(n-i+1,R1(i-1,i),m-1,R1(i,i),m-1,c,conjg(s))
-        end if
-c apply rotation to columns of Q1
-        call crot(m-1,Q1(1,i-1),1,Q1(1,i),1,c,s)
-      end do
-c the last iteration is special, as we don't have the first row of
-c R and first column of Q
-      call clartg(Q(j,1),w,c,s,rr)
-      w = rr
-      call csscal(n,c,R1(1,1),m-1)
-      call caxpy(n,-s,R(1,1),m,R1(1,1),m-1)
-c apply rotation to columns of Q1
-      call csscal(m-1,c,Q1(1,1),1)
-      if (j > 1) then
-        call caxpy(j-1,-conjg(s),Q(1,1),1,Q1(1,1),1)
-      end if
-      if (j < m) then
-        call caxpy(m-j,-conjg(s),Q(j+1,1),1,Q1(j,1),1)
-      end if
-      end 

--- a/libcruft/qrupdate/cqrinc.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,74 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqrinc(m,n,k,Q,R,R1,j,x)
-c purpose:      updates a QR factorization after inserting a new
-c               column.      
-c               i.e., given an m-by-k unitary matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               forms an m-by-(n+1) matrix R1 so that Q1 is again unitary,
-c               R1 upper trapezoidal, and 
-c               Q1*R1 = [A(:,1:j-1); Q*Q'*x; A(:,j:n-1)], where A = Q*R.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R. k <= m.
-c Q (io)        on entry, the unitary matrix Q.
-c               on exit, the updated matrix Q1
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new column in R1
-c x (in)        the column being inserted
-c
-      integer m,n,k,j
-      complex Q(m,k),R(k,n),R1(k,n+1),x(m)
-      complex w
-      external xerbla,ccopy,cqrqhv,cgemv
-      integer info,i,jj
-c quick return if possible      
-      if (m <= 0) return
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 2
-      else if (j < 1 .or. j > n+1) then
-        info = 6
-      end if
-      if (info /= 0) then
-        call xerbla('CQRINC',info)
-      end if
-c copy leading portion of R 
-      call ccopy(k*(j-1),R,1,R1,1)
-      if (j <= n) then
-        call ccopy(k*(n+1-j),R(1,j),1,R1(1,j+1),1)
-      end if
-      call cgemv('C',m,min(k,j-1),cmplx(1e0),Q,m,x,1,
-     +           cmplx(0e0),R1(1,j),1)
-      if (j < k) then
-c eliminate tail, updating Q(:,j:k) and R1(j:k,j+1:n+1)
-        jj = min(j,n)+1
-        call cqrqhv(m,n+1-j,k-j+1,Q(1,j),m,R1(j,jj),m,x,w)
-c assemble inserted column        
-        R1(j,j) = w
-        do i = j+1,k
-          R1(i,j) = 0e0
-        end do 
-      end if
-      end 

--- a/libcruft/qrupdate/cqrinr.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,73 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqrinr(m,n,Q,Q1,R,R1,j,x)
-c purpose:      updates a QR factorization after inserting a new
-c               row.      
-c               i.e., given an m-by-m unitary matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:m+1, this subroutine forms the (m+1)-by-(m+1) matrix 
-c               Q1 and an (m+1)-by-n matrix R1 so that Q1 is again 
-c               unitary, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix R. 
-c n (in)        number of columns of the matrix R
-c Q (in)        the orthogonal matrix Q
-c Q1 (out)      the updated matrix Q1
-c R (in)        the upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new row in R1
-c x (in)        the row being added
-c
-      integer m,n,j
-      complex Q(m,m),Q1(m+1,m+1),R(m,n),R1(m+1,n),x(n)
-      external xerbla,clacpy,ccopy,cqhqr
-      integer i
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 2
-      else if (j < 1 .or. j > m+1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('CQRINR',info)
-      end if
-c setup the new matrix Q1
-c permute the columns of Q1 and rows of R1 so that c the new row ends 
-c up being the topmost row.      
-      if (j > 1) then
-        call clacpy('0',j-1,m,Q(1,1),m,Q1(1,2),m+1)
-      end if
-      if (j <= m) then
-        call clacpy('0',m-j+1,m,Q(j,1),m,Q1(j+1,2),m+1)
-      end if
-c zero the rest of Q1      
-      do i = 1,m+1
-        Q1(i,1) = 0e0
-        Q1(j,i) = 0e0
-      end do
-      Q1(j,1) = 1e0
-c setup the new matrix R1
-      call ccopy(n,x,1,R1(1,1),m+1)
-      call clacpy('0',m,n,R(1,1),m,R1(2,1),m+1)
-c rotate to form proper QR      
-      call cqhqr(m+1,n,m+1,Q1,m+1,R1,m+1)
-      end 

--- a/libcruft/qrupdate/cqrqhu.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqrqhu(m,n,k,Q,ldq,R,ldr,u,rr)
-c purpose:      given an m-by-k matrix Q, an upper trapezoidal 
-c               k-by-n matrix R, and a k-vector u, 
-c               this subroutine updates the matrices Q -> Q1 and 
-c               R -> R1 so that Q1 = Q*G', R1 = G*R, u1(2:k) = 0 
-c               with G unitary, R1 upper Hessenberg, and u1 = G*u.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q and rows of R.
-c Q (io)        on entry, the unitary matrix Q.
-c               on exit, the updated matrix Q1.
-c ldq (in)      leading dimension of Q.
-c R (io)        on entry, the upper triangular matrix R.
-c               on exit, the updated upper Hessenberg matrix R1.
-c ldr (in)      leading dimension of R.
-c u (in)        the k-vector u.
-c rr (out)      the first element of Q1'*u on exit.
-c
-c               if Q is unitary, so is Q1. It is not strictly
-c               necessary, however.
-      integer m,n,k,ldq,ldr
-      complex Q(ldq,*),R(ldr,*),u(*),rr
-      real c
-      complex s,w
-      external xerbla,clartg,crot
-      integer i,info
-c quick return if possible.
-      if (k <= 0) return
-c check arguments.      
-      info = 0
-      if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('CQRQHU',info)
-      end if
-      rr = u(k)
-      do i = k-1,1,-1
-        w = rr
-        if (w /= cmplx(0e0,0e0)) then
-          call clartg(u(i),w,c,s,rr)
-c apply rotation to rows of R if necessary        
-          if (i <= n) then
-            call crot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s)
-          end if
-c apply rotation to columns of Q if necessary
-          if (m > 0) then
-            call crot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s))
-          end if
-        else
-c no rotation necessary
-          rr = u(i)
-        end if          
-      end do
-      end 
-

--- a/libcruft/qrupdate/cqrqhv.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,75 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqrqhv(m,n,k,Q,ldq,R,ldr,u,rr)
-c purpose:      given an m-by-k matrix Q, an upper trapezoidal 
-c               k-by-n matrix R, and an m-vector u, this subroutine 
-c               updates the matrices Q -> Q1 and R -> R1 so that 
-c               Q1 = Q*G', R1 = G*R, w1(2:m) = 0 with G unitary, 
-c               R1 upper Hessenberg, and w1 = Q1'*u.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q and rows of R. k <= m.
-c Q (io)        on entry, the unitary matrix Q.
-c               on exit, the updated matrix Q1.
-c ldq (in)      leading dimension of Q.
-c R (io)        on entry, the upper triangular matrix R.
-c               on exit, the updated upper Hessenberg matrix R1.
-c ldr (in)      leading dimension of R.
-c u (in)        the m-vector u.
-c rr (out)      the first element of Q1'*u on exit.
-c
-c               if Q is unitary, so is Q1. It is not strictly
-c               necessary, however.
-      integer m,n,k,ldq,ldr
-      complex Q(ldq,*),R(ldr,*),u(*),rr
-      real c
-      complex s,w,w1,cdotc
-      external xerbla,cdotc,clartg,crot
-      integer i,info
-c quick return if possible.
-      if (k <= 0) return
-c check arguments.      
-      info = 0
-      if (k > m) then
-        info = 3
-      else if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('CQRQHV',info)
-      end if
-c form each element of w = Q'*u when necessary.
-      rr = cdotc(m,Q(1,k),1,u,1)
-      do i = k-1,1,-1
-        w1 = rr
-        w = cdotc(m,Q(1,i),1,u,1)
-        call clartg(w,w1,c,s,rr)
-c apply rotation to rows of R if necessary        
-        if (i <= n) then
-          call crot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s)
-        end if
-c apply rotation to columns of Q
-        call crot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s))
-      end do
-      end 
-

--- a/libcruft/qrupdate/cqrshc.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,97 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine cqrshc(m,n,k,Q,R,i,j)
-c purpose:      updates a QR factorization after circular shift of
-c               columns.      
-c               i.e., given an m-by-k unitary matrix Q, an k-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               R -> R1 so that Q1 is again unitary, R1 upper trapezoidal, 
-c               and 
-c               Q1*R1 = A(:,p), where A = Q*R and p is the permutation
-c               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
-c               [1:j-1,shift(j:i,+1),i+1:n] if j > i.
-c               if m == 0, the matrix Q is ignored.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q, or 0 if Q is not needed.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R.
-c Q (io)        on entry, the (unitary) matrix Q.
-c               on exit, the updated matrix Q1
-c R (io)        on entry, the upper trapezoidal m-by-n matrix R.
-c               on exit, the updated matrix R1.
-c i (in)        the first index determining the range (see above)
-c j (in)        the second index determining the range (see above)
-c
-      integer m,n,k,i,j
-      complex Q(m,k),R(k,n)
-      external xerbla,cswap,cqhqr,cqrqhu
-      complex w
-      integer l,jj,kk,info
-
-c quick return if possible
-      if (k <= 0 .or. n <= 1) return
-      info = 0
-      if (m /= 0 .and. k > m) then
-        info = 3
-      else if (i < 1 .or. i > n) then
-        info = 6
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('CQRSHC',info)
-      end if
-
-      if (i < j) then
-c shift columns
-        do l = i,j-1
-          call cswap(min(k,l+1),R(1,l),1,R(1,l+1),1)
-        end do
-c retriangularize
-        if (i < k) then
-          kk = min(k,j)
-          if (m > 0) then
-            call cqhqr(m,n+1-i,kk+1-i,Q(1,i),m,R(i,i),k)
-          else
-            call cqhqr(0,n+1-i,kk+1-i,Q,1,R(i,i),k)
-          endif
-        end if
-      else if (j < i) then
-c shift columns
-        do l = i,j+1,-1
-          call cswap(min(k,i),R(1,l),1,R(1,l-1),1)
-        end do
-c retriangularize
-        if (j < k) then
-          jj = min(j+1,n)
-          kk = min(k,i)
-          if (m > 0) then
-            call cqrqhu(m,n-j,kk+1-j,Q(1,j),m,R(j,jj),k,R(j,j),w)
-          else
-            call cqrqhu(0,n-j,kk+1-j,Q,1,R(j,jj),k,R(j,j),w)
-          end if
-          R(j,j) = w
-          do jj = j+1,kk
-            R(jj,j) = 0
-          end do
-        end if
-      end if
-      end

--- a/libcruft/qrupdate/dch1dn.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dch1dn(n,R,u,w,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine downdates R -> R1 so that
-c               R1'*R1 = A - u*u' 
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n,info
-      double precision R(n,n),u(n)
-      double precision w(n)
-      external dtrsv,dlartg,dnrm2
-      double precision rho,dnrm2
-      double precision rr,ui,t
-      integer i,j
-
-c quick return if possible
-      if (n <= 0) return
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0d0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call dtrsv('U','T','N',n,R,n,u,1)
-      rho = dnrm2(n,u,1)
-c check positive definiteness      
-      rho = 1 - rho**2
-      if (rho <= 0d0) then
-        info = 1
-        return
-      end if
-      rho = sqrt(rho)
-c eliminate R' \ u
-      do i = n,1,-1
-        ui = u(i)
-c generate next rotation        
-        call dlartg(rho,ui,w(i),u(i),rr)
-        rho = rr
-      end do
-c apply rotations
-      do i = n,1,-1
-        ui = 0d0
-        do j = i,1,-1
-          t = w(j)*ui + u(j)*R(j,i)
-          R(j,i) = w(j)*R(j,i) - u(j)*ui
-          ui = t
-        end do
-      end do
-
-      info = 0
-      end

--- a/libcruft/qrupdate/dch1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,57 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine dch1up(n,R,u,w)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A + u*u'
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n
-      double precision R(n,n),u(n)
-      double precision w(n)
-      external dlartg
-      double precision rr,ui,t
-      integer i,j
-      
-      do i = 1,n
-c apply stored rotations, column-wise
-        ui = u(i)
-        do j = 1,i-1
-          t = w(j)*R(j,i) + u(j)*ui
-          ui = w(j)*ui - u(j)*R(j,i)
-          R(j,i) = t
-        end do
-c generate next rotation        
-        call dlartg(R(i,i),ui,w(i),u(i),rr)
-        R(i,i) = rr
-      end do
-      end
-

--- a/libcruft/qrupdate/dchdex.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,61 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine dchdex(n,R,R1,j)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the deleted row/column
-      integer n,j,info
-      double precision R(n,n),R1(n-1,n-1)
-      double precision Qdum,c,s,rr
-      external xerbla,dlacpy,dqhqr,dlartg
-
-c quick return if possible
-      if (n == 1) return
-
-c check arguments      
-      info = 0
-      if (n <= 0) then
-        info = 1
-      else if (j < 1 .or. j > n) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('DCHDEX',info)
-      end if
-
-c setup the new matrix R1
-      if (j > 1) then
-        call dlacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1)
-      end if
-      if (j < n) then
-        call dlacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1)
-        call dqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1)
-c eliminate R(n,n)      
-        call dlartg(R1(n-1,n-1),R(n,n),c,s,rr)
-        R1(n-1,n-1) = rr
-      endif
-      end

--- a/libcruft/qrupdate/dchinx.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,108 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine dchinx(n,R,R1,j,u,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
-c               jj = [1:j-1,j+1:n+1].
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the inserted row/column
-c u (in)        the vector (n+1) determining the rank-1 update
-c info (out)    on exit, if info = 1, the 
-c               definiteness.
-
-      integer n,j,info
-      double precision R(n,n),R1(n+1,n+1),u(n+1)
-      double precision rho,Qdum,w,dnrm2
-      external xerbla,dcopy,dlacpy,dtrsv,dnrm2,dqrqhu
-      integer jj
-
-c quick return if possible
-      if (n == 0) then
-        if (u(1) <= 0) then
-          info = 1
-          return
-        else
-          R(1,1) = sqrt(u(1))
-        end if
-      end if
-
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 1
-      else if (j < 1 .or. j > n+1) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('DCHINX',info)
-      end if
-
-c copy shifted vector
-      if (j > 1) then
-        call dcopy(j-1,u,1,R1(1,j),1)
-      end if
-      w = u(j)
-      if (j < n+1) then
-        call dcopy(n-j+1,u(j+1),1,R1(j,j),1)
-      end if
-      
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0d0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call dtrsv('U','T','N',n,R,n,R1(1,j),1)
-      rho = dnrm2(n,R1(1,j),1)
-c check positive definiteness      
-      rho = u(j) - rho**2
-      if (rho <= 0d0) then
-        info = 1
-        return
-      end if
-      R1(n+1,n+1) = sqrt(rho)
-
-c setup the new matrix R1
-      do i = 1,n+1
-        R1(n+1,i) = 0d0
-      end do
-      if (j > 1) then
-        call dlacpy('0',n,j-1,R(1,1),n,R1(1,1),n+1)
-      end if
-      if (j <= n) then
-        call dlacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1)
-c retriangularize
-        jj = min(j+1,n)
-        call dqrqhu(0,n+1-j,n-j,Qdum,1,R1(j,jj),n+1,R1(j,j),w)
-        R1(j,j) = w
-        do jj = j+1,n
-          R1(jj,j) = 0d0
-        end do
-      end if
-
-      end

--- a/libcruft/qrupdate/dqhqr.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqhqr(m,n,k,Q,ldq,R,ldr)
-c purpose:      given an k-by-n upper Hessenberg matrix R and
-c               an m-by-k matrix Q, this subroutine updates
-c               R -> R1 and Q -> Q1 so that R1 is upper 
-c               trapezoidal, R1 = G*R and Q1 = Q*G', where
-c               G is an orthogonal matrix, giving Q1*R1 = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q
-c n (in)        number of columns of the matrix R
-c k (in)        number of columns of Q and rows of R.
-c Q (io)        on entry, the orthogonal matrix Q
-c               on exit, the updated matrix Q1
-c ldq (in)      leading dimension of Q
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated upper Hessenberg matrix R1
-c ldr (in)      leading dimension of R
-c
-      integer m,n,k,ldq,ldr
-      double precision Q(ldq,*),R(ldr,*)
-      double precision c
-      double precision s,rr
-      external xerbla,dlartg,drot
-      integer info,i
-c quick return if possible.
-      if (n <= 0 .or. k <= 1) return
-c check arguments.
-      info = 0
-      if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('DQHQR',info)
-      end if
-c triangularize      
-      do i = 1,min(k-1,n)
-        call dlartg(R(i,i),R(i+1,i),c,s,rr)
-        R(i,i) = rr
-        R(i+1,i) = 0d0
-        if (i < n) then
-          call drot(n-i,R(i,i+1),ldr,R(i+1,i+1),ldr,c,s)
-        end if
-c apply rotation to Q        
-        if (m > 0) then
-          call drot(m,Q(1,i),1,Q(1,i+1),1,c,s)
-        end if
-      end do
-      end 
-

--- a/libcruft/qrupdate/dqr1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,52 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqr1up(m,n,k,Q,R,u,v)
-c purpose:      updates a QR factorization after rank-1 modification
-c               i.e., given a m-by-k orthogonal Q and m-by-n upper 
-c               trapezoidal R, an m-vector u and n-vector v, 
-c               this subroutine updates Q -> Q1 and R -> R1 so that
-c               Q1*R1 = Q*R + Q*Q'u*v', and Q1 is again orthonormal
-c               and R1 upper trapezoidal.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R. k <= m.
-c Q (io)        on entry, the orthogonal m-by-k matrix Q.
-c               on exit, the updated matrix Q1.
-c R (io)        on entry, the upper trapezoidal m-by-n matrix R..
-c               on exit, the updated matrix R1.
-c u (in)        the left m-vector.
-c v (in)        the right n-vector.
-c
-      integer m,n,k
-      double precision Q(m,k),R(k,n),u(m),v(n)
-      double precision w
-      external dqrqhv,dqhqr,daxpy
-c quick return if possible      
-      if (m <= 0 .or. n <= 0) return
-c eliminate tail of Q'*u
-      call dqrqhv(m,n,k,Q,m,R,m,u,w)
-c update R      
-
-      call daxpy(n,w,v,1,R(1,1),m)
-
-c retriangularize R
-      call dqhqr(m,n,k,Q,m,R,k)
-      end 

--- a/libcruft/qrupdate/dqrdec.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,66 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqrdec(m,n,k,Q,R,R1,j)
-c purpose:      updates a QR factorization after deleting
-c               a column.      
-c               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               forms an m-by-(n-1) matrix R1 so that Q1 remains
-c               orthogonal, R1 is upper trapezoidal, and 
-c               Q1*R1 = [A(:,1:j-1) A(:,j+1:n)], where A = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R.
-c Q (io)        on entry, the orthogonal m-by-k matrix Q.
-c               on exit, the updated matrix Q1.
-c R (in)        the original upper trapezoidal matrix R.
-c R1 (out)      the updated matrix R1.
-c j (in)        the position of the deleted column in R.
-c               1 <= j <= n.
-c
-      integer m,n,k,j
-      double precision Q(m,k),R(k,n),R1(k,n-1)
-      external xerbla,dcopy,dqhqr
-      integer info
-c quick return if possible      
-      if (m <= 0 .or. k <= 0 .or. n == 1) return
-c check arguments      
-      info = 0
-      if (n < 1) then
-        info = 2
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('DQRDEC',info)
-      end if
-c copy leading portion
-      call dcopy(k*(j-1),R,1,R1,1)
-      if (j < n) then
-c copy trailing portion of R        
-        call dcopy(k*(n-j),R(1,j+1),1,R1(1,j),1)
-c if necessary, retriangularize R1(j:k,j:n-1) and update Q(:,j:k)
-        if (j < k) then
-          call dqhqr(m,n-j,k-j+1,Q(1,j),m,R1(j,j),k)
-        end if
-      end if
-      end

--- a/libcruft/qrupdate/dqrder.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,93 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqrder(m,n,Q,Q1,R,R1,j)
-c purpose:      updates a QR factorization after deleting a row.      
-c               i.e., given an m-by-m orthogonal matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:m, this subroutine forms the (m-1)-by-(m-1) matrix 
-c               Q1 and an (m-1)-by-n matrix R1 so that Q1 is again 
-c               orthogonal, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R.
-c               (real version)
-c               
-c arguments:
-c m (in)        number of rows of the matrix R. 
-c n (in)        number of columns of the matrix R
-c Q (in)        the orthogonal matrix Q
-c Q1 (out)      the updated matrix Q1
-c R (in)        the upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new row in R1
-c
-      integer m,n,j
-      double precision Q(m,m),Q1(m-1,m-1),R(m,n),R1(m-1,n)
-      double precision c
-      double precision s,rr,w
-      external xerbla,dlacpy,dlartg,drot,dscal,daxpy
-      integer i
-c quick return if possible      
-      if (m == 1) return
-c check arguments      
-      info = 0
-      if (m < 1) then
-        info = 1
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('DQRDER',info)
-      end if
-c setup the new matrix Q1
-c permute the columns of Q and rows of R so that the deleted row ends 
-c up being the topmost row.      
-      if (j > 1) then
-        call dlacpy('0',j-1,m-1,Q(1,2),m,Q1(1,1),m-1)
-      end if
-      if (j < m) then
-        call dlacpy('0',m-j,m-1,Q(j+1,2),m,Q1(j,1),m-1)
-      end if
-c setup the new matrix R1
-      call dlacpy('0',m-1,n,R(2,1),m,R1(1,1),m-1)
-c eliminate Q(j,2:m)
-      w = Q(j,m)
-      do i = m-1,2,-1
-        call dlartg(Q(j,i),w,c,s,rr)
-        w = rr
-c apply rotation to rows of R1
-        if (i <= n) then
-          call drot(n-i+1,R1(i-1,i),m-1,R1(i,i),m-1,c,s)
-        end if
-c apply rotation to columns of Q1
-        call drot(m-1,Q1(1,i-1),1,Q1(1,i),1,c,s)
-      end do
-c the last iteration is special, as we don't have the first row of
-c R and first column of Q
-      call dlartg(Q(j,1),w,c,s,rr)
-      w = rr
-      call dscal(n,c,R1(1,1),m-1)
-      call daxpy(n,-s,R(1,1),m,R1(1,1),m-1)
-c apply rotation to columns of Q1
-      call dscal(m-1,c,Q1(1,1),1)
-      if (j > 1) then
-        call daxpy(j-1,-s,Q(1,1),1,Q1(1,1),1)
-      end if
-      if (j < m) then
-        call daxpy(m-j,-s,Q(j+1,1),1,Q1(j,1),1)
-      end if
-      end 

--- a/libcruft/qrupdate/dqrinc.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,73 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqrinc(m,n,k,Q,R,R1,j,x)
-c purpose:      updates a QR factorization after inserting a new
-c               column.      
-c               i.e., given an m-by-k orthogonal matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               forms an m-by-(n+1) matrix R1 so that Q1 is again
-c               orthogonal, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(:,1:j-1); Q*Q'*x; A(:,j:n-1)], where A = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R. k <= m.
-c Q (io)        on entry, the orthogonal matrix Q.
-c               on exit, the updated matrix Q1
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new column in R1
-c x (in)        the column being inserted
-c
-      integer m,n,k,j
-      double precision Q(m,k),R(k,n),R1(k,n+1),x(m)
-      double precision w
-      external xerbla,dcopy,dqrqhv,dgemv
-      integer info,i,jj
-c quick return if possible      
-      if (m <= 0) return
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 2
-      else if (j < 1 .or. j > n+1) then
-        info = 6
-      end if
-      if (info /= 0) then
-        call xerbla('DQRINC',info)
-      end if
-c copy leading portion of R 
-      call dcopy(k*(j-1),R,1,R1,1)
-      if (j <= n) then
-        call dcopy(k*(n+1-j),R(1,j),1,R1(1,j+1),1)
-      end if
-      call dgemv('T',m,min(k,j-1),1d0,Q,m,x,1,0d0,R1(1,j),1)
-      if (j < k) then
-c eliminate tail, updating Q(:,j:k) and R1(j:k,j+1:n+1)
-        jj = min(j,n)+1
-        call dqrqhv(m,n+1-j,k-j+1,Q(1,j),m,R1(j,jj),m,x,w)
-c assemble inserted column        
-        R1(j,j) = w
-        do i = j+1,k
-          R1(i,j) = 0d0
-        end do 
-      end if
-      end 

--- a/libcruft/qrupdate/dqrinr.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,73 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqrinr(m,n,Q,Q1,R,R1,j,x)
-c purpose:      updates a QR factorization after inserting a new
-c               row.      
-c               i.e., given an m-by-m orthogonal matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:m+1, this subroutine forms the (m+1)-by-(m+1) matrix 
-c               Q1 and an (m+1)-by-n matrix R1 so that Q1 is again 
-c               orthogonal, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix R. 
-c n (in)        number of columns of the matrix R
-c Q (in)        the orthogonal matrix Q
-c Q1 (out)      the updated matrix Q1
-c R (in)        the upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new row in R1
-c x (in)        the row being added
-c
-      integer m,n,j
-      double precision Q(m,m),Q1(m+1,m+1),R(m,n),R1(m+1,n),x(n)
-      external xerbla,dlacpy,dcopy,dqhqr
-      integer i
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 2
-      else if (j < 1 .or. j > m+1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('DQRINR',info)
-      end if
-c setup the new matrix Q1
-c permute the columns of Q1 and rows of R1 so that c the new row ends 
-c up being the topmost row.      
-      if (j > 1) then
-        call dlacpy('0',j-1,m,Q(1,1),m,Q1(1,2),m+1)
-      end if
-      if (j <= m) then
-        call dlacpy('0',m-j+1,m,Q(j,1),m,Q1(j+1,2),m+1)
-      end if
-c zero the rest of Q1      
-      do i = 1,m+1
-        Q1(i,1) = 0d0
-        Q1(j,i) = 0d0
-      end do
-      Q1(j,1) = 1d0
-c setup the new matrix R1
-      call dcopy(n,x,1,R1(1,1),m+1)
-      call dlacpy('0',m,n,R(1,1),m,R1(2,1),m+1)
-c rotate to form proper QR      
-      call dqhqr(m+1,n,m+1,Q1,m+1,R1,m+1)
-      end 

--- a/libcruft/qrupdate/dqrqhu.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqrqhu(m,n,k,Q,ldq,R,ldr,u,rr)
-c purpose:      given an m-by-k matrix Q, an upper trapezoidal 
-c               k-by-n matrix R, and a k-vector u, 
-c               this subroutine updates the matrices Q -> Q1 and 
-c               R -> R1 so that Q1 = Q*G', R1 = G*R, u1(2:k) = 0 
-c               with G orthogonal, R1 upper Hessenberg, and u1 = G*u.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q and rows of R.
-c Q (io)        on entry, the orthogonal matrix Q.
-c               on exit, the updated matrix Q1.
-c ldq (in)      leading dimension of Q.
-c R (io)        on entry, the upper triangular matrix R.
-c               on exit, the updated upper Hessenberg matrix R1.
-c ldr (in)      leading dimension of R.
-c u (in)        the k-vector u.
-c rr (out)      the first element of Q1'*u on exit.
-c
-c               if Q is orthogonal, so is Q1. It is not strictly
-c               necessary, however.
-      integer m,n,k,ldq,ldr
-      double precision Q(ldq,*),R(ldr,*),u(*),rr
-      double precision c
-      double precision s,w
-      external xerbla,dlartg,drot
-      integer i,info
-c quick return if possible.
-      if (k <= 0) return
-c check arguments.      
-      info = 0
-      if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('DQRQHU',info)
-      end if
-      rr = u(k)
-      do i = k-1,1,-1
-        w = rr
-        if (w /= 0d0) then
-          call dlartg(u(i),w,c,s,rr)
-c apply rotation to rows of R if necessary        
-          if (i <= n) then
-            call drot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s)
-          end if
-c apply rotation to columns of Q if necessary
-          if (m > 0) then
-            call drot(m,Q(1,i),1,Q(1,i+1),1,c,s)
-          end if
-        else
-c no rotation necessary
-          rr = u(i)
-        end if          
-      end do
-      end 
-

--- a/libcruft/qrupdate/dqrqhv.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,75 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqrqhv(m,n,k,Q,ldq,R,ldr,u,rr)
-c purpose:      given an m-by-k matrix Q, an upper trapezoidal 
-c               k-by-n matrix R, and an m-vector u, this subroutine 
-c               updates the matrices Q -> Q1 and R -> R1 so that 
-c               Q1 = Q*G', R1 = G*R, w1(2:m) = 0 with G orthogonal, 
-c               R1 upper Hessenberg, and w1 = Q1'*u.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q and rows of R. k <= m.
-c Q (io)        on entry, the orthogonal matrix Q.
-c               on exit, the updated matrix Q1.
-c ldq (in)      leading dimension of Q.
-c R (io)        on entry, the upper triangular matrix R.
-c               on exit, the updated upper Hessenberg matrix R1.
-c ldr (in)      leading dimension of R.
-c u (in)        the m-vector u.
-c rr (out)      the first element of Q1'*u on exit.
-c
-c               if Q is orthogonal, so is Q1. It is not strictly
-c               necessary, however.
-      integer m,n,k,ldq,ldr
-      double precision Q(ldq,*),R(ldr,*),u(*),rr
-      double precision c
-      double precision s,w,w1,ddot
-      external xerbla,ddot,dlartg,drot
-      integer i,info
-c quick return if possible.
-      if (k <= 0) return
-c check arguments.      
-      info = 0
-      if (k > m) then
-        info = 3
-      else if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('DQRQHV',info)
-      end if
-c form each element of w = Q'*u when necessary.
-      rr = ddot(m,Q(1,k),1,u,1)
-      do i = k-1,1,-1
-        w1 = rr
-        w = ddot(m,Q(1,i),1,u,1)
-        call dlartg(w,w1,c,s,rr)
-c apply rotation to rows of R if necessary        
-        if (i <= n) then
-          call drot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s)
-        end if
-c apply rotation to columns of Q
-        call drot(m,Q(1,i),1,Q(1,i+1),1,c,s)
-      end do
-      end 
-

--- a/libcruft/qrupdate/dqrshc.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,97 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine dqrshc(m,n,k,Q,R,i,j)
-c purpose:      updates a QR factorization after circular shift of
-c               columns.      
-c               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               R -> R1 so that Q1 is again orthogonal, R1 upper 
-c               trapezoidal, and 
-c               Q1*R1 = A(:,p), where A = Q*R and p is the permutation
-c               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
-c               [1:j-1,shift(j:i,+1),i+1:n] if j > i.
-c               if m == 0, the matrix Q is ignored.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q, or 0 if Q is not needed.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R.
-c Q (io)        on entry, the (orthogonal) matrix Q.
-c               on exit, the updated matrix Q1
-c R (io)        on entry, the upper trapezoidal m-by-n matrix R.
-c               on exit, the updated matrix R1.
-c i (in)        the first index determining the range (see above)
-c j (in)        the second index determining the range (see above)
-c
-      integer m,n,k,i,j
-      double precision Q(m,k),R(k,n)
-      external xerbla,dswap,dqhqr,dqrqhu
-      double precision w
-      integer l,jj,kk,info
-
-c quick return if possible
-      if (k <= 0 .or. n <= 1) return
-      info = 0
-      if (m /= 0 .and. k > m) then
-        info = 3
-      else if (i < 1 .or. i > n) then
-        info = 6
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('DQRSHC',info)
-      end if
-
-      if (i < j) then
-c shift columns
-        do l = i,j-1
-          call dswap(min(k,l+1),R(1,l),1,R(1,l+1),1)
-        end do
-c retriangularize
-        if (i < k) then
-          kk = min(k,j)
-          if (m > 0) then
-            call dqhqr(m,n+1-i,kk+1-i,Q(1,i),m,R(i,i),k)
-          else
-            call dqhqr(0,n+1-i,kk+1-i,Q,1,R(i,i),k)
-          endif
-        end if
-      else if (j < i) then
-c shift columns
-        do l = i,j+1,-1
-          call dswap(min(k,i),R(1,l),1,R(1,l-1),1)
-        end do
-c retriangularize
-        if (j < k) then
-          jj = min(j+1,n)
-          kk = min(k,i)
-          if (m > 0) then
-            call dqrqhu(m,n-j,kk+1-j,Q(1,j),m,R(j,jj),k,R(j,j),w)
-          else
-            call dqrqhu(0,n-j,kk+1-j,Q,1,R(j,jj),k,R(j,j),w)
-          end if
-          R(j,j) = w
-          do jj = j+1,kk
-            R(jj,j) = 0
-          end do
-        end if
-      end if
-      end

--- a/libcruft/qrupdate/sch1dn.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sch1dn(n,R,u,w,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine downdates R -> R1 so that
-c               R1'*R1 = A - u*u' 
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n,info
-      real R(n,n),u(n)
-      real w(n)
-      external strsv,slartg,snrm2
-      real rho,snrm2
-      real rr,ui,t
-      integer i,j
-
-c quick return if possible
-      if (n <= 0) return
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0e0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call strsv('U','T','N',n,R,n,u,1)
-      rho = snrm2(n,u,1)
-c check positive definiteness      
-      rho = 1 - rho**2
-      if (rho <= 0e0) then
-        info = 1
-        return
-      end if
-      rho = sqrt(rho)
-c eliminate R' \ u
-      do i = n,1,-1
-        ui = u(i)
-c generate next rotation        
-        call slartg(rho,ui,w(i),u(i),rr)
-        rho = rr
-      end do
-c apply rotations
-      do i = n,1,-1
-        ui = 0e0
-        do j = i,1,-1
-          t = w(j)*ui + u(j)*R(j,i)
-          R(j,i) = w(j)*R(j,i) - u(j)*ui
-          ui = t
-        end do
-      end do
-
-      info = 0
-      end

--- a/libcruft/qrupdate/sch1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,57 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine sch1up(n,R,u,w)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A + u*u'
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n
-      real R(n,n),u(n)
-      real w(n)
-      external slartg
-      real rr,ui,t
-      integer i,j
-      
-      do i = 1,n
-c apply stored rotations, column-wise
-        ui = u(i)
-        do j = 1,i-1
-          t = w(j)*R(j,i) + u(j)*ui
-          ui = w(j)*ui - u(j)*R(j,i)
-          R(j,i) = t
-        end do
-c generate next rotation        
-        call slartg(R(i,i),ui,w(i),u(i),rr)
-        R(i,i) = rr
-      end do
-      end
-

--- a/libcruft/qrupdate/schdex.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,61 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine schdex(n,R,R1,j)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the deleted row/column
-      integer n,j,info
-      real R(n,n),R1(n-1,n-1)
-      real Qdum,c,s,rr
-      external xerbla,slacpy,sqhqr,slartg
-
-c quick return if possible
-      if (n == 1) return
-
-c check arguments      
-      info = 0
-      if (n <= 0) then
-        info = 1
-      else if (j < 1 .or. j > n) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('SQRDEX',info)
-      end if
-
-c setup the new matrix R1
-      if (j > 1) then
-        call slacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1)
-      end if
-      if (j < n) then
-        call slacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1)
-        call sqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1)
-c eliminate R(n,n)      
-        call slartg(R1(n-1,n-1),R(n,n),c,s,rr)
-        R1(n-1,n-1) = rr
-      endif
-      end

--- a/libcruft/qrupdate/schinx.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,108 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine schinx(n,R,R1,j,u,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
-c               jj = [1:j-1,j+1:n+1].
-c               (real version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the inserted row/column
-c u (in)        the vector (n+1) determining the rank-1 update
-c info (out)    on exit, if info = 1, the 
-c               definiteness.
-
-      integer n,j,info
-      real R(n,n),R1(n+1,n+1),u(n+1)
-      real rho,Qdum,w,snrm2
-      external xerbla,scopy,slacpy,strsv,snrm2,sqrqhu
-      integer jj
-
-c quick return if possible
-      if (n == 0) then
-        if (u(1) <= 0) then
-          info = 1
-          return
-        else
-          R(1,1) = sqrt(u(1))
-        end if
-      end if
-
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 1
-      else if (j < 1 .or. j > n+1) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('SCHINX',info)
-      end if
-
-c copy shifted vector
-      if (j > 1) then
-        call scopy(j-1,u,1,R1(1,j),1)
-      end if
-      w = u(j)
-      if (j < n+1) then
-        call scopy(n-j+1,u(j+1),1,R1(j,j),1)
-      end if
-      
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0e0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call strsv('U','T','N',n,R,n,R1(1,j),1)
-      rho = snrm2(n,R1(1,j),1)
-c check positive definiteness      
-      rho = u(j) - rho**2
-      if (rho <= 0e0) then
-        info = 1
-        return
-      end if
-      R1(n+1,n+1) = sqrt(rho)
-
-c setup the new matrix R1
-      do i = 1,n+1
-        R1(n+1,i) = 0e0
-      end do
-      if (j > 1) then
-        call slacpy('0',n,j-1,R(1,1),n,R1(1,1),n+1)
-      end if
-      if (j <= n) then
-        call slacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1)
-c retriangularize
-        jj = min(j+1,n)
-        call sqrqhu(0,n+1-j,n-j,Qdum,1,R1(j,jj),n+1,R1(j,j),w)
-        R1(j,j) = w
-        do jj = j+1,n
-          R1(jj,j) = 0e0
-        end do
-      end if
-
-      end

--- a/libcruft/qrupdate/sqhqr.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqhqr(m,n,k,Q,ldq,R,ldr)
-c purpose:      given an k-by-n upper Hessenberg matrix R and
-c               an m-by-k matrix Q, this subroutine updates
-c               R -> R1 and Q -> Q1 so that R1 is upper 
-c               trapezoidal, R1 = G*R and Q1 = Q*G', where
-c               G is an orthogonal matrix, giving Q1*R1 = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q
-c n (in)        number of columns of the matrix R
-c k (in)        number of columns of Q and rows of R.
-c Q (io)        on entry, the orthogonal matrix Q
-c               on exit, the updated matrix Q1
-c ldq (in)      leading dimension of Q
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated upper Hessenberg matrix R1
-c ldr (in)      leading dimension of R
-c
-      integer m,n,k,ldq,ldr
-      real Q(ldq,*),R(ldr,*)
-      real c
-      real s,rr
-      external xerbla,slartg,srot
-      integer info,i
-c quick return if possible.
-      if (n <= 0 .or. k <= 1) return
-c check arguments.
-      info = 0
-      if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('SQHQR',info)
-      end if
-c triangularize      
-      do i = 1,min(k-1,n)
-        call slartg(R(i,i),R(i+1,i),c,s,rr)
-        R(i,i) = rr
-        R(i+1,i) = 0e0
-        if (i < n) then
-          call srot(n-i,R(i,i+1),ldr,R(i+1,i+1),ldr,c,s)
-        end if
-c apply rotation to Q        
-        if (m > 0) then
-          call srot(m,Q(1,i),1,Q(1,i+1),1,c,s)
-        end if
-      end do
-      end 
-

--- a/libcruft/qrupdate/sqr1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,52 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqr1up(m,n,k,Q,R,u,v)
-c purpose:      updates a QR factorization after rank-1 modification
-c               i.e., given a m-by-k orthogonal Q and m-by-n upper 
-c               trapezoidal R, an m-vector u and n-vector v, 
-c               this subroutine updates Q -> Q1 and R -> R1 so that
-c               Q1*R1 = Q*R + Q*Q'u*v', and Q1 is again orthonormal
-c               and R1 upper trapezoidal.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R. k <= m.
-c Q (io)        on entry, the orthogonal m-by-k matrix Q.
-c               on exit, the updated matrix Q1.
-c R (io)        on entry, the upper trapezoidal m-by-n matrix R..
-c               on exit, the updated matrix R1.
-c u (in)        the left m-vector.
-c v (in)        the right n-vector.
-c
-      integer m,n,k
-      real Q(m,k),R(k,n),u(m),v(n)
-      real w
-      external sqrqhv,sqhqr,saxpy
-c quick return if possible      
-      if (m <= 0 .or. n <= 0) return
-c eliminate tail of Q'*u
-      call sqrqhv(m,n,k,Q,m,R,m,u,w)
-c update R      
-
-      call saxpy(n,w,v,1,R(1,1),m)
-
-c retriangularize R
-      call sqhqr(m,n,k,Q,m,R,k)
-      end 

--- a/libcruft/qrupdate/sqrdec.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,66 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqrdec(m,n,k,Q,R,R1,j)
-c purpose:      updates a QR factorization after deleting
-c               a column.      
-c               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               forms an m-by-(n-1) matrix R1 so that Q1 remains
-c               orthogonal, R1 is upper trapezoidal, and 
-c               Q1*R1 = [A(:,1:j-1) A(:,j+1:n)], where A = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R.
-c Q (io)        on entry, the orthogonal m-by-k matrix Q.
-c               on exit, the updated matrix Q1.
-c R (in)        the original upper trapezoidal matrix R.
-c R1 (out)      the updated matrix R1.
-c j (in)        the position of the deleted column in R.
-c               1 <= j <= n.
-c
-      integer m,n,k,j
-      real Q(m,k),R(k,n),R1(k,n-1)
-      external xerbla,scopy,sqhqr
-      integer info
-c quick return if possible      
-      if (m <= 0 .or. k <= 0 .or. n == 1) return
-c check arguments      
-      info = 0
-      if (n < 1) then
-        info = 2
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('SQRDEC',info)
-      end if
-c copy leading portion
-      call scopy(k*(j-1),R,1,R1,1)
-      if (j < n) then
-c copy trailing portion of R        
-        call scopy(k*(n-j),R(1,j+1),1,R1(1,j),1)
-c if necessary, retriangularize R1(j:k,j:n-1) and update Q(:,j:k)
-        if (j < k) then
-          call sqhqr(m,n-j,k-j+1,Q(1,j),m,R1(j,j),k)
-        end if
-      end if
-      end

--- a/libcruft/qrupdate/sqrder.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,93 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqrder(m,n,Q,Q1,R,R1,j)
-c purpose:      updates a QR factorization after deleting a row.      
-c               i.e., given an m-by-m orthogonal matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:m, this subroutine forms the (m-1)-by-(m-1) matrix 
-c               Q1 and an (m-1)-by-n matrix R1 so that Q1 is again 
-c               orthogonal, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R.
-c               (real version)
-c               
-c arguments:
-c m (in)        number of rows of the matrix R. 
-c n (in)        number of columns of the matrix R
-c Q (in)        the orthogonal matrix Q
-c Q1 (out)      the updated matrix Q1
-c R (in)        the upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new row in R1
-c
-      integer m,n,j
-      real Q(m,m),Q1(m-1,m-1),R(m,n),R1(m-1,n)
-      real c
-      real s,rr,w
-      external xerbla,slacpy,slartg,srot,sscal,saxpy
-      integer i
-c quick return if possible      
-      if (m == 1) return
-c check arguments      
-      info = 0
-      if (m < 1) then
-        info = 1
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('SQRDER',info)
-      end if
-c setup the new matrix Q1
-c permute the columns of Q and rows of R so that the deleted row ends 
-c up being the topmost row.      
-      if (j > 1) then
-        call slacpy('0',j-1,m-1,Q(1,2),m,Q1(1,1),m-1)
-      end if
-      if (j < m) then
-        call slacpy('0',m-j,m-1,Q(j+1,2),m,Q1(j,1),m-1)
-      end if
-c setup the new matrix R1
-      call slacpy('0',m-1,n,R(2,1),m,R1(1,1),m-1)
-c eliminate Q(j,2:m)
-      w = Q(j,m)
-      do i = m-1,2,-1
-        call slartg(Q(j,i),w,c,s,rr)
-        w = rr
-c apply rotation to rows of R1
-        if (i <= n) then
-          call srot(n-i+1,R1(i-1,i),m-1,R1(i,i),m-1,c,s)
-        end if
-c apply rotation to columns of Q1
-        call srot(m-1,Q1(1,i-1),1,Q1(1,i),1,c,s)
-      end do
-c the last iteration is special, as we don't have the first row of
-c R and first column of Q
-      call slartg(Q(j,1),w,c,s,rr)
-      w = rr
-      call sscal(n,c,R1(1,1),m-1)
-      call saxpy(n,-s,R(1,1),m,R1(1,1),m-1)
-c apply rotation to columns of Q1
-      call sscal(m-1,c,Q1(1,1),1)
-      if (j > 1) then
-        call saxpy(j-1,-s,Q(1,1),1,Q1(1,1),1)
-      end if
-      if (j < m) then
-        call saxpy(m-j,-s,Q(j+1,1),1,Q1(j,1),1)
-      end if
-      end 

--- a/libcruft/qrupdate/sqrinc.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,75 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqrinc(m,n,k,Q,R,R1,j,x)
-c purpose:      updates a QR factorization after inserting a new
-c               column.      
-c               i.e., given an m-by-k orthogonal matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               forms an m-by-(n+1) matrix R1 so that Q1 is again
-c               orthogonal, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(:,1:j-1); Q*Q'*x; A(:,j:n-1)], where A = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R. k <= m.
-c Q (io)        on entry, the orthogonal matrix Q.
-c               on exit, the updated matrix Q1
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new column in R1
-c x (in)        the column being inserted
-c
-      integer m,n,k,j
-      real Q(m,k),R(k,n),R1(k,n+1),x(m)
-
-
-      real w
-      external xerbla,scopy,sqrqhv,sgemv
-      integer info,i,jj
-c quick return if possible      
-      if (m <= 0) return
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 2
-      else if (j < 1 .or. j > n+1) then
-        info = 6
-      end if
-      if (info /= 0) then
-        call xerbla('SQRINC',info)
-      end if
-c copy leading portion of R 
-      call scopy(k*(j-1),R,1,R1,1)
-      if (j <= n) then
-        call scopy(k*(n+1-j),R(1,j),1,R1(1,j+1),1)
-      end if
-      call sgemv('T',m,min(k,j-1),1e0,Q,m,x,1,0e0,R1(1,j),1)
-      if (j < k) then
-c eliminate tail, updating Q(:,j:k) and R1(j:k,j+1:n+1)
-        jj = min(j,n)+1
-        call sqrqhv(m,n+1-j,k-j+1,Q(1,j),m,R1(j,jj),m,x,w)
-c assemble inserted column        
-        R1(j,j) = w
-        do i = j+1,k
-          R1(i,j) = 0e0
-        end do 
-      end if
-      end 

--- a/libcruft/qrupdate/sqrinr.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,73 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqrinr(m,n,Q,Q1,R,R1,j,x)
-c purpose:      updates a QR factorization after inserting a new
-c               row.      
-c               i.e., given an m-by-m orthogonal matrix Q, an m-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:m+1, this subroutine forms the (m+1)-by-(m+1) matrix 
-c               Q1 and an (m+1)-by-n matrix R1 so that Q1 is again 
-c               orthogonal, R1 upper trapezoidal, and 
-c               Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix R. 
-c n (in)        number of columns of the matrix R
-c Q (in)        the orthogonal matrix Q
-c Q1 (out)      the updated matrix Q1
-c R (in)        the upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the new row in R1
-c x (in)        the row being added
-c
-      integer m,n,j
-      real Q(m,m),Q1(m+1,m+1),R(m,n),R1(m+1,n),x(n)
-      external xerbla,slacpy,scopy,sqhqr
-      integer i
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 2
-      else if (j < 1 .or. j > m+1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('SQRINR',info)
-      end if
-c setup the new matrix Q1
-c permute the columns of Q1 and rows of R1 so that c the new row ends 
-c up being the topmost row.      
-      if (j > 1) then
-        call slacpy('0',j-1,m,Q(1,1),m,Q1(1,2),m+1)
-      end if
-      if (j <= m) then
-        call slacpy('0',m-j+1,m,Q(j,1),m,Q1(j+1,2),m+1)
-      end if
-c zero the rest of Q1      
-      do i = 1,m+1
-        Q1(i,1) = 0e0
-        Q1(j,i) = 0e0
-      end do
-      Q1(j,1) = 1e0
-c setup the new matrix R1
-      call scopy(n,x,1,R1(1,1),m+1)
-      call slacpy('0',m,n,R(1,1),m,R1(2,1),m+1)
-c rotate to form proper QR      
-      call sqhqr(m+1,n,m+1,Q1,m+1,R1,m+1)
-      end 

--- a/libcruft/qrupdate/sqrqhu.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqrqhu(m,n,k,Q,ldq,R,ldr,u,rr)
-c purpose:      given an m-by-k matrix Q, an upper trapezoidal 
-c               k-by-n matrix R, and a k-vector u, 
-c               this subroutine updates the matrices Q -> Q1 and 
-c               R -> R1 so that Q1 = Q*G', R1 = G*R, u1(2:k) = 0 
-c               with G orthogonal, R1 upper Hessenberg, and u1 = G*u.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q and rows of R.
-c Q (io)        on entry, the orthogonal matrix Q.
-c               on exit, the updated matrix Q1.
-c ldq (in)      leading dimension of Q.
-c R (io)        on entry, the upper triangular matrix R.
-c               on exit, the updated upper Hessenberg matrix R1.
-c ldr (in)      leading dimension of R.
-c u (in)        the k-vector u.
-c rr (out)      the first element of Q1'*u on exit.
-c
-c               if Q is orthogonal, so is Q1. It is not strictly
-c               necessary, however.
-      integer m,n,k,ldq,ldr
-      real Q(ldq,*),R(ldr,*),u(*),rr
-      real c
-      real s,w
-      external xerbla,slartg,srot
-      integer i,info
-c quick return if possible.
-      if (k <= 0) return
-c check arguments.      
-      info = 0
-      if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('SQRQHU',info)
-      end if
-      rr = u(k)
-      do i = k-1,1,-1
-        w = rr
-        if (w /= 0e0) then
-          call slartg(u(i),w,c,s,rr)
-c apply rotation to rows of R if necessary        
-          if (i <= n) then
-            call srot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s)
-          end if
-c apply rotation to columns of Q if necessary
-          if (m > 0) then
-            call srot(m,Q(1,i),1,Q(1,i+1),1,c,s)
-          end if
-        else
-c no rotation necessary
-          rr = u(i)
-        end if          
-      end do
-      end 
-

--- a/libcruft/qrupdate/sqrqhv.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,75 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqrqhv(m,n,k,Q,ldq,R,ldr,u,rr)
-c purpose:      given an m-by-k matrix Q, an upper trapezoidal 
-c               k-by-n matrix R, and an m-vector u, this subroutine 
-c               updates the matrices Q -> Q1 and R -> R1 so that 
-c               Q1 = Q*G', R1 = G*R, w1(2:m) = 0 with G orthogonal, 
-c               R1 upper Hessenberg, and w1 = Q1'*u.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q and rows of R. k <= m.
-c Q (io)        on entry, the orthogonal matrix Q.
-c               on exit, the updated matrix Q1.
-c ldq (in)      leading dimension of Q.
-c R (io)        on entry, the upper triangular matrix R.
-c               on exit, the updated upper Hessenberg matrix R1.
-c ldr (in)      leading dimension of R.
-c u (in)        the m-vector u.
-c rr (out)      the first element of Q1'*u on exit.
-c
-c               if Q is orthogonal, so is Q1. It is not strictly
-c               necessary, however.
-      integer m,n,k,ldq,ldr
-      real Q(ldq,*),R(ldr,*),u(*),rr
-      real c
-      real s,w,w1,sdot
-      external xerbla,sdot,slartg,srot
-      integer i,info
-c quick return if possible.
-      if (k <= 0) return
-c check arguments.      
-      info = 0
-      if (k > m) then
-        info = 3
-      else if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('SQRQHV',info)
-      end if
-c form each element of w = Q'*u when necessary.
-      rr = sdot(m,Q(1,k),1,u,1)
-      do i = k-1,1,-1
-        w1 = rr
-        w = sdot(m,Q(1,i),1,u,1)
-        call slartg(w,w1,c,s,rr)
-c apply rotation to rows of R if necessary        
-        if (i <= n) then
-          call srot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s)
-        end if
-c apply rotation to columns of Q
-        call srot(m,Q(1,i),1,Q(1,i+1),1,c,s)
-      end do
-      end 
-

--- a/libcruft/qrupdate/sqrshc.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,97 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine sqrshc(m,n,k,Q,R,i,j)
-c purpose:      updates a QR factorization after circular shift of
-c               columns.      
-c               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
-c               upper trapezoidal matrix R and index j in the range 
-c               1:n+1, this subroutine updates the matrix Q -> Q1 and 
-c               R -> R1 so that Q1 is again orthogonal, R1 upper
-c               trapezoidal, and 
-c               Q1*R1 = A(:,p), where A = Q*R and p is the permutation
-c               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
-c               [1:j-1,shift(j:i,+1),i+1:n] if j > i.
-c               if m == 0, the matrix Q is ignored.
-c               (real version)
-c arguments:
-c m (in)        number of rows of the matrix Q, or 0 if Q is not needed.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R.
-c Q (io)        on entry, the (orthogonal) matrix Q.
-c               on exit, the updated matrix Q1
-c R (io)        on entry, the upper trapezoidal m-by-n matrix R.
-c               on exit, the updated matrix R1.
-c i (in)        the first index determining the range (see above)
-c j (in)        the second index determining the range (see above)
-c
-      integer m,n,k,i,j
-      real Q(m,k),R(k,n)
-      external xerbla,sswap,sqhqr,sqrqhu
-      real w
-      integer l,jj,kk,info
-
-c quick return if possible
-      if (k <= 0 .or. n <= 1) return
-      info = 0
-      if (m /= 0 .and. k > m) then
-        info = 3
-      else if (i < 1 .or. i > n) then
-        info = 6
-      else if (j < 1 .or. j > n) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('SQRSHC',info)
-      end if
-
-      if (i < j) then
-c shift columns
-        do l = i,j-1
-          call sswap(min(k,l+1),R(1,l),1,R(1,l+1),1)
-        end do
-c retriangularize
-        if (i < k) then
-          kk = min(k,j)
-          if (m > 0) then
-            call sqhqr(m,n+1-i,kk+1-i,Q(1,i),m,R(i,i),k)
-          else
-            call sqhqr(0,n+1-i,kk+1-i,Q,1,R(i,i),k)
-          endif
-        end if
-      else if (j < i) then
-c shift columns
-        do l = i,j+1,-1
-          call sswap(min(k,i),R(1,l),1,R(1,l-1),1)
-        end do
-c retriangularize
-        if (j < k) then
-          jj = min(j+1,n)
-          kk = min(k,i)
-          if (m > 0) then
-            call sqrqhu(m,n-j,kk+1-j,Q(1,j),m,R(j,jj),k,R(j,j),w)
-          else
-            call sqrqhu(0,n-j,kk+1-j,Q,1,R(j,jj),k,R(j,j),w)
-          end if
-          R(j,j) = w
-          do jj = j+1,kk
-            R(jj,j) = 0
-          end do
-        end if
-      end if
-      end

--- a/libcruft/qrupdate/zch1dn.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine zch1dn(n,R,u,w,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a hermitian positive definite matrix A, i.e.
-c               A = R'*R, this subroutine downdates R -> R1 so that
-c               R1'*R1 = A - u*u' 
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n,info
-      double complex R(n,n),u(n)
-      double precision w(n)
-      external ztrsv,zlartg,dznrm2
-      double precision rho,dznrm2
-      double complex crho,rr,ui,t
-      integer i,j
-
-c quick return if possible
-      if (n <= 0) return
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0d0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call ztrsv('U','C','N',n,R,n,u,1)
-      rho = dznrm2(n,u,1)
-c check positive definiteness      
-      rho = 1 - rho**2
-      if (rho <= 0d0) then
-        info = 1
-        return
-      end if
-      crho = sqrt(rho)
-c eliminate R' \ u
-      do i = n,1,-1
-        ui = u(i)
-c generate next rotation        
-        call zlartg(crho,ui,w(i),u(i),rr)
-        crho = rr
-      end do
-c apply rotations
-      do i = n,1,-1
-        ui = 0d0
-        do j = i,1,-1
-          t = w(j)*ui + u(j)*R(j,i)
-          R(j,i) = w(j)*R(j,i) - conjg(u(j))*ui
-          ui = t
-        end do
-      end do
-
-      info = 0
-      end

--- a/libcruft/qrupdate/zch1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,56 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine zch1up(n,R,u,w)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a hermitian positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A + u*u' or A - u*u'
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated matrix R1
-c u (io)        the vector determining the rank-1 update
-c               on exit, u is destroyed.
-c w (w)         a real workspace vector of size n
-c 
-c NOTE: the workspace vector is used to store the rotations
-c       so that R does not need to be traversed by rows.
-c
-      integer n
-      double complex R(n,n),u(n)
-      double precision w(n)
-      external zlartg
-      double complex rr,ui,t
-      integer i,j
-      
-      do i = 1,n
-c apply stored rotations, column-wise
-        ui = conjg(u(i))
-        do j = 1,i-1
-          t = w(j)*R(j,i) + u(j)*ui
-          ui = w(j)*ui - conjg(u(j))*R(j,i)
-          R(j,i) = t
-        end do
-c generate next rotation        
-        call zlartg(R(i,i),ui,w(i),u(i),rr)
-        R(i,i) = rr
-      end do
-      end
-

--- a/libcruft/qrupdate/zchdex.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,62 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine zchdex(n,R,R1,j)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the deleted row/column
-      integer n,j,info
-      double complex R(n,n),R1(n-1,n-1)
-      double precision c
-      double complex Qdum,s,rr
-      external xerbla,zlacpy,zqhqr,zlartg
-
-c quick return if possible
-      if (n == 1) return
-
-c check arguments      
-      info = 0
-      if (n <= 0) then
-        info = 1
-      else if (j < 1 .or. j > n) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('ZCHDEX',info)
-      end if
-
-c setup the new matrix R1
-      if (j > 1) then
-        call zlacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1)
-      end if
-      if (j < n) then
-        call zlacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1)
-        call zqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1)
-c eliminate R(n,n)      
-        call zlartg(R1(n-1,n-1),R(n,n),c,s,rr)
-        R1(n-1,n-1) = rr
-      endif
-      end

--- a/libcruft/qrupdate/zchinx.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,109 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-
-      subroutine zchinx(n,R,R1,j,u,info)
-c purpose:      given an upper triangular matrix R that is a Cholesky
-c               factor of a symmetric positive definite matrix A, i.e.
-c               A = R'*R, this subroutine updates R -> R1 so that
-c               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
-c               jj = [1:j-1,j+1:n+1].
-c               (complex version)
-c arguments:
-c n (in)        the order of matrix R
-c R (in)        the original upper trapezoidal matrix R
-c R1 (out)      the updated matrix R1
-c j (in)        the position of the inserted row/column
-c u (in)        the vector (n+1) determining the rank-1 update
-c info (out)    on exit, if info = 1, the 
-c               definiteness.
-
-      integer n,j,info
-      double complex R(n,n),R1(n+1,n+1),u(n+1)
-      double precision rho,dznrm2
-      double complex Qdum,w
-      external xerbla,zcopy,zlacpy,ztrsv,dznrm2,zqrqhu
-      integer jj
-
-c quick return if possible
-      if (n == 0) then
-        if (dble(u(1)) <= 0) then
-          info = 1
-          return
-        else
-          R(1,1) = sqrt(dble(u(1)))
-        end if
-      end if
-
-c check arguments      
-      info = 0
-      if (n < 0) then
-        info = 1
-      else if (j < 1 .or. j > n+1) then
-        info = 4
-      end if
-      if (info /= 0) then
-        call xerbla('ZCHINX',info)
-      end if
-
-c copy shifted vector
-      if (j > 1) then
-        call zcopy(j-1,u,1,R1(1,j),1)
-      end if
-      w = u(j)
-      if (j < n+1) then
-        call zcopy(n-j+1,u(j+1),1,R1(j,j),1)
-      end if
-      
-c check for singularity of R
-      do i = 1,n
-        if (R(i,i) == 0d0) then
-          info = 2
-          return
-        end if
-      end do
-c form R' \ u
-      call ztrsv('U','T','N',n,R,n,R1(1,j),1)
-      rho = dznrm2(n,R1(1,j),1)
-c check positive definiteness      
-      rho = u(j) - rho**2
-      if (rho <= 0d0) then
-        info = 1
-        return
-      end if
-      R1(n+1,n+1) = sqrt(rho)
-
-c setup the new matrix R1
-      do i = 1,n+1
-        R1(n+1,i) = 0d0
-      end do
-      if (j > 1) then
-        call zlacpy('0',j-1,n,R(1,1),n,R1(1,1),n+1)
-      end if
-      if (j <= n) then
-        call zlacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1)
-c retriangularize
-        jj = min(j+1,n)
-        call zqrqhu(0,n+1-j,n-j,Qdum,1,R1(j,jj),n+1,R1(j,j),w)
-        R1(j,j) = w
-        do jj = j+1,n
-          R1(jj,j) = 0d0
-        end do
-      end if
-
-      end

--- a/libcruft/qrupdate/zqhqr.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine zqhqr(m,n,k,Q,ldq,R,ldr)
-c purpose:      given an k-by-n upper Hessenberg matrix R and
-c               an m-by-k matrix Q, this subroutine updates
-c               R -> R1 and Q -> Q1 so that R1 is upper 
-c               trapezoidal, R1 = G*R and Q1 = Q*G', where
-c               G is an unitary matrix, giving Q1*R1 = Q*R.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q
-c n (in)        number of columns of the matrix R
-c k (in)        number of columns of Q and rows of R.
-c Q (io)        on entry, the unitary matrix Q
-c               on exit, the updated matrix Q1
-c ldq (in)      leading dimension of Q
-c R (io)        on entry, the upper triangular matrix R
-c               on exit, the updated upper Hessenberg matrix R1
-c ldr (in)      leading dimension of R
-c
-      integer m,n,k,ldq,ldr
-      double complex Q(ldq,*),R(ldr,*)
-      double precision c
-      double complex s,rr
-      external xerbla,zlartg,zrot
-      integer info,i
-c quick return if possible.
-      if (n <= 0 .or. k <= 1) return
-c check arguments.
-      info = 0
-      if (ldq < 1) then
-        info = 5
-      else if (ldr < 1) then
-        info = 7
-      end if
-      if (info /= 0) then
-        call xerbla('ZQHQR',info)
-      end if
-c triangularize      
-      do i = 1,min(k-1,n)
-        call zlartg(R(i,i),R(i+1,i),c,s,rr)
-        R(i,i) = rr
-        R(i+1,i) = 0d0
-        if (i < n) then
-          call zrot(n-i,R(i,i+1),ldr,R(i+1,i+1),ldr,c,s)
-        end if
-c apply rotation to Q        
-        if (m > 0) then
-          call zrot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s))
-        end if
-      end do
-      end 
-

--- a/libcruft/qrupdate/zqr1up.f	Tue Jan 20 21:15:17 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,53 +0,0 @@
-c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
-c 
-c Author: Jaroslav Hajek <highegg@gmail.com>
-c 
-c This source is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c 
-c This program is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-c GNU General Public License for more details.
-c 
-c You should have received a copy of the GNU General Public License
-c along with this software; see the file COPYING.  If not, see
-c <http://www.gnu.org/licenses/>.
-c 
-      subroutine zqr1up(m,n,k,Q,R,u,v)
-c purpose:      updates a QR factorization after rank-1 modification
-c               i.e., given a m-by-k unitary Q and m-by-n upper 
-c               trapezoidal R, an m-vector u and n-vector v, 
-c               this subroutine updates Q -> Q1 and R -> R1 so that
-c               Q1*R1 = Q*R + Q*Q'u*v', and Q1 is again unitary
-c               and R1 upper trapezoidal.
-c               (complex version)
-c arguments:
-c m (in)        number of rows of the matrix Q.
-c n (in)        number of columns of the matrix R.
-c k (in)        number of columns of Q, and rows of R. k <= m.
-c Q (io)        on entry, the unitary m-by-k matrix Q.
-c               on exit, the updated matrix Q1.
-c R (io)        on entry, the upper trapezoidal m-by-n matrix R.
-c               on exit, the updated matrix R1.
-c u (in)        the left m-vector.
-c v (in)        the right n-vector.
-c
-      integer m,n,k
-      double complex Q(m,k),R(k,n),u(m),v(n)
-      double complex w
-      external zqrqhv,zqhqr
-      integer i
-c quick return if possible      
-      if (m <= 0 .or. n <= 0) return
-c eliminate tail of Q'*u
-      call zqrqhv(m,n,k,Q,m,R,m,u,w)
-c update R      
-      do i = 1,n
-        R(1,i) = R(1,i) + w*conjg(v(i))
-      end do
-c retriangularize R
-      call zqhqr(m,n,k,Q,m,R,k)
-      end