# HG changeset patch
# User Jaroslav Hajek
# Date 1204469199 -3600
# Node ID 320345866761954a71bf7724787c21bdb6023d88
# Parent 205ddddc2990371c55c2a3f092019e902fc177b6
implement Cholesky factorization updating
diff -r 205ddddc2990 -r 320345866761 libcruft/ChangeLog
--- a/libcruft/ChangeLog Sun Feb 24 13:27:33 2008 +0100
+++ b/libcruft/ChangeLog Sun Mar 02 15:46:39 2008 +0100
@@ -1,3 +1,9 @@ 2008-02-24 Jaroslav Hajek
+
+ * qrupdate/dch1dn.f, qrupdate/dch1up.f,
+ qrupdate/zch1dn.f, qrupdate/zch1up.f: add new sources
+ * Makefile.in: include new sources
+
2008-02-24 Jaroslav Hajek
* qrupdate/: add new library
diff -r 205ddddc2990 -r 320345866761 libcruft/qrupdate/Makefile.in
--- a/libcruft/qrupdate/Makefile.in Sun Feb 24 13:27:33 2008 +0100
+++ b/libcruft/qrupdate/Makefile.in Sun Mar 02 15:46:39 2008 +0100
@@ -33,6 +33,8 @@ FSRC = dqr1up.f zqr1up.f \
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
diff -r 205ddddc2990 -r 320345866761 libcruft/qrupdate/dch1dn.f
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/dch1dn.f Sun Mar 02 15:46:39 2008 +0100
@@ -0,0 +1,79 @@
+c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic
+c
+c Author: Jaroslav Hajek
+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 .
+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
diff -r 205ddddc2990 -r 320345866761 libcruft/qrupdate/dch1up.f
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/dch1up.f Sun Mar 02 15:46:39 2008 +0100
@@ -0,0 +1,57 @@
+c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic
+c
+c Author: Jaroslav Hajek
+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 .
+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
+
diff -r 205ddddc2990 -r 320345866761 libcruft/qrupdate/zch1dn.f
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/zch1dn.f Sun Mar 02 15:46:39 2008 +0100
@@ -0,0 +1,79 @@
+c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic
+c
+c Author: Jaroslav Hajek
+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 .
+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
diff -r 205ddddc2990 -r 320345866761 libcruft/qrupdate/zch1up.f
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/qrupdate/zch1up.f Sun Mar 02 15:46:39 2008 +0100
@@ -0,0 +1,56 @@
+c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic
+c
+c Author: Jaroslav Hajek
+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 .
+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
+
diff -r 205ddddc2990 -r 320345866761 liboctave/ChangeLog
--- a/liboctave/ChangeLog Sun Feb 24 13:27:33 2008 +0100
+++ b/liboctave/ChangeLog Sun Mar 02 15:46:39 2008 +0100
@@ -1,3 +1,8 @@ 2008-02-24 Jaroslav Hajek
+
+ * dbleCHOL.h, dbleCHOL.cc, CmplxCHOL.h, CmplxCHOL.cc:
+ Implement rank-1 updating and downdating methods.
+
2008-02-24 Jaroslav Hajek
* dbleQR.h, dbleQR.cc, CmplxQR.h, CmplxQR.cc:
diff -r 205ddddc2990 -r 320345866761 liboctave/CmplxCHOL.cc
--- a/liboctave/CmplxCHOL.cc Sun Feb 24 13:27:33 2008 +0100
+++ b/liboctave/CmplxCHOL.cc Sun Mar 02 15:46:39 2008 +0100
@@ -21,10 +21,13 @@ along with Octave; see the file COPYING.
*/
+// updating/downdating by Jaroslav Hajek 2008
+
#ifdef HAVE_CONFIG_H
#include
#endif
+#include
#include "dMatrix.h"
#include "dRowVector.h"
#include "CmplxCHOL.h"
@@ -47,6 +50,13 @@ extern "C"
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*,
+ int &);
+
}
octave_idx_type
@@ -151,6 +161,47 @@ ComplexCHOL::inverse (void) const
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 (ComplexMatrix u)
+{
+ int n = chol_mat.rows ();
+ if (u.length () != n)
+ {
+ (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+ return;
+ }
+
+ OCTAVE_LOCAL_BUFFER(double, w, n);
+
+ F77_XFCN(zch1up, ZCH1UP, (n, chol_mat.fortran_vec (),
+ u.fortran_vec (), w));
+}
+
+int ComplexCHOL::downdate (ComplexMatrix u)
+{
+ int n = chol_mat.rows ();
+ if (u.length () != n)
+ {
+ (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+ return -1;
+ }
+
+ OCTAVE_LOCAL_BUFFER(double, w, n);
+ int info;
+
+ F77_XFCN(zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (),
+ u.fortran_vec (), w, info));
+
+ return info;
+}
+
ComplexMatrix
chol2inv (const ComplexMatrix& r)
{
diff -r 205ddddc2990 -r 320345866761 liboctave/CmplxCHOL.h
--- a/liboctave/CmplxCHOL.h Sun Feb 24 13:27:33 2008 +0100
+++ b/liboctave/CmplxCHOL.h Sun Mar 02 15:46:39 2008 +0100
@@ -20,6 +20,8 @@ along with Octave; see the file COPYING.
.
*/
+
+// updating/downdating by Jaroslav Hajek 2008
#if !defined (octave_ComplexCHOL_h)
#define octave_ComplexCHOL_h 1
@@ -63,6 +65,12 @@ public:
ComplexMatrix inverse (void) const;
+ void set (const ComplexMatrix& R);
+
+ void update (ComplexMatrix u);
+
+ int downdate (ComplexMatrix u);
+
friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexCHOL& a);
private:
diff -r 205ddddc2990 -r 320345866761 liboctave/CmplxQR.h
--- a/liboctave/CmplxQR.h Sun Feb 24 13:27:33 2008 +0100
+++ b/liboctave/CmplxQR.h Sun Mar 02 15:46:39 2008 +0100
@@ -20,6 +20,8 @@ along with Octave; see the file COPYING.
.
*/
+
+// updating/downdating by Jaroslav Hajek 2008
#if !defined (octave_ComplexQR_h)
#define octave_ComplexQR_h 1
diff -r 205ddddc2990 -r 320345866761 liboctave/dbleCHOL.cc
--- a/liboctave/dbleCHOL.cc Sun Feb 24 13:27:33 2008 +0100
+++ b/liboctave/dbleCHOL.cc Sun Mar 02 15:46:39 2008 +0100
@@ -21,10 +21,13 @@ along with Octave; see the file COPYING.
*/
+// updating/downdating by Jaroslav Hajek 2008
+
#ifdef HAVE_CONFIG_H
#include
#endif
+#include
#include "dRowVector.h"
#include "dbleCHOL.h"
#include "f77-fcn.h"
@@ -47,6 +50,13 @@ extern "C"
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,47 @@ CHOL::inverse (void) const
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 (Matrix u)
+{
+ int n = chol_mat.rows ();
+ if (u.length () != n)
+ {
+ (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+ return;
+ }
+
+ OCTAVE_LOCAL_BUFFER(double, w, n);
+
+ F77_XFCN(dch1up, DCH1UP, (n, chol_mat.fortran_vec (),
+ u.fortran_vec (), w));
+}
+
+int CHOL::downdate (Matrix u)
+{
+ int n = chol_mat.rows ();
+ if (u.length () != n)
+ {
+ (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+ return -1;
+ }
+
+ OCTAVE_LOCAL_BUFFER(double, w, n);
+ int info;
+
+ F77_XFCN(dch1dn, DCH1DN, (n, chol_mat.fortran_vec (),
+ u.fortran_vec (), w, info));
+
+ return info;
+}
+
Matrix
chol2inv (const Matrix& r)
{
diff -r 205ddddc2990 -r 320345866761 liboctave/dbleCHOL.h
--- a/liboctave/dbleCHOL.h Sun Feb 24 13:27:33 2008 +0100
+++ b/liboctave/dbleCHOL.h Sun Mar 02 15:46:39 2008 +0100
@@ -20,6 +20,8 @@ along with Octave; see the file COPYING.
.
*/
+
+// updating/downdating by Jaroslav Hajek 2008
#if !defined (octave_CHOL_h)
#define octave_CHOL_h 1
@@ -60,6 +62,12 @@ public:
// Compute the inverse of a matrix using the Cholesky factorization.
Matrix inverse (void) const;
+ void set(const Matrix& R);
+
+ void update (Matrix u);
+
+ int downdate (Matrix u);
+
friend OCTAVE_API std::ostream& operator << (std::ostream& os, const CHOL& a);
private:
diff -r 205ddddc2990 -r 320345866761 liboctave/dbleQR.cc
--- a/liboctave/dbleQR.cc Sun Feb 24 13:27:33 2008 +0100
+++ b/liboctave/dbleQR.cc Sun Mar 02 15:46:39 2008 +0100
@@ -20,8 +20,6 @@ along with Octave; see the file COPYING.
.
*/
-
-// updating/downdating by Jaroslav Hajek 2008
#ifdef HAVE_CONFIG_H
#include
diff -r 205ddddc2990 -r 320345866761 src/ChangeLog
--- a/src/ChangeLog Sun Feb 24 13:27:33 2008 +0100
+++ b/src/ChangeLog Sun Mar 02 15:46:39 2008 +0100
@@ -1,3 +1,8 @@ 2008-02-24 Jaroslav Hajek
+
+ * DLD-FUNCTIONS/cholupdate.cc: add new source
+ * Makefile.in: include new source in build
+
2008-02-24 Jaroslav Hajek
* DLD-FUNCTIONS/qrupdate.cc, DLD-FUNCTIONS/qrinsert.cc,
diff -r 205ddddc2990 -r 320345866761 src/DLD-FUNCTIONS/cholupdate.cc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/cholupdate.cc Sun Mar 02 15:46:39 2008 +0100
@@ -0,0 +1,182 @@
+/* Copyright (C) 2008 VZLU Prague, a.s., Czech Republic
+ *
+ * Author: Jaroslav Hajek
+ *
+ * This source 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 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program 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 this software; see the file COPYING. If not, see
+ * .
+ */
+
+#ifdef HAVE_CONFIG_H
+#include
+#endif
+
+#include "dbleCHOL.h"
+#include "CmplxCHOL.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "oct-obj.h"
+#include "utils.h"
+
+DEFUN_DLD (cholupdate, args, nargout,
+"-*- texinfo -*-\n\
address@hidden Function}\
+{[R1,info]} = cholupdate(R,u,op)\n\
+Updates or downdates a Cholesky factorization. Given an upper triangular\n\
+matrix @var{R} and a column vector @var{u},\n\
+this function tries to determine another upper triangular\n\
+matrix @var{R1} so that\n\
address@hidden @bullet\n\
address@hidden
address@hidden'address@hidden = @var{R}'address@hidden + @address@hidden'\n\
+if @var{op} = \"+\"\n\
address@hidden
address@hidden'address@hidden = @var{R}'address@hidden - @address@hidden'\n\
+if @var{op} = \"-\"\n\
address@hidden itemize\n\
+\n\
+If @var{op} = \"-\", @var{info} is set to\n\
address@hidden
address@hidden 0 if the downdate was successful,\n\
address@hidden 1 if @var{R}'address@hidden - @address@hidden' is not positive definite,\n\
address@hidden 2 if @var{R} is singular.\n\
address@hidden itemize\n\
+\n\
+If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
address@hidden,qrupdate}\n\
address@hidden deftypefn")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ octave_value argR,argu,argop;
+
+ if ((nargin == 3 || nargin == 2) && nargout <= 2
+ && (argR = args (0),argR.is_matrix_type ())
+ && (argu = args (1),argu.is_matrix_type ())
+ && (nargin < 3 || (argop = args(2), argop.is_string ())))
+ {
+ int 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);
+
+ retval (0) = fact.chol_matrix ();
+
+ if (nargout > 1)
+ retval (1) = err;
+ else if (err)
+ error ("cholupdate: downdate violates positiveness");
+ }
+ else
+ {
+ // real 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);
+
+ retval (0) = fact.chol_matrix ();
+
+ if (nargout > 1)
+ retval (1) = err;
+ else if (err)
+ error ("cholupdate: downdate violates positiveness");
+ }
+
+ }
+ 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)
+*/
diff -r 205ddddc2990 -r 320345866761 src/Makefile.in
--- a/src/Makefile.in Sun Feb 24 13:27:33 2008 +0100
+++ b/src/Makefile.in Sun Mar 02 15:46:39 2008 +0100
@@ -69,7 +69,7 @@ DLD_XSRC := balance.cc besselj.cc betain
gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
givens.cc hess.cc inv.cc kron.cc lsode.cc \
lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \
- qrdelete.cc qrinsert.cc qrupdate.cc \
+ qrdelete.cc qrinsert.cc qrupdate.cc cholupdate.cc \
quad.cc qz.cc rand.cc regexp.cc schur.cc sparse.cc \
spparms.cc sqrtm.cc svd.cc syl.cc symrcm.cc symbfact.cc \
time.cc tsearch.cc typecast.cc \