Mercurial > octave-nkf

changeset 7554:40574114c514

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement Cholesky factorization updating
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 04 Mar 2008 22:25:50 -0500
parents 56be6f31dd4e
children df583cd2f21e
files libcruft/ChangeLog libcruft/qrupdate/Makefile.in libcruft/qrupdate/dch1dn.f libcruft/qrupdate/dch1up.f libcruft/qrupdate/zch1dn.f libcruft/qrupdate/zch1up.f liboctave/ChangeLog liboctave/CmplxCHOL.cc liboctave/CmplxCHOL.h liboctave/CmplxQR.cc liboctave/CmplxQR.h liboctave/dbleCHOL.cc liboctave/dbleCHOL.h liboctave/dbleQR.cc src/ChangeLog src/DLD-FUNCTIONS/chol.cc src/DLD-FUNCTIONS/qr.cc
diffstat 17 files changed, 583 insertions(+), 18 deletions(-) [+]
line wrap: on
line diff
--- a/libcruft/ChangeLog	Tue Mar 04 21:47:11 2008 -0500
+++ b/libcruft/ChangeLog	Tue Mar 04 22:25:50 2008 -0500
@@ -1,4 +1,8 @@
-2008-02-24  Jaroslav Hajek <highegg@gmail.com>
+2008-03-04 Jaroslav Hajek <highegg@gmail.com>
+
+	* qrupdate/dch1dn.f, qrupdate/dch1up.f, 
+	qrupdate/zch1dn.f, qrupdate/zch1up.f: New files.
+	* qrupdate/Makefile.in (FSRC): Add them to the list.
 
 	* qrupdate/Makefile.in, qrupdate/dqhqr.f, qrupdate/dqr1up.f,
 	qrupdate/dqrdec.f, qrupdate/dqrder.f, qrupdate/dqrinc.f,
--- a/libcruft/qrupdate/Makefile.in	Tue Mar 04 21:47:11 2008 -0500
+++ b/libcruft/qrupdate/Makefile.in	Tue Mar 04 22:25:50 2008 -0500
@@ -31,8 +31,10 @@
        dqrdec.f zqrdec.f \
        dqrinr.f zqrinr.f \
        dqrder.f zqrder.f \
+       dch1up.f zch1up.f \
+       dch1dn.f zch1dn.f \
        dqrqhv.f zqrqhv.f \
-       dqhqr.f	zqhqr.f 
+       dqhqr.f zqhqr.f 
 
 include $(TOPDIR)/Makeconf
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/dch1dn.f	Tue Mar 04 22:25:50 2008 -0500
@@ -0,0 +1,79 @@
+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,dcopy,dlartg,dnrm2
+      double precision rho,dnrm2
+      double precision rr,ui,t
+      integer i,j
+
+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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/dch1up.f	Tue Mar 04 22:25:50 2008 -0500
@@ -0,0 +1,57 @@
+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,dlartv
+      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
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/zch1dn.f	Tue Mar 04 22:25:50 2008 -0500
@@ -0,0 +1,79 @@
+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,zcopy,zlartg,dznrm2
+      double precision rho,dznrm2
+      double complex crho,rr,ui,t
+      integer i,j
+
+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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/zch1up.f	Tue Mar 04 22:25:50 2008 -0500
@@ -0,0 +1,56 @@
+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/liboctave/ChangeLog	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/ChangeLog	Tue Mar 04 22:25:50 2008 -0500
@@ -1,5 +1,12 @@
 2008-03-04  Jaroslav Hajek  <highegg@gmail.com>
 
+	* dbleCHOL.cc (CHOL::set, CHOL::update, CHOL::downdate):
+	New functions.
+	* dbleCHOL.h: Provide decls.
+	* CmplxCHOL.cc (ComplexCHOL::set, ComplexCHOL::update,
+	ComplexCHOL::downdate): New functions.
+	* CmplxCHOL.h: Provide decls.
+
 	* dbleQR.cc (QR::update, QR::insert_col, QR::delete_col,
 	QR::insert_row, QR::delete_row): New methods.
 	(QR::QR (const Matrix&, const MAtrix&)): New constructor.
--- a/liboctave/CmplxCHOL.cc	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/CmplxCHOL.cc	Tue Mar 04 22:25:50 2008 -0500
@@ -21,10 +21,14 @@
 
 */
 
+// updating/downdating by Jaroslav Hajek 2008
+
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
 
+#include <vector>
+
 #include "dMatrix.h"
 #include "dRowVector.h"
 #include "CmplxCHOL.h"
@@ -47,6 +51,12 @@
 			     Complex*, const octave_idx_type&, const double&,
 			     double&, Complex*, double*, 
 			     octave_idx_type& F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, Complex*, double*);
+
+  F77_RET_T
+  F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, Complex*, double*, 
+                             octave_idx_type&);
 }
 
 octave_idx_type
