# HG changeset patch # User Jaroslav Hajek # Date 1204687550 18000 # Node ID 40574114c514ac020915e11f6bd6c11131e3ed64 # Parent 56be6f31dd4ecbd22c83acf952af9d2731ce5b95 implement Cholesky factorization updating diff -r 56be6f31dd4e -r 40574114c514 libcruft/ChangeLog --- 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 +2008-03-04 Jaroslav Hajek + + * 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, diff -r 56be6f31dd4e -r 40574114c514 libcruft/qrupdate/Makefile.in --- 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 diff -r 56be6f31dd4e -r 40574114c514 libcruft/qrupdate/dch1dn.f --- /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 +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 56be6f31dd4e -r 40574114c514 libcruft/qrupdate/dch1up.f --- /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 +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 56be6f31dd4e -r 40574114c514 libcruft/qrupdate/zch1dn.f --- /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 +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 56be6f31dd4e -r 40574114c514 libcruft/qrupdate/zch1up.f --- /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 +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 56be6f31dd4e -r 40574114c514 liboctave/ChangeLog --- 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 + * 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. diff -r 56be6f31dd4e -r 40574114c514 liboctave/CmplxCHOL.cc --- 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 #endif +#include + #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) { diff -r 56be6f31dd4e -r 40574114c514 liboctave/CmplxCHOL.h --- 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: diff -r 56be6f31dd4e -r 40574114c514 liboctave/CmplxQR.cc --- 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 } /* diff -r 56be6f31dd4e -r 40574114c514 liboctave/CmplxQR.h --- 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 diff -r 56be6f31dd4e -r 40574114c514 liboctave/dbleCHOL.cc --- 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 #endif +#include + #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) { diff -r 56be6f31dd4e -r 40574114c514 liboctave/dbleCHOL.h --- 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: diff -r 56be6f31dd4e -r 40574114c514 liboctave/dbleQR.cc --- 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 #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 } diff -r 56be6f31dd4e -r 40574114c514 src/ChangeLog --- 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 +2008-03-04 Jaroslav Hajek + + * DLD-FUNCTIONS/chol.cc (Fcholupdate): New function. * DLD-FUNCTIONS/qr.cc (Fqrupdate, Fqrinsert, Fqrdelete): New functions. diff -r 56be6f31dd4e -r 40574114c514 src/DLD-FUNCTIONS/chol.cc --- 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 +// , Copyright (C) 2008 VZLU Prague, a.s., Czech +// Republic. + #ifdef HAVE_CONFIG_H #include #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++ *** diff -r 56be6f31dd4e -r 40574114c514 src/DLD-FUNCTIONS/qr.cc --- 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\