# 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 \