@@ -151,6 +161,55 @@
   return chol2inv_internal (chol_mat);
 }
 
+void
+ComplexCHOL::set (const ComplexMatrix& R)
+{
+  if (R.is_square ()) 
+    chol_mat = R;
+  else
+    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
+}
+
+void
+ComplexCHOL::update (const ComplexMatrix& u)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      ComplexMatrix tmp = u;
+
+      OCTAVE_LOCAL_BUFFER (double, w, n);
+
+      F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (),
+				 tmp.fortran_vec (), w));
+    }
+  else
+    (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+}
+
+octave_idx_type
+ComplexCHOL::downdate (const ComplexMatrix& u)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      ComplexMatrix tmp = u;
+
+      OCTAVE_LOCAL_BUFFER (double, w, n);
+
+      F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (),
+				 tmp.fortran_vec (), w, info));
+    }
+  else
+    (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+
+  return info;
+}
+
 ComplexMatrix
 chol2inv (const ComplexMatrix& r)
 {
--- a/liboctave/CmplxCHOL.h	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/CmplxCHOL.h	Tue Mar 04 22:25:50 2008 -0500
@@ -21,6 +21,8 @@
 
 */
 
+// updating/downdating by Jaroslav Hajek 2008
+
 #if !defined (octave_ComplexCHOL_h)
 #define octave_ComplexCHOL_h 1
 
@@ -63,6 +65,12 @@
 
   ComplexMatrix inverse (void) const;
 
+  void set (const ComplexMatrix& R);
+
+  void update (const ComplexMatrix& u);
+
+  octave_idx_type downdate (const ComplexMatrix& u);
+
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexCHOL& a);
 
 private:
--- a/liboctave/CmplxQR.cc	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/CmplxQR.cc	Tue Mar 04 22:25:50 2008 -0500
@@ -275,10 +275,6 @@
 void
 ComplexQR::economize (void)
 {
-  idx_vector c (':'), i (Range (1, r.columns ()));
-  q = ComplexMatrix (q.index (c, i));
-  r = ComplexMatrix (r.index (i, c));
-#if 0
   octave_idx_type r_nc = r.columns ();
 
   if (r.rows () > r_nc)
@@ -286,7 +282,6 @@
       q.resize (q.rows (), r_nc);
       r.resize (r_nc, r_nc);
     }
-#endif
 }
 
 /*
--- a/liboctave/CmplxQR.h	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/CmplxQR.h	Tue Mar 04 22:25:50 2008 -0500
@@ -21,6 +21,8 @@
 
 */
 
+// updating/downdating by Jaroslav Hajek 2008
+
 #if !defined (octave_ComplexQR_h)
 #define octave_ComplexQR_h 1
 
--- a/liboctave/dbleCHOL.cc	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/dbleCHOL.cc	Tue Mar 04 22:25:50 2008 -0500
@@ -21,10 +21,14 @@
 
 */
 
+// updating/downdating by Jaroslav Hajek 2008
+
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
 
+#include <vector>
+
 #include "dRowVector.h"
 #include "dbleCHOL.h"
 #include "f77-fcn.h"
@@ -47,6 +51,12 @@
 			     double*, const octave_idx_type&, const double&,
 			     double&, double*, octave_idx_type*, 
 			     octave_idx_type& F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, double*, double*);
+
+  F77_RET_T
+  F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, 
+                             int &);
 }
 
 octave_idx_type
@@ -155,6 +165,55 @@
   return chol2inv_internal (chol_mat);
 }
 
+void
+CHOL::set (const Matrix& R)
+{
+  if (R.is_square ()) 
+    chol_mat = R;
+  else
+    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
+}
+
+void
+CHOL::update (const Matrix& u)
+{
+  int n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      Matrix tmp = u;
+
+      OCTAVE_LOCAL_BUFFER (double, w, n);
+
+      F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (),
+				 tmp.fortran_vec (), w));
+    }
+  else
+    (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+}
+
+octave_idx_type
+CHOL::downdate (const Matrix& u)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      Matrix tmp = u;
+
+      OCTAVE_LOCAL_BUFFER (double, w, n);
+
+      F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (),
+				 tmp.fortran_vec (), w, info));
+    }
+  else
+    (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+
+  return info;
+}
+
 Matrix
 chol2inv (const Matrix& r)
 {
--- a/liboctave/dbleCHOL.h	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/dbleCHOL.h	Tue Mar 04 22:25:50 2008 -0500
@@ -21,6 +21,8 @@
 
 */
 
+// updating/downdating by Jaroslav Hajek 2008
+
 #if !defined (octave_CHOL_h)
 #define octave_CHOL_h 1
 
@@ -60,6 +62,12 @@
   // Compute the inverse of a matrix using the Cholesky factorization.
   Matrix inverse (void) const;
 
+  void set (const Matrix& R);
+
+  void update (const Matrix& u);
+
+  octave_idx_type downdate (const Matrix& u);
+
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const CHOL& a);
 
 private:
