Mercurial > octave-nkf

changeset 8562:a6edd5c23cb5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use replacement methods if qrupdate is not available
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 22 Jan 2009 11:10:47 +0100
parents 66165de2cc42
children 3a3421a9f0bb
files ChangeLog configure.in liboctave/ChangeLog liboctave/CmplxCHOL.cc liboctave/CmplxCHOL.h liboctave/CmplxQR.cc liboctave/CmplxQR.h liboctave/dbleCHOL.cc liboctave/dbleCHOL.h liboctave/dbleQR.cc liboctave/dbleQR.h liboctave/fCmplxCHOL.cc liboctave/fCmplxCHOL.h liboctave/fCmplxQR.cc liboctave/fCmplxQR.h liboctave/floatCHOL.cc liboctave/floatCHOL.h liboctave/floatQR.cc liboctave/floatQR.h scripts/ChangeLog scripts/optimization/fsolve.m src/ChangeLog src/DLD-FUNCTIONS/chol.cc src/DLD-FUNCTIONS/qr.cc
diffstat 24 files changed, 1668 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog	Wed Jan 21 21:48:24 2009 -0500
+++ b/ChangeLog	Thu Jan 22 11:10:47 2009 +0100
@@ -1,3 +1,7 @@
+2009-01-22  Jaroslav Hajek  <highegg@gmail.com>
+
+	* configure.in: Fix qrupdate warning message.
+
 2009-01-21  John W. Eaton  <jwe@octave.org>
 
 	* Makeconf.in: Substitute X11_INCFLAGS and X11_LIBS.
--- a/configure.in	Wed Jan 21 21:48:24 2009 -0500
+++ b/configure.in	Thu Jan 22 11:10:47 2009 +0100
@@ -885,7 +885,7 @@
      [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."
+warn_qrupdate="qrupdate not found. The QR & Cholesky updating functions will be slow."
 if test "$with_qrupdate" = yes; then
   with_qrupdate=no
   if $have_fortran_compiler; then 
--- a/liboctave/ChangeLog	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/ChangeLog	Thu Jan 22 11:10:47 2009 +0100
@@ -1,3 +1,16 @@
+2009-01-22  Jaroslav Hajek  <highegg@gmail.com>
+
+	* dbleQR.h: Optionally declare warn_qrupdate_once.
+	* dbleQR.cc: Define it.
+	* (CmplxQR.h, dbleQR.h, fCmplxQR.h, floatQR.h): Declare replacement
+	methods unconditionally.
+	* (CmplxQR.cc, dbleQR.cc, fCmplxQR.cc, floatQR.cc): Define
+	updating replacement methods.
+	* (CmplxCHOL.h, dbleCHOL.h, fCmplxCHOL.h, floatCHOL.h): Declare
+	replacement methods unconditionally.
+	* (CmplxCHOL.cc, dbleCHOL.cc, fCmplxCHOL.cc, floatCHOL.cc): Define
+	updating replacement methods.
+
 2009-01-21  Jaroslav Hajek  <highegg@gmail.com>
 
 	* Range.cc ( operator + (double x, const Range& r),
--- a/liboctave/CmplxCHOL.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/CmplxCHOL.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -34,6 +34,9 @@
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
+#ifndef HAVE_QRUPDATE
+#include "dbleQR.h"
+#endif
 
 extern "C"
 {
@@ -291,6 +294,146 @@
     }
 }
 
+#else
+
+void
+ComplexCHOL::update (const ComplexColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      init (chol_mat.hermitian () * chol_mat 
+            + ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false);
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+}
+
+static bool
+singular (const ComplexMatrix& a)
+{
+  for (octave_idx_type i = 0; i < a.rows (); i++)
+    if (a(i,i) == 0.0) return true;
+  return false;
+}
+
+octave_idx_type
+ComplexCHOL::downdate (const ComplexColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else
+        {
+          info = init (chol_mat.hermitian () * chol_mat 
+                       - ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false);
+          if (info) info = 1;
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  return info;
+}
+
+octave_idx_type
+ComplexCHOL::insert_sym (const ComplexColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+  
+  if (u.length () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  else if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+  else
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else if (u(j).imag () != 0.0)
+        info = 3;
+      else
+        {
+          ComplexMatrix a = chol_mat.hermitian () * chol_mat;
+          ComplexMatrix a1 (n+1, n+1);
+          for (octave_idx_type k = 0; k < n+1; k++)
+            for (octave_idx_type l = 0; l < n+1; l++)
+              {
+                if (l == j)
+                  a1(k, l) = u(k);
+                else if (k == j)
+                  a1(k, l) = std::conj (u(l));
+                else
+                  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
+              }
+          info = init (a1, false);
+          if (info) info = 1;
+        }
+    }
+
+  return info;
+}
+
+void
+ComplexCHOL::delete_sym (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+  
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+  else
+    {
+      ComplexMatrix a = chol_mat.hermitian () * chol_mat;
+      a.delete_elements (1, idx_vector (j));
+      a.delete_elements (0, idx_vector (j));
+      init (a, false);
+    }
+}
+
+void
+ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+  
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+  else
+    {
+      ComplexMatrix a = chol_mat.hermitian () * chol_mat;
+      Array<octave_idx_type> p (n);
+      for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+      if (i < j)
+        {
+          for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+          p(j) = i;
+        }
+      else if (j < i)
+        {
+          p(j) = i;
+          for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+        }
+
+      init (a.index (idx_vector (p), idx_vector (p)), false);
+    }
+}
+
 #endif
 
 ComplexMatrix
--- a/liboctave/CmplxCHOL.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/CmplxCHOL.h	Thu Jan 22 11:10:47 2009 +0100
@@ -67,8 +67,6 @@
 
   void set (const ComplexMatrix& R);
 
-#ifdef HAVE_QRUPDATE
-
   void update (const ComplexColumnVector& u);
 
   octave_idx_type downdate (const ComplexColumnVector& u);
@@ -79,8 +77,6 @@
 
   void shift_sym (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexCHOL& a);
 
 private:
--- a/liboctave/CmplxQR.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/CmplxQR.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -168,14 +168,28 @@
 
 ComplexQR::ComplexQR (const ComplexMatrix& q_arg, const ComplexMatrix& r_arg)
 {
-  if (q_arg.columns () != r_arg.rows ()) 
+  octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
+  octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
+  if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
     {
-      (*current_liboctave_error_handler) ("QR dimensions mismatch");
-      return;
+      q = q_arg;
+      r = r_arg;
     }
+  else
+    (*current_liboctave_error_handler) ("QR dimensions mismatch");
+}
 
-  this->q = q_arg;
-  this->r = r_arg;
+QR::type
+ComplexQR::get_type (void) const
+{
+  QR::type retval;
+  if (!q.is_empty () && q.is_square ())
+    retval = QR::std;
+  else if (q.rows () > q.columns () && r.is_square ())
+    retval = QR::economy;
+  else
+    retval = QR::raw;
+  return retval;
 }
 
 #ifdef HAVE_QRUPDATE
@@ -196,7 +210,7 @@
                                  utmp.fortran_vec (), vtmp.fortran_vec (), w, rw));
     }
   else