--- a/liboctave/dbleQR.cc	Tue Mar 04 21:47:11 2008 -0500
+++ b/liboctave/dbleQR.cc	Tue Mar 04 22:25:50 2008 -0500
@@ -21,8 +21,6 @@
 
 */
 
-// updating/downdating by Jaroslav Hajek 2008
-
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
@@ -265,10 +263,6 @@
 void
 QR::economize (void)
 {
-  idx_vector c (':'), i (Range (1, r.columns ()));
-  q = Matrix (q.index (c, i));
-  r = Matrix (r.index (i, c));
-#if 0
   octave_idx_type r_nc = r.columns ();
 
   if (r.rows () > r_nc)
@@ -276,7 +270,6 @@
       q.resize (q.rows (), r_nc);
       r.resize (r_nc, r_nc);
     }
-#endif
 }
 
 
--- a/src/ChangeLog	Tue Mar 04 21:47:11 2008 -0500
+++ b/src/ChangeLog	Tue Mar 04 22:25:50 2008 -0500
@@ -1,4 +1,6 @@
-2008-02-24  Jaroslav Hajek <highegg@gmail.com>
+2008-03-04  Jaroslav Hajek <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/chol.cc (Fcholupdate): New function.
 
 	* DLD-FUNCTIONS/qr.cc (Fqrupdate, Fqrinsert, Fqrdelete):
 	New functions.
--- a/src/DLD-FUNCTIONS/chol.cc	Tue Mar 04 21:47:11 2008 -0500
+++ b/src/DLD-FUNCTIONS/chol.cc	Tue Mar 04 22:25:50 2008 -0500
@@ -21,6 +21,10 @@
 
 */
 
+// The cholupdate function was written by Jaroslav Hajek
+// <highegg@gmail.com>, Copyright (C) 2008  VZLU Prague, a.s., Czech
+// Republic.
+
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
@@ -441,6 +445,157 @@
   return retval;
 }
 