-    (*current_liboctave_error_handler) ("QR update dimensions mismatch");
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
 }
 
 void
@@ -431,6 +445,249 @@
     }
 }
 
+#else
+
+// Replacement update methods.
+
+void
+ComplexQR::update (const ComplexColumnVector& u, const ComplexColumnVector& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () == m && v.length () == n)
+    {
+      init(q*r + ComplexMatrix (u) * ComplexMatrix (v).hermitian (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+void
+ComplexQR::update (const ComplexMatrix& u, const ComplexMatrix& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
+    {
+      init(q*r + u * v.hermitian (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+static
+ComplexMatrix insert_col (const ComplexMatrix& a, octave_idx_type i,
+                          const ComplexColumnVector& x)
+{
+  ComplexMatrix retval (a.rows (), a.columns () + 1);
+  retval.assign (idx_vector::colon, idx_vector (0, i),
+                 a.index (idx_vector::colon, idx_vector (0, i)));
+  retval.assign (idx_vector::colon, idx_vector (i), x);
+  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
+                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
+  return retval;
+}
+
+static
+ComplexMatrix insert_row (const ComplexMatrix& a, octave_idx_type i,
+                          const ComplexRowVector& x)
+{
+  ComplexMatrix retval (a.rows () + 1, a.columns ());
+  retval.assign (idx_vector (0, i), idx_vector::colon,
+                 a.index (idx_vector (0, i), idx_vector::colon));
+  retval.assign (idx_vector (i), idx_vector::colon, x);
+  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
+                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
+  return retval;
+}
+
+static
+ComplexMatrix delete_col (const ComplexMatrix& a, octave_idx_type i)
+{
+  ComplexMatrix retval = a;
+  retval.delete_elements (1, idx_vector (i));
+  return retval;
+}
+
+static
+ComplexMatrix delete_row (const ComplexMatrix& a, octave_idx_type i)
+{
+  ComplexMatrix retval = a;
+  retval.delete_elements (0, idx_vector (i));
+  return retval;
+}
+
+static
+ComplexMatrix shift_cols (const ComplexMatrix& a, 
+                          octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = a.columns ();
+  Array<octave_idx_type> p (n);
+  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+  if (i < j)
+    {
+      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+      p(j) = i;
+    }
+  else if (j < i)
+    {
+      p(j) = i;
+      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+    }
+
+  return a.index (idx_vector::colon, idx_vector (p));
+}
+
+void
+ComplexQR::insert_col (const ComplexColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > n) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_col (q*r, j, u), get_type ());
+    }
+}
+
+void
+ComplexQR::insert_col (const ComplexMatrix& u, const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, ASCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (u.length () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      ComplexMatrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::insert_col (a, js(i), u.column (i));
+      init (a, get_type ());
+    }
+}
+
+void
+ComplexQR::delete_col (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_col (q*r, j), get_type ());
+    }
+}
+
+void
+ComplexQR::delete_col (const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, DESCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      ComplexMatrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::delete_col (a, js(i));
+      init (a, get_type ());
+    }
+}
+
+void
+ComplexQR::insert_row (const ComplexRowVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square () || u.length () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > m) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_row (q*r, j, u), get_type ());
+    }
+}
+
+void
+ComplexQR::delete_row (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  else if (j < 0 || j > m-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_row (q*r, j), get_type ());
+    }
+}
+
+void
+ComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+  else
+    {
+      init (::shift_cols (q*r, i, j), get_type ());
+    }
+}
+
 #endif
 
 /*
--- a/liboctave/CmplxQR.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/CmplxQR.h	Thu Jan 22 11:10:47 2009 +0100
@@ -65,7 +65,7 @@
 
   ComplexMatrix R (void) const { return r; }
 
-#ifdef HAVE_QRUPDATE
+  QR::type get_type (void) const;
 
   void update (const ComplexColumnVector& u, const ComplexColumnVector& v);
 
@@ -85,8 +85,6 @@
 
   void shift_cols (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend std::ostream&  operator << (std::ostream&, const ComplexQR&);
 
 protected:
--- a/liboctave/dbleCHOL.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/dbleCHOL.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -33,6 +33,9 @@
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
+#ifndef HAVE_QRUPDATE
+#include "dbleQR.h"
+#endif
 
 extern "C"
 {
@@ -294,6 +297,144 @@
     }
 }
 
+#else
+
+void
+CHOL::update (const ColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      init (chol_mat.transpose () * chol_mat 
+            + Matrix (u) * Matrix (u).transpose (), false);
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+}
+
+static bool
+singular (const Matrix& a)
+{
+  for (octave_idx_type i = 0; i < a.rows (); i++)
+    if (a(i,i) == 0.0) return true;
+  return false;
+}
+
+octave_idx_type
+CHOL::downdate (const ColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else
+        {
+          info = init (chol_mat.transpose () * chol_mat 
+                - Matrix (u) * Matrix (u).transpose (), false);
+          if (info) info = 1;
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  return info;
+}
+
+octave_idx_type
+CHOL::insert_sym (const ColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  else if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+  else
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else
+        {
+          Matrix a = chol_mat.transpose () * chol_mat;
+          Matrix a1 (n+1, n+1);
+          for (octave_idx_type k = 0; k < n+1; k++)
+            for (octave_idx_type l = 0; l < n+1; l++)
+              {
+                if (l == j)
+                  a1(k, l) = u(k);
+                else if (k == j)
+                  a1(k, l) = u(l);
+                else
+                  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
+              }
+          info = init (a1, false);
+          if (info) info = 1;
+        }
+    }
+
+  return info;
+}
+
+void
+CHOL::delete_sym (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+  else
+    {
+      Matrix a = chol_mat.transpose () * chol_mat;
+      a.delete_elements (1, idx_vector (j));
+      a.delete_elements (0, idx_vector (j));
+      init (a, false);
+    }
+}
+
+void
+CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+  else
+    {
+      Matrix a = chol_mat.transpose () * chol_mat;
+      Array<octave_idx_type> p (n);
+      for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+      if (i < j)
+        {
+          for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+          p(j) = i;
+        }
+      else if (j < i)
+        {
+          p(j) = i;
+          for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+        }
+
+      init (a.index (idx_vector (p), idx_vector (p)), false);
+    }
+}
+
 #endif
 
 Matrix
--- a/liboctave/dbleCHOL.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/dbleCHOL.h	Thu Jan 22 11:10:47 2009 +0100
@@ -64,8 +64,6 @@
 
   void set (const Matrix& R);
 
-#ifdef HAVE_QRUPDATE
-
   void update (const ColumnVector& u);
 
   octave_idx_type downdate (const ColumnVector& u);
@@ -76,8 +74,6 @@
 
   void shift_sym (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const CHOL& a);
 
 private:
--- a/liboctave/dbleQR.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/dbleQR.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -159,14 +159,28 @@
 
 QR::QR (const Matrix& q_arg, const Matrix& r_arg)
 {
-  if (q_arg.columns () != r_arg.rows ()) 
+  octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
+  octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
+  if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
     {
-      (*current_liboctave_error_handler) ("QR dimensions mismatch");
-      return;
+      q = q_arg;
+      r = r_arg;
     }
+  else
+    (*current_liboctave_error_handler) ("QR dimensions mismatch");
+}
 
-  this->q = q_arg;
-  this->r = r_arg;
+QR::type
+QR::get_type (void) const
+{
+  QR::type retval;
+  if (!q.is_empty () && q.is_square ())
+    retval = QR::std;
+  else if (q.rows () > q.columns () && r.is_square ())
+    retval = QR::economy;
+  else
+    retval = QR::raw;
+  return retval;
 }
 
 #ifdef HAVE_QRUPDATE
@@ -186,7 +200,7 @@
                                  utmp.fortran_vec (), vtmp.fortran_vec (), w));
     }
   else
-    (*current_liboctave_error_handler) ("QR update dimensions mismatch");
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
 }
 
 void
@@ -418,6 +432,261 @@
     }
 }
 
+#else
+
+// Replacement update methods.
+
+void
+QR::update (const ColumnVector& u, const ColumnVector& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () == m && v.length () == n)
+    {
+      init(q*r + Matrix (u) * Matrix (v).transpose (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+void
+QR::update (const Matrix& u, const Matrix& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
+    {
+      init(q*r + u * v.transpose (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+static
+Matrix insert_col (const Matrix& a, octave_idx_type i,
+                        const ColumnVector& x)
+{
+  Matrix retval (a.rows (), a.columns () + 1);
+  retval.assign (idx_vector::colon, idx_vector (0, i),
+                 a.index (idx_vector::colon, idx_vector (0, i)));
+  retval.assign (idx_vector::colon, idx_vector (i), x);
+  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
+                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
+  return retval;
+}
+
+static
+Matrix insert_row (const Matrix& a, octave_idx_type i,
+                        const RowVector& x)
+{
+  Matrix retval (a.rows () + 1, a.columns ());
+  retval.assign (idx_vector (0, i), idx_vector::colon,
+                 a.index (idx_vector (0, i), idx_vector::colon));
+  retval.assign (idx_vector (i), idx_vector::colon, x);
+  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
+                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
+  return retval;
+}
+
+static
+Matrix delete_col (const Matrix& a, octave_idx_type i)
+{
+  Matrix retval = a;
+  retval.delete_elements (1, idx_vector (i));
+  return retval;
+}
+
+static
+Matrix delete_row (const Matrix& a, octave_idx_type i)
+{
+  Matrix retval = a;
+  retval.delete_elements (0, idx_vector (i));
+  return retval;
+}
+
+static
+Matrix shift_cols (const Matrix& a, 
+                        octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = a.columns ();
+  Array<octave_idx_type> p (n);
+  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+  if (i < j)
+    {
+      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+      p(j) = i;
+    }
+  else if (j < i)
+    {
+      p(j) = i;
+      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+    }
+
+  return a.index (idx_vector::colon, idx_vector (p));
+}
+
+void
+QR::insert_col (const ColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > n) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_col (q*r, j, u), get_type ());
+    }
+}
+
+void
+QR::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, ASCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (u.length () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      Matrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::insert_col (a, js(i), u.column (i));
+      init (a, get_type ());
+    }
+}
+
+void
+QR::delete_col (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_col (q*r, j), get_type ());
+    }
+}
+
+void
+QR::delete_col (const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, DESCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      Matrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::delete_col (a, js(i));
+      init (a, get_type ());
+    }
+}
+
+void
+QR::insert_row (const RowVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square () || u.length () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > m) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_row (q*r, j, u), get_type ());
+    }
+}
+
+void
+QR::delete_row (octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  else if (j < 0 || j > m-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_row (q*r, j), get_type ());
+    }
+}
+
+void
+QR::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+  else
+    {
+      init (::shift_cols (q*r, i, j), get_type ());
+    }
+}
+
+void warn_qrupdate_once (void)
+{
+  static bool warned = false;
+  if (! warned)
+    {
+      (*current_liboctave_warning_handler)
+        ("In this version of Octave, QR & Cholesky updating routines\n"
+         "simply update the matrix and recalculate factorizations.\n"
+         "To use fast algorithms, link Octave with the qrupdate library.\n"
+         "See <http://sourceforge.net/projects/qrupdate>.\n");
+      warned = true;
+    }
+}
+
 #endif
 
 /*
--- a/liboctave/dbleQR.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/dbleQR.h	Thu Jan 22 11:10:47 2009 +0100
@@ -70,7 +70,7 @@
 
   Matrix R (void) const { return r; }
 
-#ifdef HAVE_QRUPDATE
+  QR::type get_type (void) const;
 
   void update (const ColumnVector& u, const ColumnVector& v);
 
@@ -90,8 +90,6 @@
 
   void shift_cols (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend std::ostream&  operator << (std::ostream&, const QR&);
 
 protected:
@@ -100,6 +98,10 @@
   Matrix r;
 };
 
+#ifndef HAVE_QRUPDATE
+void warn_qrupdate_once (void);
+#endif
+
 #endif
 
 /*
--- a/liboctave/fCmplxCHOL.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/fCmplxCHOL.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -34,6 +34,9 @@
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
+#ifndef HAVE_QRUPDATE
+#include "dbleQR.h"
+#endif
 
 extern "C"
 {
@@ -291,6 +294,146 @@
     }
 }
 
+#else
+
+void
+FloatComplexCHOL::update (const FloatComplexColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      init (chol_mat.hermitian () * chol_mat 
+            + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false);
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+}
+
+static bool
+singular (const FloatComplexMatrix& a)
+{
+  for (octave_idx_type i = 0; i < a.rows (); i++)
+    if (a(i,i) == 0.0f) return true;
+  return false;
+}
+
+octave_idx_type
+FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else
+        {
+          info = init (chol_mat.hermitian () * chol_mat 
+                       - FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false);
+          if (info) info = 1;
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  return info;
+}
+
+octave_idx_type
+FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  else if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+  else
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else if (u(j).imag () != 0.0)
+        info = 3;
+      else
+        {
+          FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
+          FloatComplexMatrix a1 (n+1, n+1);
+          for (octave_idx_type k = 0; k < n+1; k++)
+            for (octave_idx_type l = 0; l < n+1; l++)
+              {
+                if (l == j)
+                  a1(k, l) = u(k);
+                else if (k == j)
+                  a1(k, l) = std::conj (u(l));
+                else
+                  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
+              }
+          info = init (a1, false);
+          if (info) info = 1;
+        }
+    }
+
+  return info;
+}
+
+void
+FloatComplexCHOL::delete_sym (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+  else
+    {
+      FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
+      a.delete_elements (1, idx_vector (j));
+      a.delete_elements (0, idx_vector (j));
+      init (a, false);
+    }
+}
+
+void
+FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+  else
+    {
+      FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
+      Array<octave_idx_type> p (n);
+      for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+      if (i < j)
+        {
+          for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+          p(j) = i;
+        }
+      else if (j < i)
+        {
+          p(j) = i;
+          for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+        }
+
+      init (a.index (idx_vector (p), idx_vector (p)), false);
+    }
+}
+
 #endif
 
 FloatComplexMatrix
--- a/liboctave/fCmplxCHOL.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/fCmplxCHOL.h	Thu Jan 22 11:10:47 2009 +0100
@@ -67,8 +67,6 @@
 
   void set (const FloatComplexMatrix& R);
 
-#ifdef HAVE_QRUPDATE
-
   void update (const FloatComplexColumnVector& u);
 
   octave_idx_type downdate (const FloatComplexColumnVector& u);
@@ -79,8 +77,6 @@
 
   void shift_sym (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatComplexCHOL& a);
 
 private:
--- a/liboctave/fCmplxQR.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/fCmplxQR.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -168,14 +168,28 @@
 
 FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& q_arg, const FloatComplexMatrix& r_arg)
 {
-  if (q_arg.columns () != r_arg.rows ()) 
+  octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
+  octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
+  if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
     {
-      (*current_liboctave_error_handler) ("QR dimensions mismatch");
-      return;
+      q = q_arg;
+      r = r_arg;
     }
+  else
+    (*current_liboctave_error_handler) ("QR dimensions mismatch");
+}
 
-  this->q = q_arg;
-  this->r = r_arg;
+QR::type
+FloatComplexQR::get_type (void) const
+{
+  QR::type retval;
+  if (!q.is_empty () && q.is_square ())
+    retval = QR::std;
+  else if (q.rows () > q.columns () && r.is_square ())
+    retval = QR::economy;
+  else
+    retval = QR::raw;
+  return retval;
 }
 
 #ifdef HAVE_QRUPDATE
@@ -196,7 +210,7 @@
                                  utmp.fortran_vec (), vtmp.fortran_vec (), w, rw));
     }
   else
-    (*current_liboctave_error_handler) ("QR update dimensions mismatch");
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
 }
 
 void
@@ -430,6 +444,249 @@
     }
 }
 
+#else
+
+// Replacement update methods.
+
+void
+FloatComplexQR::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () == m && v.length () == n)
+    {
+      init(q*r + FloatComplexMatrix (u) * FloatComplexMatrix (v).hermitian (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+void
+FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
+    {
+      init(q*r + u * v.hermitian (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+static
+FloatComplexMatrix insert_col (const FloatComplexMatrix& a, octave_idx_type i,
+                               const FloatComplexColumnVector& x)
+{
+  FloatComplexMatrix retval (a.rows (), a.columns () + 1);
+  retval.assign (idx_vector::colon, idx_vector (0, i),
+                 a.index (idx_vector::colon, idx_vector (0, i)));
+  retval.assign (idx_vector::colon, idx_vector (i), x);
+  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
+                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
+  return retval;
+}
+
+static
+FloatComplexMatrix insert_row (const FloatComplexMatrix& a, octave_idx_type i,
+                               const FloatComplexRowVector& x)
+{
+  FloatComplexMatrix retval (a.rows () + 1, a.columns ());
+  retval.assign (idx_vector (0, i), idx_vector::colon,
+                 a.index (idx_vector (0, i), idx_vector::colon));
+  retval.assign (idx_vector (i), idx_vector::colon, x);
+  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
+                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
+  return retval;
+}
+
+static
+FloatComplexMatrix delete_col (const FloatComplexMatrix& a, octave_idx_type i)
+{
+  FloatComplexMatrix retval = a;
+  retval.delete_elements (1, idx_vector (i));
+  return retval;
+}
+
+static
+FloatComplexMatrix delete_row (const FloatComplexMatrix& a, octave_idx_type i)
+{
+  FloatComplexMatrix retval = a;
+  retval.delete_elements (0, idx_vector (i));
+  return retval;
+}
+
+static
+FloatComplexMatrix shift_cols (const FloatComplexMatrix& a, 
+                               octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = a.columns ();
+  Array<octave_idx_type> p (n);
+  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+  if (i < j)
+    {
+      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+      p(j) = i;
+    }
+  else if (j < i)
+    {
+      p(j) = i;
+      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+    }
+
+  return a.index (idx_vector::colon, idx_vector (p));
+}
+
+void
+FloatComplexQR::insert_col (const FloatComplexColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > n) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_col (q*r, j, u), get_type ());
+    }
+}
+
+void
+FloatComplexQR::insert_col (const FloatComplexMatrix& u, const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, ASCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (u.length () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      FloatComplexMatrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::insert_col (a, js(i), u.column (i));
+      init (a, get_type ());
+    }
+}
+
+void
+FloatComplexQR::delete_col (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_col (q*r, j), get_type ());
+    }
+}
+
+void
+FloatComplexQR::delete_col (const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, DESCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      FloatComplexMatrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::delete_col (a, js(i));
+      init (a, get_type ());
+    }
+}
+
+void
+FloatComplexQR::insert_row (const FloatComplexRowVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square () || u.length () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > m) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_row (q*r, j, u), get_type ());
+    }
+}
+
+void
+FloatComplexQR::delete_row (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  else if (j < 0 || j > m-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_row (q*r, j), get_type ());
+    }
+}
+
+void
+FloatComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+  else
+    {
+      init (::shift_cols (q*r, i, j), get_type ());
+    }
+}
+
 #endif
 
 /*
--- a/liboctave/fCmplxQR.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/fCmplxQR.h	Thu Jan 22 11:10:47 2009 +0100
@@ -66,7 +66,7 @@
 
   FloatComplexMatrix R (void) const { return r; }
 
-#ifdef HAVE_QRUPDATE
+  QR::type get_type (void) const;
 
   void update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v);
 
@@ -86,8 +86,6 @@
 
   void shift_cols (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend std::ostream&  operator << (std::ostream&, const FloatComplexQR&);
 
 protected:
--- a/liboctave/floatCHOL.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/floatCHOL.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -33,6 +33,9 @@
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
+#ifndef HAVE_QRUPDATE
+#include "dbleQR.h"
+#endif
 
 extern "C"
 {
@@ -294,6 +297,144 @@
     }
 }
 
+#else
+
+void
+FloatCHOL::update (const FloatColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      init (chol_mat.transpose () * chol_mat 
+            + FloatMatrix (u) * FloatMatrix (u).transpose (), false);
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+}
+
+static bool
+singular (const FloatMatrix& a)
+{
+  for (octave_idx_type i = 0; i < a.rows (); i++)
+    if (a(i,i) == 0.0f) return true;
+  return false;
+}
+
+octave_idx_type
+FloatCHOL::downdate (const FloatColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else
+        {
+          info = init (chol_mat.transpose () * chol_mat 
+                - FloatMatrix (u) * FloatMatrix (u).transpose (), false);
+          if (info) info = 1;
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  return info;
+}
+
+octave_idx_type
+FloatCHOL::insert_sym (const FloatColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  else if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+  else
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else
+        {
+          FloatMatrix a = chol_mat.transpose () * chol_mat;
+          FloatMatrix a1 (n+1, n+1);
+          for (octave_idx_type k = 0; k < n+1; k++)
+            for (octave_idx_type l = 0; l < n+1; l++)
+              {
+                if (l == j)
+                  a1(k, l) = u(k);
+                else if (k == j)
+                  a1(k, l) = u(l);
+                else
+                  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
+              }
+          info = init (a1, false);
+          if (info) info = 1;
+        }
+    }
+
+  return info;
+}
+
+void
+FloatCHOL::delete_sym (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+  else
+    {
+      FloatMatrix a = chol_mat.transpose () * chol_mat;
+      a.delete_elements (1, idx_vector (j));
+      a.delete_elements (0, idx_vector (j));
+      init (a, false);
+    }
+}
+
+void
+FloatCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+  else
+    {
+      FloatMatrix a = chol_mat.transpose () * chol_mat;
+      Array<octave_idx_type> p (n);
+      for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+      if (i < j)
+        {
+          for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+          p(j) = i;
+        }
+      else if (j < i)
+        {
+          p(j) = i;
+          for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+        }
+
+      init (a.index (idx_vector (p), idx_vector (p)), false);
+    }
+}
+
 #endif
 
 FloatMatrix
--- a/liboctave/floatCHOL.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/floatCHOL.h	Thu Jan 22 11:10:47 2009 +0100
@@ -64,8 +64,6 @@
 
   void set (const FloatMatrix& R);
 
-#ifdef HAVE_QRUPDATE
-
   void update (const FloatColumnVector& u);
 
   octave_idx_type downdate (const FloatColumnVector& u);
@@ -76,8 +74,6 @@
 
   void shift_sym (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatCHOL& a);
 
 private:
--- a/liboctave/floatQR.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/floatQR.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -159,14 +159,28 @@
 
 FloatQR::FloatQR (const FloatMatrix& q_arg, const FloatMatrix& r_arg)
 {
-  if (q_arg.columns () != r_arg.rows ()) 
+  octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
+  octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
+  if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
     {
-      (*current_liboctave_error_handler) ("QR dimensions mismatch");
-      return;
+      q = q_arg;
+      r = r_arg;
     }
+  else
+    (*current_liboctave_error_handler) ("QR dimensions mismatch");
+}
 
-  this->q = q_arg;
-  this->r = r_arg;
+QR::type
+FloatQR::get_type (void) const
+{
+  QR::type retval;
+  if (!q.is_empty () && q.is_square ())
+    retval = QR::std;
+  else if (q.rows () > q.columns () && r.is_square ())
+    retval = QR::economy;
+  else
+    retval = QR::raw;
+  return retval;
 }
 
 #ifdef HAVE_QRUPDATE
@@ -186,7 +200,7 @@
                                  utmp.fortran_vec (), vtmp.fortran_vec (), w));
     }
   else
-    (*current_liboctave_error_handler) ("QR update dimensions mismatch");
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
 }
 
 void
@@ -418,6 +432,249 @@
     }
 }
 
+#else
+
+// Replacement update methods.
+
+void
+FloatQR::update (const FloatColumnVector& u, const FloatColumnVector& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () == m && v.length () == n)
+    {
+      init(q*r + FloatMatrix (u) * FloatMatrix (v).transpose (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+void
+FloatQR::update (const FloatMatrix& u, const FloatMatrix& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
+    {
+      init(q*r + u * v.transpose (), get_type ());
+    }
+  else
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+}
+
+static
+FloatMatrix insert_col (const FloatMatrix& a, octave_idx_type i,
+                        const FloatColumnVector& x)
+{
+  FloatMatrix retval (a.rows (), a.columns () + 1);
+  retval.assign (idx_vector::colon, idx_vector (0, i),
+                 a.index (idx_vector::colon, idx_vector (0, i)));
+  retval.assign (idx_vector::colon, idx_vector (i), x);
+  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
+                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
+  return retval;
+}
+
+static
+FloatMatrix insert_row (const FloatMatrix& a, octave_idx_type i,
+                        const FloatRowVector& x)
+{
+  FloatMatrix retval (a.rows () + 1, a.columns ());
+  retval.assign (idx_vector (0, i), idx_vector::colon,
+                 a.index (idx_vector (0, i), idx_vector::colon));
+  retval.assign (idx_vector (i), idx_vector::colon, x);
+  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
+                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
+  return retval;
+}
+
+static
+FloatMatrix delete_col (const FloatMatrix& a, octave_idx_type i)
+{
+  FloatMatrix retval = a;
+  retval.delete_elements (1, idx_vector (i));
+  return retval;
+}
+
+static
+FloatMatrix delete_row (const FloatMatrix& a, octave_idx_type i)
+{
+  FloatMatrix retval = a;
+  retval.delete_elements (0, idx_vector (i));
+  return retval;
+}
+
+static
+FloatMatrix shift_cols (const FloatMatrix& a, 
+                        octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = a.columns ();
+  Array<octave_idx_type> p (n);
+  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+  if (i < j)
+    {
+      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+      p(j) = i;
+    }
+  else if (j < i)
+    {
+      p(j) = i;
+      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+    }
+
+  return a.index (idx_vector::colon, idx_vector (p));
+}
+
+void
+FloatQR::insert_col (const FloatColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.length () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > n) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_col (q*r, j, u), get_type ());
+    }
+}
+
+void
+FloatQR::insert_col (const FloatMatrix& u, const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, ASCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (u.length () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      FloatMatrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::insert_col (a, js(i), u.column (i));
+      init (a, get_type ());
+    }
+}
+
+void
+FloatQR::delete_col (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_col (q*r, j), get_type ());
+    }
+}
+
+void
+FloatQR::delete_col (const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, DESCENDING);
+  octave_idx_type nj = js.length ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else if (nj > 0)
+    {
+      FloatMatrix a = q*r;
+      for (octave_idx_type i = 0; i < js.length (); i++)
+        a = ::delete_col (a, js(i));
+      init (a, get_type ());
+    }
+}
+
+void
+FloatQR::insert_row (const FloatRowVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square () || u.length () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  else if (j < 0 || j > m) 
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+  else
+    {
+      init (::insert_row (q*r, j, u), get_type ());
+    }
+}
+
+void
+FloatQR::delete_row (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  else if (j < 0 || j > m-1) 
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+  else
+    {
+      init (::delete_row (q*r, j), get_type ());
+    }
+}
+
+void
+FloatQR::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+  else
+    {
+      init (::shift_cols (q*r, i, j), get_type ());
+    }
+}
+
 #endif
 
 
--- a/liboctave/floatQR.h	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/floatQR.h	Thu Jan 22 11:10:47 2009 +0100
@@ -64,7 +64,7 @@
 
   FloatMatrix R (void) const { return r; }
 
-#ifdef HAVE_QRUPDATE
+  QR::type get_type (void) const;
 
   void update (const FloatColumnVector& u, const FloatColumnVector& v);
 
@@ -84,8 +84,6 @@
 
   void shift_cols (octave_idx_type i, octave_idx_type j);
 
-#endif
-
   friend std::ostream&  operator << (std::ostream&, const FloatQR&);
 
 protected:
--- a/scripts/ChangeLog	Wed Jan 21 21:48:24 2009 -0500
+++ b/scripts/ChangeLog	Thu Jan 22 11:10:47 2009 +0100
@@ -1,3 +1,7 @@
+2009-01-22  Jaroslav Hajek  <highegg@gmail.com>
+
+	* optimization/fsolve.m: Undo the last change.
+
 2009-01-18  Thorsten Meyer  <thorsten.meyier@gmx.de>
 
 	* miscellaneous/doc.m: Add test for existence of info file.
--- a/scripts/optimization/fsolve.m	Wed Jan 21 21:48:24 2009 -0500
+++ b/scripts/optimization/fsolve.m	Thu Jan 22 11:10:47 2009 +0100
@@ -67,8 +67,6 @@
 
 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options)
 
-  persistent have_qrupdate = exist ('qrupdate') == 5;
-
   if (nargin < 3)
     options = struct ();
   endif
@@ -268,11 +266,7 @@
       endif
 
       ## Update the QR factorization.
-      if (have_qrupdate)
-        [q, r] = qrupdate (q, r, u, v);
-      else
-        [q, r] = qr (q*r + u*v');
-      endif
+      [q, r] = qrupdate (q, r, u, v);
 
     endwhile
   endwhile
--- a/src/ChangeLog	Wed Jan 21 21:48:24 2009 -0500
+++ b/src/ChangeLog	Thu Jan 22 11:10:47 2009 +0100
@@ -1,3 +1,8 @@
+2009-01-22  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/qr.cc: Remove HAVE_QRUPDATE check.
+	* DLD-FUNCTIONS/chol.cc: Dtto.
+
 2009-01-21  John W. Eaton  <jwe@octave.org>
 
 	* Makefile.in (display.o): Add X11_INCFLAGS to CPPFLAGS.
--- a/src/DLD-FUNCTIONS/chol.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/src/DLD-FUNCTIONS/chol.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -576,8 +576,6 @@
   return retval;
 }
 
-#ifdef HAVE_QRUPDATE
-
 DEFUN_DLD (cholupdate, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\
@@ -1276,8 +1274,6 @@
 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
 */
 
-#endif
-
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/src/DLD-FUNCTIONS/qr.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/src/DLD-FUNCTIONS/qr.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -739,8 +739,6 @@
 
 */
 
-#ifdef HAVE_QRUPDATE
-
 static
 bool check_qr_dims (const octave_value& q, const octave_value& r,
                     bool allow_ecf = false)
@@ -1580,8 +1578,6 @@
 %! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
 */
 
-#endif
-
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***