+DEFUN_DLD (cholupdate, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{R1}, @var{info}]} = cholupdate (@var{R}, @var{u}, @var{op})\n\
+Update or downdate a Cholesky factorization.  Given an upper triangular\n\
+matrix @var{R} and a column vector @var{u}, attempt to determine another\n\
+upper triangular matrix @var{R1} such that\n\
+@itemize @bullet\n\
+@item\n\
+@var{R1}'*@var{R1} = @var{R}'*@var{R} + @var{u}*@var{u}'\n\
+if @var{op} is \"+\"\n\
+@item\n\
+@var{R1}'*@var{R1} = @var{R}'*@var{R} - @var{u}*@var{u}'\n\
+if @var{op} is \"-\"\n\
+@end itemize\n\
+\n\
+If @var{op} is \"-\", @var{info} is set to\n\
+@itemize\n\
+@item 0 if the downdate was successful,\n\
+@item 1 if @var{R}'*@var{R} - @var{u}*@var{u}' is not positive definite,\n\
+@item 2 if @var{R} is singular.\n\
+@end itemize\n\
+\n\
+If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
+@seealso{chol, qrupdate}\n\
+@end deftypefn")
+{
+  int nargin = args.length ();
+
+  octave_value_list retval;
+
+  octave_value argR,argu,argop;
+
+  if ((nargin == 3 || nargin == 2)
+      && (argR = args(0), argR.is_matrix_type ())
+      && (argu = args(1), argu.is_matrix_type ())
+      && (nargin < 3 || (argop = args(2), argop.is_string ())))
+    {
+      octave_idx_type n = argR.rows ();
+
+      std::string op = (nargin < 3) ? "+" : argop.string_value();
+
+      bool down = false;
+
+      if (nargin < 3 || (op == "+") || (down = op == "-"))
+        if (argR.columns () == n 
+            && argu.rows () == n && argu.columns () == 1)
+          {
+            if (argR.is_real_matrix () && argu.is_real_matrix ())
+              {
+                // real case
+                Matrix R = argR.matrix_value ();
+                Matrix u = argu.matrix_value ();
+
+                CHOL fact;
+                fact.set (R);
+                int err = 0;
+
+                if (down)
+                  err = fact.downdate (u);
+                else
+                  fact.update (u);
+
+                if (nargout > 1)
+                  retval(1) = err;
+                else if (err)
+                  error ("cholupdate: downdate violates positiveness");
+
+                retval(0) = fact.chol_matrix ();
+              }
+            else
+              {
+                // complex case
+                ComplexMatrix R = argR.complex_matrix_value ();
+                ComplexMatrix u = argu.complex_matrix_value ();
+
+                ComplexCHOL fact;
+                fact.set (R);
+                int err = 0;
+
+                if (down)
+                  err = fact.downdate (u);
+                else
+                  fact.update (u);
+
+                if (nargout > 1)
+                  retval(1) = err;
+                else if (err)
+                  error ("cholupdate: downdate violates positiveness");
+
+                retval(0) = fact.chol_matrix ();
+              }
+          }
+        else
+          error ("cholupdate: dimension mismatch");
+      else
+        error ("cholupdate: op must be \"+\" or \"-\"");
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
+
+/*
+%!test
+%! A = [  0.436997  -0.131721   0.124120  -0.061673 ;
+%!       -0.131721   0.738529   0.019851  -0.140295 ;
+%!        0.124120   0.019851   0.354879  -0.059472 ;
+%!       -0.061673  -0.140295  -0.059472   0.600939 ];
+%! 
+%! u = [  0.98950 ;
+%!        0.39844 ;
+%!        0.63484 ;
+%!        0.13351 ];
+%! 
+%! R = chol(A);
+%! 
+%! R1 = cholupdate(R,u);
+%! 
+%! assert(norm(triu(R1)-R1,Inf) == 0)
+%! assert(norm(R1'*R1 - R'*R - u*u',Inf) < 1e1*eps)
+%! 
+%! R1 = cholupdate(R1,u,"-");
+%! 
+%! assert(norm(triu(R1)-R1,Inf) == 0)
+%! assert(norm(R1 - R,Inf) < 1e1*eps)
+%! 
+%!test
+%! A = [  0.5585528 + 0.0000000i  -0.1662088 - 0.0315341i   0.0107873 + 0.0236411i  -0.0276775 - 0.0186073i ;
+%!       -0.1662088 + 0.0315341i   0.6760061 + 0.0000000i   0.0011452 - 0.0475528i   0.0145967 + 0.0247641i ;
+%!        0.0107873 - 0.0236411i   0.0011452 + 0.0475528i   0.6263149 - 0.0000000i  -0.1585837 - 0.0719763i ;
+%!       -0.0276775 + 0.0186073i   0.0145967 - 0.0247641i  -0.1585837 + 0.0719763i   0.6034234 - 0.0000000i ];
+%! 
+%! u = [ 0.54267 + 0.91519i ;
+%!       0.99647 + 0.43141i ;
+%!       0.83760 + 0.68977i ;
+%!       0.39160 + 0.90378i ];
+%! 
+%! R = chol(A);
+%! 
+%! R1 = cholupdate(R,u);
+%! 
+%! assert(norm(triu(R1)-R1,Inf) == 0)
+%! assert(norm(R1'*R1 - R'*R - u*u',Inf) < 1e1*eps)
+%! 
+%! R1 = cholupdate(R1,u,"-");
+%! 
+%! assert(norm(triu(R1)-R1,Inf) == 0)
+%! assert(norm(R1 - R,Inf) < 1e1*eps)
+*/
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/src/DLD-FUNCTIONS/qr.cc	Tue Mar 04 21:47:11 2008 -0500
+++ b/src/DLD-FUNCTIONS/qr.cc	Tue Mar 04 22:25:50 2008 -0500
@@ -433,7 +433,7 @@
 
 DEFUN_DLD (qrupdate, args, ,
   "-*- texinfo -*-\n\
-@deftypefn{Loadable Function} {[@var{Q1}, @var{R1}]} = qrupdate (@var{Q}, @var{R}, @var{u}, @var{v})\n\
+@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}]} = qrupdate (@var{Q}, @var{R}, @var{u}, @var{v})\n\
 Given a QR@tie{}factorization of a real or complex matrix\n\
 @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\
 @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization\n\
@@ -554,7 +554,7 @@
 
 DEFUN_DLD (qrinsert, args, ,
   "-*- texinfo -*-\n\
-@deftypefn{Loadable Function} {[@var{Q1}, @var{R1}]} = qrinsert (@var{Q}, @var{R}, @var{j}, @var{x}, @var{orient})\n\
+@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}]} = qrinsert (@var{Q}, @var{R}, @var{j}, @var{x}, @var{orient})\n\
 Given a QR@tie{}factorization of a real or complex matrix\n\
 @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\
 @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of\n\
@@ -729,7 +729,7 @@
 
 DEFUN_DLD (qrdelete, args, ,
   "-*- texinfo -*-\n\
-@deftypefn{Loadable Function} {[@var{Q1}, @var{R1}]} = qrdelete (@var{Q}, @var{R}, @var{j}, @var{orient})\n\
+@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}]} = qrdelete (@var{Q}, @var{R}, @var{j}, @var{orient})\n\
 Given a QR@tie{}factorization of a real or complex matrix\n\
 @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\
 @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of\n\