# HG changeset patch # User Jaroslav Hajek # Date 1203856053 -3600 # Node ID 205ddddc2990371c55c2a3f092019e902fc177b6 # Parent 8a6965a011764c14a5d33e15d498453c3613fe3b implementation of QR factorization updating diff -r 8a6965a01176 -r 205ddddc2990 configure.in --- a/configure.in Sun Feb 24 03:32:49 2008 -0500 +++ b/configure.in Sun Feb 24 13:27:33 2008 +0100 @@ -1834,6 +1834,7 @@ AC_CONFIG_FILES([Makefile octMakefile Ma libcruft/lapack/Makefile libcruft/minpack/Makefile libcruft/misc/Makefile libcruft/odepack/Makefile libcruft/ordered-qz/Makefile libcruft/quadpack/Makefile + libcruft/qrupdate/Makefile libcruft/ranlib/Makefile libcruft/slatec-fn/Makefile libcruft/slatec-err/Makefile libcruft/villad/Makefile libcruft/blas-xtra/Makefile libcruft/lapack-xtra/Makefile]) diff -r 8a6965a01176 -r 205ddddc2990 libcruft/ChangeLog --- a/libcruft/ChangeLog Sun Feb 24 03:32:49 2008 -0500 +++ b/libcruft/ChangeLog Sun Feb 24 13:27:33 2008 +0100 @@ -1,3 +1,15 @@ 2008-02-14 John W. Eaton + + * qrupdate/: add new library + * qrupdate/dqhqr.f, qrupdate/dqr1up.f, qrupdate/dqrdec.f, + qrupdate/dqrder.f, qrupdate/dqrinc.f, qrupdate/dqrinr.f, + qrupdate/dqrqhv.f, qrupdate/zqhqr.f, qrupdate/zqr1up.f, + qrupdate/zqrdec.f, qrupdate/zqrder.f, qrupdate/zqrinc.f, + qrupdate/zqrinr.f, qrupdate/zqrqhv.f: add new sources + + * qrupdate/Makefile.in: add configuration file + * Makefile.in: include qrupdate in build process + 2008-02-14 John W. Eaton * misc/f77-fcn.h (F77_XFCN): Call octave_rethrow_exception here diff -r 8a6965a01176 -r 205ddddc2990 libcruft/Makefile.in --- a/libcruft/Makefile.in Sun Feb 24 03:32:49 2008 -0500 +++ b/libcruft/Makefile.in Sun Feb 24 13:27:33 2008 +0100 @@ -43,7 +43,7 @@ INSTALL_DATA = @INSTALL_DATA@ CRUFT_DIRS = amos @BLAS_DIR@ blas-xtra daspk dasrt dassl \ @FFT_DIR@ @LAPACK_DIR@ lapack-xtra minpack \ - misc odepack ordered-qz quadpack ranlib \ + misc odepack ordered-qz quadpack qrupdate ranlib \ slatec-err slatec-fn villad SUBDIRS = $(CRUFT_DIRS) diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/Makefile.in --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/Makefile.in Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,41 @@ +# Makefile for octave's libcruft/factupd directory +# +# Copyright (C) 2008 VZLU, a.s. +# +# This file is part of Octave. +# +# Octave 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 3 of the License, or (at +# your option) any later version. +# +# Octave 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 Octave; see the file COPYING. If not, see +# . +# +# Author: Jaroslav Hajek + +TOPDIR = ../.. + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ + +EXTERNAL_DISTFILES = $(DISTFILES) + +FSRC = dqr1up.f zqr1up.f \ + dqrinc.f zqrinc.f \ + dqrdec.f zqrdec.f \ + dqrinr.f zqrinr.f \ + dqrder.f zqrder.f \ + dqrqhv.f zqrqhv.f \ + dqhqr.f zqhqr.f + +include $(TOPDIR)/Makeconf + +include ../Makerules diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/dqhqr.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqhqr.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,69 @@ +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 dqhqr(m,n,k,Q,ldq,R,ldr) +c purpose: given an k-by-n upper Hessenberg matrix R and +c an m-by-k matrix Q, this subroutine updates +c R -> R1 and Q -> Q1 so that R1 is upper +c trapezoidal, R1 = G*R and Q1 = Q*G', where +c G is an orthogonal matrix, giving Q1*R1 = Q*R. +c (real version) +c arguments: +c m (in) number of rows of the matrix Q +c n (in) number of columns of the matrix R +c k (in) number of columns of Q and rows of R. +c Q (io) on entry, the orthogonal matrix Q +c on exit, the updated matrix Q1 +c ldq (in) leading dimension of Q +c R (io) on entry, the upper triangular matrix R +c on exit, the updated upper Hessenberg matrix R1 +c ldr (in) leading dimension of R +c + integer m,n,k,ldq,ldr + double precision Q(ldq,*),R(ldr,*) + double precision c + double precision s,rr + external dlartg,drot + integer info,i +c quick return if possible. + if (n <= 0 .or. k <= 1) return +c check arguments. + info = 0 + if (ldq < 1) then + info = 5 + else if (ldr < 1) then + info = 7 + end if + if (info /= 0) then + call xerbla('DQHQR',info) + end if +c triangularize + do i = 1,min(k-1,n) + call dlartg(R(i,i),R(i+1,i),c,s,rr) + R(i,i) = rr + R(i+1,i) = 0d0 + if (i < n) then + call drot(n-i,R(i,i+1),ldr,R(i+1,i+1),ldr,c,s) + end if +c apply rotation to Q + if (m > 0) then + call drot(m,Q(1,i),1,Q(1,i+1),1,c,s) + end if + end do + end + diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/dqr1up.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqr1up.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,52 @@ +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 dqr1up(m,n,k,Q,R,u,v) +c purpose: updates a QR factorization after rank-1 modification +c i.e., given a m-by-k orthogonal Q and m-by-n upper +c trapezoidal R, an m-vector u and n-vector v, +c this subroutine updates Q -> Q1 and R -> R1 so that +c Q1*R1 = Q*R + Q*Q'u*v', and Q1 is again orthonormal +c and R1 upper trapezoidal. +c (real version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q, and rows of R. k <= m. +c Q (io) on entry, the orthogonal m-by-k matrix Q. +c on exit, the updated matrix Q1. +c R (io) on entry, the upper trapezoidal m-by-n matrix R.. +c on exit, the updated matrix R1. +c u (in) the left m-vector. +c v (in) the right n-vector. +c + integer m,n,k + double precision Q(m,k),R(k,n),u(m),v(n) + double precision w + external dqrqhv,dqhqr,daxpy +c quick return if possible + if (m <= 0 .or. n <= 0) return +c eliminate tail of Q'*u + call dqrqhv(m,n,k,Q,m,R,m,u,w) +c update R + + call daxpy(n,w,v,1,R(1,1),m) + +c retriangularize R + call dqhqr(m,n,k,Q,m,R,k) + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/dqrdec.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqrdec.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,66 @@ +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 dqrdec(m,n,k,Q,R,R1,j) +c purpose: updates a QR factorization after deleting +c a column. +c i.e., given an m-by-k orthogonal matrix Q, an k-by-n +c upper trapezoidal matrix R and index j in the range +c 1:n+1, this subroutine updates the matrix Q -> Q1 and +c forms an m-by-(n-1) matrix R1 so that Q1 remains +c orthogonal, R1 is upper trapezoidal, and +c Q1*R1 = [A(:,1:j-1) A(:,j+1:n)], where A = Q*R. +c (real version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q, and rows of R. +c Q (io) on entry, the orthogonal m-by-k matrix Q. +c on exit, the updated matrix Q1. +c R (in) the original upper trapezoidal matrix R. +c R1 (out) the updated matrix R1. +c j (in) the position of the deleted column in R. +c 1 <= j <= n. +c + integer m,n,k,j + double precision Q(m,k),R(k,n),R1(k,n-1) + external dcopy,dqhqr + integer info +c quick return if possible + if (m <= 0 .or. k <= 0 .or. n == 1) return +c check arguments + info = 0 + if (n < 1) then + info = 2 + else if (j < 1 .or. j > n) then + info = 7 + end if + if (info /= 0) then + call xerbla('DQRDEC',info) + end if +c copy leading portion + call dcopy(k*(j-1),R,1,R1,1) + if (j < n) then +c copy trailing portion of R + call dcopy(k*(n-j),R(1,j+1),1,R1(1,j),1) +c if necessary, retriangularize R1(j:k,j:n-1) and update Q(:,j:k) + if (j < k) then + call dqhqr(m,n-j,k-j+1,Q(1,j),m,R1(j,j),k) + end if + end if + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/dqrder.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqrder.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,93 @@ +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 dqrder(m,n,Q,Q1,R,R1,j) +c purpose: updates a QR factorization after deleting a row. +c i.e., given an m-by-m orthogonal matrix Q, an m-by-n +c upper trapezoidal matrix R and index j in the range +c 1:m, this subroutine forms the (m-1)-by-(m-1) matrix +c Q1 and an (m-1)-by-n matrix R1 so that Q1 is again +c orthogonal, R1 upper trapezoidal, and +c Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R. +c (real version) +c +c arguments: +c m (in) number of rows of the matrix R. +c n (in) number of columns of the matrix R +c Q (in) the orthogonal matrix Q +c Q1 (out) the updated matrix Q1 +c R (in) the upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the new row in R1 +c + integer m,n,j + double precision Q(m,m),Q1(m-1,m-1),R(m,n),R1(m-1,n) + double precision c + double precision s,rr,w + external xerbla,dlacpy,dlartg,drot,dscal,daxpy + integer i +c quick return if possible + if (m == 1) return +c check arguments + info = 0 + if (m < 1) then + info = 1 + else if (j < 1 .or. j > n) then + info = 7 + end if + if (info /= 0) then + call xerbla('DQRDER',info) + end if +c setup the new matrix Q1 +c permute the columns of Q and rows of R so that the deleted row ends +c up being the topmost row. + if (j > 1) then + call dlacpy('0',j-1,m-1,Q(1,2),m,Q1(1,1),m-1) + end if + if (j < m) then + call dlacpy('0',m-j,m-1,Q(j+1,2),m,Q1(j,1),m-1) + end if +c setup the new matrix R1 + call dlacpy('0',m-1,n,R(2,1),m,R1(1,1),m-1) +c eliminate Q(j,2:m) + w = Q(j,m) + do i = m-1,2,-1 + call dlartg(Q(j,i),w,c,s,rr) + w = rr +c apply rotation to rows of R1 + if (i <= n) then + call drot(n-i+1,R1(i-1,i),m-1,R1(i,i),m-1,c,s) + end if +c apply rotation to columns of Q1 + call drot(m-1,Q1(1,i-1),1,Q1(1,i),1,c,s) + end do +c the last iteration is special, as we don't have the first row of +c R and first column of Q + call dlartg(Q(j,1),w,c,s,rr) + w = rr + call dscal(n,c,R1(1,1),m-1) + call daxpy(n,-s,R(1,1),m,R1(1,1),m-1) +c apply rotation to columns of Q1 + call dscal(m-1,c,Q1(1,1),1) + if (j > 1) then + call daxpy(j-1,-s,Q(1,1),1,Q1(1,1),1) + end if + if (j < m) then + call daxpy(m-j,-s,Q(j+1,1),1,Q1(j,1),1) + end if + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/dqrinc.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqrinc.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,73 @@ +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 dqrinc(m,n,k,Q,R,R1,j,x) +c purpose: updates a QR factorization after inserting a new +c column. +c i.e., given an m-by-k orthogonal matrix Q, an m-by-n +c upper trapezoidal matrix R and index j in the range +c 1:n+1, this subroutine updates the matrix Q -> Q1 and +c forms an m-by-(n+1) matrix R1 so that Q1 is again unitary, +c R1 upper trapezoidal, and +c Q1*R1 = [A(:,1:j-1); Q*Q'*x; A(:,j:n-1)], where A = Q*R. +c (real version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q, and rows of R. k <= m. +c Q (io) on entry, the orthogonal matrix Q. +c on exit, the updated matrix Q1 +c R (in) the original upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the new column in R1 +c x (in) the column being inserted +c + integer m,n,k,j + double precision Q(m,k),R(k,n),R1(k,n+1),x(m) + double precision w + external xerbla,dcopy,dqrqhv + integer info,i,jj +c quick return if possible + if (m <= 0) return +c check arguments + info = 0 + if (n < 0) then + info = 2 + else if (j < 1 .or. j > n+1) then + info = 6 + end if + if (info /= 0) then + call xerbla('DQRDER',info) + end if +c copy leading portion of R + call dcopy(k*(j-1),R,1,R1,1) + if (j <= n) then + call dcopy(k*(n+1-j),R(1,j),1,R1(1,j+1),1) + end if + call dgemv('T',m,min(k,j-1),1d0,Q,m,x,1,0d0,R1(1,j),1) + if (j < k) then +c eliminate tail, updating Q(:,j:k) and R1(j:k,j+1:n+1) + jj = min(j,n)+1 + call dqrqhv(m,n+1-j,k-j+1,Q(1,j),m,R1(j,jj),m,x,w) +c assemble inserted column + R1(j,j) = w + do i = j+1,k + R1(i,j) = 0d0 + end do + end if + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/dqrinr.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqrinr.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,73 @@ +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 dqrinr(m,n,Q,Q1,R,R1,j,x) +c purpose: updates a QR factorization after inserting a new +c row. +c i.e., given an m-by-m orthogonal matrix Q, an m-by-n +c upper trapezoidal matrix R and index j in the range +c 1:m+1, this subroutine forms the (m+1)-by-(m+1) matrix +c Q1 and an (m+1)-by-n matrix R1 so that Q1 is again +c orthogonal, R1 upper trapezoidal, and +c Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R. +c (real version) +c arguments: +c m (in) number of rows of the matrix R. +c n (in) number of columns of the matrix R +c Q (in) the orthogonal matrix Q +c Q1 (out) the updated matrix Q1 +c R (in) the upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the new row in R1 +c x (in) the row being added +c + integer m,n,j + double precision Q(m,m),Q1(m+1,m+1),R(m,n),R1(m+1,n),x(n) + external xerbla,dlacpy,dcopy,dqhqr + integer i +c check arguments + info = 0 + if (n < 0) then + info = 2 + else if (j < 1 .or. j > m+1) then + info = 7 + end if + if (info /= 0) then + call xerbla('DQRINR',info) + end if +c setup the new matrix Q1 +c permute the columns of Q1 and rows of R1 so that c the new row ends +c up being the topmost row. + if (j > 1) then + call dlacpy('0',j-1,m,Q(1,1),m,Q1(1,2),m+1) + end if + if (j <= m) then + call dlacpy('0',m-j+1,m,Q(j,1),m,Q1(j+1,2),m+1) + end if +c zero the rest of Q1 + do i = 1,m+1 + Q1(i,1) = 0d0 + Q1(j,i) = 0d0 + end do + Q1(j,1) = 1d0 +c setup the new matrix R1 + call dcopy(n,x,1,R1(1,1),m+1) + call dlacpy('0',m,n,R(1,1),m,R1(2,1),m+1) +c rotate to form proper QR + call dqhqr(m+1,n,m+1,Q1,m+1,R1,m+1) + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/dqrqhv.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqrqhv.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,75 @@ +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 dqrqhv(m,n,k,Q,ldq,R,ldr,u,rr) +c purpose: given an m-by-k matrix Q, an upper trapezoidal +c k-by-n matrix R, and an m-vector u, this subroutine +c updates the matrices Q -> Q1 and R -> R1 so that +c Q1 = Q*G', R1 = G*R, w1(2:m) = 0 with G orthogonal, +c R1 upper Hessenberg, and w1 = Q1'*u. +c (real version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q and rows of R. k <= m. +c Q (io) on entry, the orthogonal matrix Q. +c on exit, the updated matrix Q1. +c ldq (in) leading dimension of Q. +c R (io) on entry, the upper triangular matrix R. +c on exit, the updated upper Hessenberg matrix R1. +c ldr (in) leading dimension of R. +c u (in) the m-vector u. +c rr (out) the first element of Q1'*u on exit. +c +c if Q is orthogonal, so is Q1. It is not strictly +c necessary, however. + integer m,n,k,ldq,ldr + double precision Q(ldq,*),R(ldr,*),u(*),rr + double precision c + double precision s,w,w1,ddot + external ddot,dlartg,drot + integer i,info +c quick return if possible. + if (k <= 0) return +c check arguments. + info = 0 + if (k > m) then + info = 3 + else if (ldq < 1) then + info = 5 + else if (ldr < 1) then + info = 7 + end if + if (info /= 0) then + call xerbla('DQRQHV',info) + end if +c form each element of w = Q'*u when necessary. + rr = ddot(m,Q(1,k),1,u,1) + do i = k-1,1,-1 + w1 = rr + w = ddot(m,Q(1,i),1,u,1) + call dlartg(w,w1,c,s,rr) +c apply rotation to rows of R if necessary + if (i <= n) then + call drot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s) + end if +c apply rotation to columns of Q + call drot(m,Q(1,i),1,Q(1,i+1),1,c,s) + end do + end + diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/zqhqr.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqhqr.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,69 @@ +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 zqhqr(m,n,k,Q,ldq,R,ldr) +c purpose: given an k-by-n upper Hessenberg matrix R and +c an m-by-k matrix Q, this subroutine updates +c R -> R1 and Q -> Q1 so that R1 is upper +c trapezoidal, R1 = G*R and Q1 = Q*G', where +c G is an unitary matrix, giving Q1*R1 = Q*R. +c (complex version) +c arguments: +c m (in) number of rows of the matrix Q +c n (in) number of columns of the matrix R +c k (in) number of columns of Q and rows of R. +c Q (io) on entry, the unitary matrix Q +c on exit, the updated matrix Q1 +c ldq (in) leading dimension of Q +c R (io) on entry, the upper triangular matrix R +c on exit, the updated upper Hessenberg matrix R1 +c ldr (in) leading dimension of R +c + integer m,n,k,ldq,ldr + double complex Q(ldq,*),R(ldr,*) + double precision c + double complex s,rr + external zlartg,zrot + integer info,i +c quick return if possible. + if (n <= 0 .or. k <= 1) return +c check arguments. + info = 0 + if (ldq < 1) then + info = 5 + else if (ldr < 1) then + info = 7 + end if + if (info /= 0) then + call xerbla('ZQHQR',info) + end if +c triangularize + do i = 1,min(k-1,n) + call zlartg(R(i,i),R(i+1,i),c,s,rr) + R(i,i) = rr + R(i+1,i) = 0d0 + if (i < n) then + call zrot(n-i,R(i,i+1),ldr,R(i+1,i+1),ldr,c,s) + end if +c apply rotation to Q + if (m > 0) then + call zrot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s)) + end if + end do + end + diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/zqr1up.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqr1up.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,53 @@ +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 zqr1up(m,n,k,Q,R,u,v) +c purpose: updates a QR factorization after rank-1 modification +c i.e., given a m-by-k unitary Q and m-by-n upper +c trapezoidal R, an m-vector u and n-vector v, +c this subroutine updates Q -> Q1 and R -> R1 so that +c Q1*R1 = Q*R + Q*Q'u*v', and Q1 is again unitary +c and R1 upper trapezoidal. +c (complex version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q, and rows of R. k <= m. +c Q (io) on entry, the unitary m-by-k matrix Q. +c on exit, the updated matrix Q1. +c R (io) on entry, the upper trapezoidal m-by-n matrix R. +c on exit, the updated matrix R1. +c u (in) the left m-vector. +c v (in) the right n-vector. +c + integer m,n,k + double complex Q(m,k),R(k,n),u(m),v(n) + double complex w + external zqrqhv,zqhqr + integer i +c quick return if possible + if (m <= 0 .or. n <= 0) return +c eliminate tail of Q'*u + call zqrqhv(m,n,k,Q,m,R,m,u,w) +c update R + do i = 1,n + R(1,i) = R(1,i) + w*conjg(v(i)) + end do +c retriangularize R + call zqhqr(m,n,k,Q,m,R,k) + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/zqrdec.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqrdec.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,66 @@ +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 zqrdec(m,n,k,Q,R,R1,j) +c purpose: updates a QR factorization after deleting +c a column. +c i.e., given an m-by-k unitary matrix Q, an k-by-n +c upper trapezoidal matrix R and index j in the range +c 1:n+1, this subroutine updates the matrix Q -> Q1 and +c forms an m-by-(n-1) matrix R1 so that Q1 remains +c unitary, R1 is upper trapezoidal, and +c Q1*R1 = [A(:,1:j-1) A(:,j+1:n)], where A = Q*R. +c (complex version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q, and rows of R. +c Q (io) on entry, the unitary m-by-k matrix Q. +c on exit, the updated matrix Q1. +c R (in) the original upper trapezoidal matrix R. +c R1 (out) the updated matrix R1. +c j (in) the position of the deleted column in R. +c 1 <= j <= n. +c + integer m,n,k,j + double complex Q(m,k),R(k,n),R1(k,n-1) + external zcopy,zqhqr + integer info +c quick return if possible + if (m <= 0 .or. k <= 0 .or. n == 1) return +c check arguments + info = 0 + if (n < 1) then + info = 2 + else if (j < 1 .or. j > n) then + info = 7 + end if + if (info /= 0) then + call xerbla('DQRDEC',info) + end if +c copy leading portion + call zcopy(k*(j-1),R,1,R1,1) + if (j < n) then +c copy trailing portion of R + call zcopy(k*(n-j),R(1,j+1),1,R1(1,j),1) +c if necessary, retriangularize R1(j:k,j:n-1) and update Q(:,j:k) + if (j < k) then + call zqhqr(m,n-j,k-j+1,Q(1,j),m,R1(j,j),k) + end if + end if + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/zqrder.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqrder.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,93 @@ +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 zqrder(m,n,Q,Q1,R,R1,j) +c purpose: updates a QR factorization after deleting a row. +c i.e., given an m-by-m unitary matrix Q, an m-by-n +c upper trapezoidal matrix R and index j in the range +c 1:m, this subroutine forms the (m-1)-by-(m-1) matrix +c Q1 and an (m-1)-by-n matrix R1 so that Q1 is again +c unitary, R1 upper trapezoidal, and +c Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R. +c (complex version) +c +c arguments: +c m (in) number of rows of the matrix R. +c n (in) number of columns of the matrix R +c Q (in) the unitary matrix Q +c Q1 (out) the updated matrix Q1 +c R (in) the upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the new row in R1 +c + integer m,n,j + double complex Q(m,m),Q1(m-1,m-1),R(m,n),R1(m-1,n) + double precision c + double complex s,rr,w + external xerbla,zlacpy,zcopy,zlartg,zrot,zdscal,zaxpy + integer i +c quick return if possible + if (m == 1) return +c check arguments + info = 0 + if (m < 1) then + info = 1 + else if (j < 1 .or. j > n) then + info = 7 + end if + if (info /= 0) then + call xerbla('ZQRDER',info) + end if +c setup the new matrix Q1 +c permute the columns of Q and rows of R so that the deleted row ends +c up being the topmost row. + if (j > 1) then + call zlacpy('0',j-1,m-1,Q(1,2),m,Q1(1,1),m-1) + end if + if (j < m) then + call zlacpy('0',m-j,m-1,Q(j+1,2),m,Q1(j,1),m-1) + end if +c setup the new matrix R1 + call zlacpy('0',m-1,n,R(2,1),m,R1(1,1),m-1) +c eliminate Q(j,2:m) + w = Q(j,m) + do i = m-1,2,-1 + call zlartg(Q(j,i),w,c,s,rr) + w = rr +c apply rotation to rows of R1 + if (i <= n) then + call zrot(n-i+1,R1(i-1,i),m-1,R1(i,i),m-1,c,conjg(s)) + end if +c apply rotation to columns of Q1 + call zrot(m-1,Q1(1,i-1),1,Q1(1,i),1,c,s) + end do +c the last iteration is special, as we don't have the first row of +c R and first column of Q + call zlartg(Q(j,1),w,c,s,rr) + w = rr + call zdscal(n,c,R1(1,1),m-1) + call zaxpy(n,-s,R(1,1),m,R1(1,1),m-1) +c apply rotation to columns of Q1 + call zdscal(m-1,c,Q1(1,1),1) + if (j > 1) then + call zaxpy(j-1,-conjg(s),Q(1,1),1,Q1(1,1),1) + end if + if (j < m) then + call zaxpy(m-j,-conjg(s),Q(j+1,1),1,Q1(j,1),1) + end if + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/zqrinc.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqrinc.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,73 @@ +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 zqrinc(m,n,k,Q,R,R1,j,x) +c purpose: updates a QR factorization after inserting a new +c column. +c i.e., given an m-by-k unitary matrix Q, an m-by-n +c upper trapezoidal matrix R and index j in the range +c 1:n+1, this subroutine updates the matrix Q -> Q1 and +c forms an m-by-(n+1) matrix R1 so that Q1 is again unitary, +c R1 upper trapezoidal, and +c Q1*R1 = [A(:,1:j-1); Q*Q'*x; A(:,j:n-1)], where A = Q*R. +c (complex version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q, and rows of R. k <= m. +c Q (io) on entry, the unitary matrix Q. +c on exit, the updated matrix Q1 +c R (in) the original upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the new column in R1 +c x (in) the column being inserted +c + integer m,n,k,j + double complex Q(m,k),R(k,n),R1(k,n+1),x(m) + double complex w + external xerbla,zcopy,zqrqhv + integer info,i,jj +c quick return if possible + if (m <= 0) return +c check arguments + info = 0 + if (n < 0) then + info = 2 + else if (j < 1 .or. j > n+1) then + info = 6 + end if + if (info /= 0) then + call xerbla('ZQRDER',info) + end if +c copy leading portion of R + call zcopy(k*(j-1),R,1,R1,1) + if (j <= n) then + call zcopy(k*(n+1-j),R(1,j),1,R1(1,j+1),1) + end if + call zgemv('C',m,min(k,j-1),1d0,Q,m,x,1,0d0,R1(1,j),1) + if (j < k) then +c eliminate tail, updating Q(:,j:k) and R1(j:k,j+1:n+1) + jj = min(j,n)+1 + call zqrqhv(m,n+1-j,k-j+1,Q(1,j),m,R1(j,jj),m,x,w) +c assemble inserted column + R1(j,j) = w + do i = j+1,k + R1(i,j) = 0d0 + end do + end if + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/zqrinr.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqrinr.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,73 @@ +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 zqrinr(m,n,Q,Q1,R,R1,j,x) +c purpose: updates a QR factorization after inserting a new +c row. +c i.e., given an m-by-m unitary matrix Q, an m-by-n +c upper trapezoidal matrix R and index j in the range +c 1:m+1, this subroutine forms the (m+1)-by-(m+1) matrix +c Q1 and an (m+1)-by-n matrix R1 so that Q1 is again +c unitary, R1 upper trapezoidal, and +c Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R. +c (complex version) +c arguments: +c m (in) number of rows of the matrix R. +c n (in) number of columns of the matrix R +c Q (in) the orthogonal matrix Q +c Q1 (out) the updated matrix Q1 +c R (in) the upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the new row in R1 +c x (in) the row being added +c + integer m,n,j + double complex Q(m,m),Q1(m+1,m+1),R(m,n),R1(m+1,n),x(n) + external xerbla,zlacpy,dcopy,dqhqr + integer i +c check arguments + info = 0 + if (n < 0) then + info = 2 + else if (j < 1 .or. j > m+1) then + info = 7 + end if + if (info /= 0) then + call xerbla('ZQRINR',info) + end if +c setup the new matrix Q1 +c permute the columns of Q1 and rows of R1 so that c the new row ends +c up being the topmost row. + if (j > 1) then + call zlacpy('0',j-1,m,Q(1,1),m,Q1(1,2),m+1) + end if + if (j <= m) then + call zlacpy('0',m-j+1,m,Q(j,1),m,Q1(j+1,2),m+1) + end if +c zero the rest of Q1 + do i = 1,m+1 + Q1(i,1) = 0d0 + Q1(j,i) = 0d0 + end do + Q1(j,1) = 1d0 +c setup the new matrix R1 + call zcopy(n,x,1,R1(1,1),m+1) + call zlacpy('0',m,n,R(1,1),m,R1(2,1),m+1) +c rotate to form proper QR + call zqhqr(m+1,n,m+1,Q1,m+1,R1,m+1) + end diff -r 8a6965a01176 -r 205ddddc2990 libcruft/qrupdate/zqrqhv.f --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqrqhv.f Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,75 @@ +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 zqrqhv(m,n,k,Q,ldq,R,ldr,u,rr) +c purpose: given an m-by-k matrix Q, an upper trapezoidal +c k-by-n matrix R, and an m-vector u, this subroutine +c updates the matrices Q -> Q1 and R -> R1 so that +c Q1 = Q*G', R1 = G*R, w1(2:m) = 0 with G unitary, +c R1 upper Hessenberg, and w1 = Q1'*u. +c (complex version) +c arguments: +c m (in) number of rows of the matrix Q. +c n (in) number of columns of the matrix R. +c k (in) number of columns of Q and rows of R. k <= m. +c Q (io) on entry, the orthogonal matrix Q. +c on exit, the updated matrix Q1. +c ldq (in) leading dimension of Q. +c R (io) on entry, the upper triangular matrix R. +c on exit, the updated upper Hessenberg matrix R1. +c ldr (in) leading dimension of R. +c u (in) the m-vector u. +c rr (out) the first element of Q1'*u on exit. +c +c if Q is orthogonal, so is Q1. It is not strictly +c necessary, however. + integer m,n,k,ldq,ldr + double complex Q(ldq,*),R(ldr,*),u(*),rr + double precision c + double complex s,w,w1,zdotc + external zdotc,zlartg,zrot + integer i,info +c quick return if possible. + if (k <= 0) return +c check arguments. + info = 0 + if (k > m) then + info = 3 + else if (ldq < 1) then + info = 5 + else if (ldr < 1) then + info = 7 + end if + if (info /= 0) then + call xerbla('ZQRQHV',info) + end if +c form each element of w = Q'*u when necessary. + rr = zdotc(m,Q(1,k),1,u,1) + do i = k-1,1,-1 + w1 = rr + w = zdotc(m,Q(1,i),1,u,1) + call zlartg(w,w1,c,s,rr) +c apply rotation to rows of R if necessary + if (i <= n) then + call zrot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s) + end if +c apply rotation to columns of Q + call zrot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s)) + end do + end + diff -r 8a6965a01176 -r 205ddddc2990 liboctave/ChangeLog --- a/liboctave/ChangeLog Sun Feb 24 03:32:49 2008 -0500 +++ b/liboctave/ChangeLog Sun Feb 24 13:27:33 2008 +0100 @@ -1,3 +1,8 @@ 2008-02-24 John W. Eaton + + * dbleQR.h, dbleQR.cc, CmplxQR.h, CmplxQR.cc: + Implement rank-1 updating methods for QR classes. + 2008-02-24 John W. Eaton * oct-inttypes.h (octave_int_helper): New class. Provide diff -r 8a6965a01176 -r 205ddddc2990 liboctave/CmplxQR.cc --- a/liboctave/CmplxQR.cc Sun Feb 24 03:32:49 2008 -0500 +++ b/liboctave/CmplxQR.cc Sun Feb 24 13:27:33 2008 +0100 @@ -21,6 +21,8 @@ along with Octave; see the file COPYING. */ +// updating/downdating by Jaroslav Hajek 2008 + #ifdef HAVE_CONFIG_H #include #endif @@ -28,6 +30,8 @@ along with Octave; see the file COPYING. #include "CmplxQR.h" #include "f77-fcn.h" #include "lo-error.h" +#include "Range.h" +#include "idx-vector.h" extern "C" { @@ -40,6 +44,30 @@ extern "C" F77_FUNC (zungqr, ZUNGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, Complex*, const octave_idx_type&, Complex*, Complex*, const octave_idx_type&, octave_idx_type&); + + // these come from qrupdate + + F77_RET_T + F77_FUNC (zqr1up, ZQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + Complex*, Complex*, const Complex*, const Complex*); + + F77_RET_T + F77_FUNC (zqrinc, ZQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + Complex*, const Complex*, Complex*, const octave_idx_type&, const Complex*); + + F77_RET_T + F77_FUNC (zqrdec, ZQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + Complex*, const Complex*, Complex*, const octave_idx_type&); + + F77_RET_T + F77_FUNC (zqrinr, ZQRINR) (const octave_idx_type&, const octave_idx_type&, + const Complex*, Complex*, const Complex*, Complex*, + const octave_idx_type&, const Complex*); + + F77_RET_T + F77_FUNC (zqrder, ZQRDER) (const octave_idx_type&, const octave_idx_type&, + const Complex*, Complex*, const Complex*, Complex *, + const octave_idx_type&); } ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type) @@ -127,6 +155,142 @@ ComplexQR::init (const ComplexMatrix& a, } } + +ComplexQR::ComplexQR (const ComplexMatrix &q, const ComplexMatrix& r) +{ + if (q.columns () != r.rows ()) + { + (*current_liboctave_error_handler) ("QR dimensions mismatch"); + return; + } + + this->q = q; + this->r = r; +} + +void ComplexQR::update (const ComplexMatrix &u, const ComplexMatrix &v) +{ + int m = q.rows (), n = r.columns (), k = q.columns (); + + if (u.length () != m || v.length () != n) + { + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + return; + } + + F77_XFCN(zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), + u.data (), v.data () )); + +} + +void ComplexQR::insert_col (const ComplexMatrix &u, octave_idx_type j) +{ + + int m = q.rows (), n = r.columns (), k = q.columns (); + + if (u.length () != m) + { + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + return; + } + + if (j < 1 || j > n+1) + { + (*current_liboctave_error_handler) ("QR insert index out of range"); + return; + } + + ComplexMatrix r1 (m,n+1); + + F77_XFCN(zqrinc, ZQRINC, (m, n, k, q.fortran_vec (), r.data (), r1.fortran_vec (), + j, u.data () )); + + r = r1; +} + +void ComplexQR::delete_col (octave_idx_type j) +{ + int m = q.rows (), k = r.rows (), n = r.columns (); + + if (k < m && k < n) + { + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + return; + } + + if (j < 1 || j > n) + { + (*current_liboctave_error_handler) ("QR delete index out of range"); + return; + } + + ComplexMatrix r1 (k,n-1); + + F77_XFCN(zqrdec, ZQRDEC, (m, n, k, q.fortran_vec (), r.data (), r1.fortran_vec (), + j )); + + r = r1; +} + +void ComplexQR::insert_row (const ComplexMatrix &u, octave_idx_type j) +{ + int m = r.rows (), n = r.columns (); + + if (!q.is_square () || u.length () != n) + { + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + return; + } + + if (j < 1 || j > m+1) + { + (*current_liboctave_error_handler) ("QR insert index out of range"); + return; + } + + ComplexMatrix q1 (m+1,m+1); + ComplexMatrix r1 (m+1,n); + + F77_XFCN(zqrinr, ZQRINR, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j, u.data () )); + + q = q1; + r = r1; +} + +void ComplexQR::delete_row (octave_idx_type j) +{ + int m = r.rows (), n = r.columns (); + + if (!q.is_square ()) + { + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + return; + } + + if (j < 1 || j > n) + { + (*current_liboctave_error_handler) ("QR delete index out of range"); + return; + } + + ComplexMatrix q1 (m-1,m-1); + ComplexMatrix r1 (m-1,n); + + F77_XFCN(zqrder, ZQRDER, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j )); + + q = q1; + r = r1; +} + +void ComplexQR::economize () +{ + idx_vector c (':'), i (Range (1, r.columns ())); + q = ComplexMatrix (q.index (c, i)); + r = ComplexMatrix (r.index (i, c)); +} + /* ;;; Local Variables: *** ;;; mode: C++ *** diff -r 8a6965a01176 -r 205ddddc2990 liboctave/CmplxQR.h --- a/liboctave/CmplxQR.h Sun Feb 24 03:32:49 2008 -0500 +++ b/liboctave/CmplxQR.h Sun Feb 24 13:27:33 2008 +0100 @@ -27,7 +27,10 @@ along with Octave; see the file COPYING. #include #include "CMatrix.h" +#include "CColVector.h" +#include "CRowVector.h" #include "dbleQR.h" + class OCTAVE_API @@ -38,6 +41,8 @@ public: ComplexQR (void) : q (), r () { } ComplexQR (const ComplexMatrix&, QR::type = QR::std); + + ComplexQR (const ComplexMatrix& q, const ComplexMatrix& r); ComplexQR (const ComplexQR& a) : q (a.q), r (a.r) { } @@ -59,6 +64,18 @@ public: ComplexMatrix R (void) const { return r; } + void update (const ComplexMatrix &u, const ComplexMatrix &v); + + void insert_col (const ComplexMatrix &u, octave_idx_type j); + + void delete_col (octave_idx_type j); + + void insert_row (const ComplexMatrix &u, octave_idx_type j); + + void delete_row (octave_idx_type j); + + void economize(); + friend std::ostream& operator << (std::ostream&, const ComplexQR&); protected: diff -r 8a6965a01176 -r 205ddddc2990 liboctave/dbleQR.cc --- a/liboctave/dbleQR.cc Sun Feb 24 03:32:49 2008 -0500 +++ b/liboctave/dbleQR.cc Sun Feb 24 13:27:33 2008 +0100 @@ -21,6 +21,8 @@ along with Octave; see the file COPYING. */ +// updating/downdating by Jaroslav Hajek 2008 + #ifdef HAVE_CONFIG_H #include #endif @@ -28,6 +30,8 @@ along with Octave; see the file COPYING. #include "dbleQR.h" #include "f77-fcn.h" #include "lo-error.h" +#include "Range.h" +#include "idx-vector.h" extern "C" { @@ -38,6 +42,30 @@ extern "C" F77_RET_T F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, double*, const octave_idx_type&, double*, double*, const octave_idx_type&, octave_idx_type&); + + // these come from qrupdate + + F77_RET_T + F77_FUNC (dqr1up, DQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + double*, double*, const double*, const double*); + + F77_RET_T + F77_FUNC (dqrinc, DQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + double*, const double*, double*, const octave_idx_type&, const double*); + + F77_RET_T + F77_FUNC (dqrdec, DQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + double*, const double*, double*, const octave_idx_type&); + + F77_RET_T + F77_FUNC (dqrinr, DQRINR) (const octave_idx_type&, const octave_idx_type&, + const double*, double*, const double*, double*, + const octave_idx_type&, const double*); + + F77_RET_T + F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&, + const double*, double*, const double*, double *, + const octave_idx_type&); } QR::QR (const Matrix& a, QR::type qr_type) @@ -118,6 +146,142 @@ QR::init (const Matrix& a, QR::type qr_t } } +QR::QR (const Matrix &q, const Matrix& r) +{ + if (q.columns () != r.rows ()) + { + (*current_liboctave_error_handler) ("QR dimensions mismatch"); + return; + } + + this->q = q; + this->r = r; +} + +void QR::update (const Matrix &u, const Matrix &v) +{ + int m = q.rows (), n = r.columns (), k = q.columns (); + + if (u.length () != m || v.length () != n) + { + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + return; + } + + F77_XFCN(dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), + u.data (), v.data () )); + +} + +void QR::insert_col (const Matrix &u, octave_idx_type j) +{ + + int m = q.rows (), n = r.columns (), k = q.columns (); + + if (u.length () != m) + { + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + return; + } + + if (j < 1 || j > n+1) + { + (*current_liboctave_error_handler) ("QR insert index out of range"); + return; + } + + Matrix r1 (m,n+1); + + F77_XFCN(dqrinc, DQRINC, (m, n, k, q.fortran_vec (), r.data (), r1.fortran_vec (), + j, u.data () )); + + r = r1; +} + +void QR::delete_col (octave_idx_type j) +{ + int m = q.rows (), k = r.rows (), n = r.columns (); + + if (k < m && k < n) + { + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + return; + } + + if (j < 1 || j > n) + { + (*current_liboctave_error_handler) ("QR delete index out of range"); + return; + } + + Matrix r1 (k,n-1); + + F77_XFCN(dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), r.data (), r1.fortran_vec (), + j )); + + r = r1; +} + +void QR::insert_row (const Matrix &u, octave_idx_type j) +{ + int m = r.rows (), n = r.columns (); + + if (!q.is_square () || u.length () != n) + { + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + return; + } + + if (j < 1 || j > m+1) + { + (*current_liboctave_error_handler) ("QR insert index out of range"); + return; + } + + Matrix q1 (m+1,m+1); + Matrix r1 (m+1,n); + + F77_XFCN(dqrinr, DQRINR, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j, u.data () )); + + q = q1; + r = r1; +} + +void QR::delete_row (octave_idx_type j) +{ + int m = r.rows (), n = r.columns (); + + if (!q.is_square ()) + { + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + return; + } + + if (j < 1 || j > n) + { + (*current_liboctave_error_handler) ("QR delete index out of range"); + return; + } + + Matrix q1 (m-1,m-1); + Matrix r1 (m-1,n); + + F77_XFCN(dqrder, DQRDER, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j )); + + q = q1; + r = r1; +} + +void QR::economize () +{ + idx_vector c (':'), i (Range (1, r.columns ())); + q = Matrix (q.index (c, i)); + r = Matrix (r.index (i, c)); +} + + /* ;;; Local Variables: *** ;;; mode: C++ *** diff -r 8a6965a01176 -r 205ddddc2990 liboctave/dbleQR.h --- a/liboctave/dbleQR.h Sun Feb 24 03:32:49 2008 -0500 +++ b/liboctave/dbleQR.h Sun Feb 24 13:27:33 2008 +0100 @@ -21,12 +21,16 @@ along with Octave; see the file COPYING. */ +// updating/downdating by Jaroslav Hajek 2008 + #if !defined (octave_QR_h) #define octave_QR_h 1 #include #include "dMatrix.h" +#include "dColVector.h" +#include "dRowVector.h" class OCTAVE_API @@ -44,6 +48,8 @@ public: QR (void) : q (), r () { } QR (const Matrix&, QR::type = QR::std); + + QR (const Matrix& q, const Matrix& r); QR (const QR& a) : q (a.q), r (a.r) { } @@ -65,6 +71,18 @@ public: Matrix R (void) const { return r; } + void update (const Matrix &u, const Matrix &v); + + void insert_col (const Matrix &u, octave_idx_type j); + + void delete_col (octave_idx_type j); + + void insert_row (const Matrix &u, octave_idx_type j); + + void delete_row (octave_idx_type j); + + void economize(); + friend std::ostream& operator << (std::ostream&, const QR&); protected: diff -r 8a6965a01176 -r 205ddddc2990 src/ChangeLog --- a/src/ChangeLog Sun Feb 24 03:32:49 2008 -0500 +++ b/src/ChangeLog Sun Feb 24 13:27:33 2008 +0100 @@ -1,3 +1,9 @@ 2008-02-22 John W. Eaton + + * DLD-FUNCTIONS/qrupdate.cc, DLD-FUNCTIONS/qrinsert.cc, + DLD-FUNCTIONS/qrdelete.cc: add new sources + * Makefile.in: include the new sources in build process + 2008-02-22 John W. Eaton * DLD-FUNCTIONS/ccolamd.cc, DLD-FUNCTIONS/colamd.cc, diff -r 8a6965a01176 -r 205ddddc2990 src/DLD-FUNCTIONS/qrdelete.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/qrdelete.cc Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,189 @@ +/* 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 "dbleQR.h" +#include "CmplxQR.h" + +#include "defun-dld.h" +#include "error.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (qrdelete, args, nargout, +"-*- texinfo -*-\n\ address@hidden Function} \ address@hidden,@var{R1}]} = qrdelete(@var{Q},@var{R},@var{j},@var{orient})\n\ +Given a address@hidden of a real or complex matrix @address@hidden = @address@hidden,\n\ address@hidden@tie{}unitary and @address@hidden trapezoidal,\n\ +this function returns the address@hidden of @w{[A(:,1:j-1) A(:,j+1:n)]},\n\ +i.e. A with one column deleted (if @var{orient} == 'col'),\n\ +or the address@hidden of @w{[A(1:j-1,:);A(:,j+1:n)]},\n\ +i.e. A with one row deleted (if @var{orient} == 'row').\n\ address@hidden defaults to 'col'.\n\ +\n\ +If @var{orient} == 'col',\n\ +the matrix @var{Q} is not required to be square.\n\ +For Matlab compatibility, if @var{Q} is nonsquare on input, the updated\n\ +factorization is always stripped to the economical form, i.e.\n\ address@hidden(Q,2) == size(R,1) = size(R,2)}.\n\ +To get the less intelligent but more natural behaviour when Q retains it shape\n\ +and R loses one column, use 'col+' is @var{orient} instead.\n\ +If @var{orient} == 'row', @var{Q} must be square.\n\ address@hidden, qrinsert, qrupdate}\n\ address@hidden deftypefn") +{ + int nargin = args.length (); + octave_value_list retval; + + octave_value argQ,argR,argj,argor; + + if ((nargin == 3 || nargin == 4) && nargout == 2 + && (argQ = args (0), argQ.is_matrix_type ()) + && (argR = args (1), argR.is_matrix_type ()) + && (argj = args (2), argj.is_scalar_type ()) + && (nargin < 4 || (argor = args (3), argor.is_string ()))) + { + int m = argQ.rows (), k = argQ.columns (), n = argR.columns (); + bool row = false, colp = false; + std::string orient = (nargin < 4) ? "col" : argor.string_value (); + + if (nargin < 4 || (row = orient == "row") + || (orient == "col") || (colp = orient == "col+")) + if (argR.rows() == k) + { + int j = argj.scalar_value (); + if (j >= 1 && j <= (row ? n : m)+1) + { + if ( argQ.is_real_matrix () + && argR.is_real_matrix () ) + { + // real case + Matrix Q = argQ.matrix_value (); + Matrix R = argR.matrix_value (); + + QR fact (Q, R); + + if (row) + fact.delete_row(j); + else + { + fact.delete_col(j); + + if (!colp && k < m) + fact.economize (); + } + + retval (0) = fact.Q (); + retval (1) = fact.R (); + } + else + { + // complex case + ComplexMatrix Q = argQ.complex_matrix_value (); + ComplexMatrix R = argR.complex_matrix_value (); + + ComplexQR fact (Q, R); + + if (row) + fact.delete_row (j); + else + { + fact.delete_col(j); + + if (!colp && k < m) + fact.economize (); + } + + retval (0) = fact.Q (); + retval (1) = fact.R (); + } + + } + else + error ("qrdelete: index j out of range"); + } + else + error ("qrdelete: dimension mismatch"); + + else + error ("qrdelete: orient must be \"col\" or \"row\""); + } + else + print_usage (); + + return retval; +} + +/* +%!test +%! A = [0.091364 0.613038 0.027504 0.999083; +%! 0.594638 0.425302 0.562834 0.603537; +%! 0.383594 0.291238 0.742073 0.085574; +%! 0.265712 0.268003 0.783553 0.238409; +%! 0.669966 0.743851 0.457255 0.445057 ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrdelete(Q,R,3); +%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(:,1:2) A(:,4)]),Inf) < norm(A)*1e1*eps) +%! +%!test +%! A = [0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i; +%! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i; +%! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i; +%! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i; +%! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ] * I; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrdelete(Q,R,3); +%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(:,1:2) A(:,4)]),Inf) < norm(A)*1e1*eps) +%! +%!test +%! A = [0.091364 0.613038 0.027504 0.999083; +%! 0.594638 0.425302 0.562834 0.603537; +%! 0.383594 0.291238 0.742073 0.085574; +%! 0.265712 0.268003 0.783553 0.238409; +%! 0.669966 0.743851 0.457255 0.445057 ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrdelete(Q,R,3,'row'); +%! assert(norm(vec(Q'*Q - eye(4)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(1:2,:);A(4:5,:)]),Inf) < norm(A)*1e1*eps) +%! +%!test +%! A = [0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i; +%! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i; +%! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i; +%! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i; +%! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ] * I; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrdelete(Q,R,3,'row'); +%! assert(norm(vec(Q'*Q - eye(4)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(1:2,:);A(4:5,:)]),Inf) < norm(A)*1e1*eps) +*/ diff -r 8a6965a01176 -r 205ddddc2990 src/DLD-FUNCTIONS/qrinsert.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/qrinsert.cc Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,197 @@ +/* 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 "dbleQR.h" +#include "CmplxQR.h" + +#include "defun-dld.h" +#include "error.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (qrinsert, args, nargout, +"-*- texinfo -*-\n\ address@hidden Function} \ address@hidden,@var{R1}]} = qrinsert(@var{Q},@var{R},@var{j},@var{x},@var{orient})\n\ +Given a address@hidden of a real or complex matrix @address@hidden = @address@hidden,\n\ address@hidden@tie{}unitary and @address@hidden trapezoidal,\n\ +this function returns the address@hidden of @w{[A(:,1:j-1) x A(:,j:n)]},\n\ +where @var{u} is a column vector to be inserted into A (if @var{orient} == 'col'),\n\ +or the address@hidden of @w{[A(1:j-1,:);x;A(:,j:n)]},\n\ +where @var{x} is a row vector to be inserted into A (if @var{orient} == 'row').\n\ address@hidden defaults to 'col'.\n\ +\n\ +If @var{orient} == 'col' and the matrix @var{Q} is not square, then what gets inserted\n\ +is the projection of @var{u} onto the space spanned by columns of @var{Q},\n\ +i.e. Q*Q'*u.\n\ +If @var{orient} == 'row', @var{Q} must be square.\n\ address@hidden, qrupdate, qrdelete}\n\ address@hidden deftypefn") +{ + int nargin = args.length (); + octave_value_list retval; + + octave_value argQ,argR,argj,argx,argor; + + if ((nargin == 4 || nargin == 5) && nargout == 2 + && (argQ = args (0), argQ.is_matrix_type ()) + && (argR = args (1), argR.is_matrix_type ()) + && (argj = args (2), argj.is_scalar_type ()) + && (argx = args (3), argx.is_matrix_type ()) + && (nargin < 5 || (argor = args (4), argor.is_string ()))) + { + int m = argQ.rows (), n = argR.columns (), k = argQ.columns(); + bool row = false; + std::string orient = (nargin < 5) ? "col" : argor.string_value (); + + if (nargin < 5 || (row = orient == "row") || (orient == "col")) + if (argR.rows () == k + && (!row || m == k) + && argx.rows () == (row ? 1 : m) + && argx.columns () == (row ? n : 1)) + { + int j = argj.scalar_value (); + if (j >= 1 && j <= (row ? n : m)+1) + { + if ( argQ.is_real_matrix () + && argR.is_real_matrix () + && argx.is_real_matrix () ) + { + // real case + Matrix Q = argQ.matrix_value (); + Matrix R = argR.matrix_value (); + Matrix x = argx.matrix_value (); + + QR fact (Q, R); + + if (row) + fact.insert_row(x, j); + else + fact.insert_col(x, j); + + retval (0) = fact.Q (); + retval (1) = fact.R (); + } + else + { + // complex case + ComplexMatrix Q = argQ.complex_matrix_value (); + ComplexMatrix R = argR.complex_matrix_value (); + ComplexMatrix x = argx.complex_matrix_value (); + + ComplexQR fact (Q, R); + + if (row) + fact.insert_row (x, j); + else + fact.insert_col (x, j); + + retval (0) = fact.Q (); + retval (1) = fact.R (); + } + + } + else + error ("qrinsert: index j out of range"); + } + else + error ("qrinsert: dimension mismatch"); + + else + error ("qrinsert: orient must be \"col\" or \"row\""); + } + else + print_usage (); + + return retval; +} + +/* +%!test +%! A = [0.091364 0.613038 0.999083; +%! 0.594638 0.425302 0.603537; +%! 0.383594 0.291238 0.085574; +%! 0.265712 0.268003 0.238409; +%! 0.669966 0.743851 0.445057 ]; +%! +%! x = [0.85082; +%! 0.76426; +%! 0.42883; +%! 0.53010; +%! 0.80683 ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrinsert(Q,R,3,x); +%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf) < norm(A)*1e1*eps) +%! +%!test +%! A = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i; +%! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i; +%! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i; +%! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i; +%! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ]; +%! +%! x = [0.20351 + 0.05401i; +%! 0.13141 + 0.43708i; +%! 0.29808 + 0.08789i; +%! 0.69821 + 0.38844i; +%! 0.74871 + 0.25821i ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrinsert(Q,R,3,x); +%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf) < norm(A)*1e1*eps) +%! +%!test +%! A = [0.091364 0.613038 0.999083; +%! 0.594638 0.425302 0.603537; +%! 0.383594 0.291238 0.085574; +%! 0.265712 0.268003 0.238409; +%! 0.669966 0.743851 0.445057 ]; +%! +%! x = [0.85082 0.76426 0.42883 ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrinsert(Q,R,3,x,'row'); +%! assert(norm(vec(Q'*Q - eye(6)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(1:2,:);x;A(3:5,:)]),Inf) < norm(A)*1e1*eps) +%! +%!test +%! A = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i; +%! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i; +%! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i; +%! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i; +%! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ]; +%! +%! x = [0.20351 + 0.05401i 0.13141 + 0.43708i 0.29808 + 0.08789i ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrinsert(Q,R,3,x,'row'); +%! assert(norm(vec(Q'*Q - eye(6)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - [A(1:2,:);x;A(3:5,:)]),Inf) < norm(A)*1e1*eps) +*/ diff -r 8a6965a01176 -r 205ddddc2990 src/DLD-FUNCTIONS/qrupdate.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/qrupdate.cc Sun Feb 24 13:27:33 2008 +0100 @@ -0,0 +1,152 @@ +/* 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 "dbleQR.h" +#include "CmplxQR.h" + +#include "defun-dld.h" +#include "error.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (qrupdate, args, nargout, +"-*- texinfo -*-\n\ address@hidden Function}\ address@hidden,@var{R1}]} = qrupdate(@var{Q},@var{R},@var{u},@var{v})\n\ +Given a address@hidden of a real or complex matrix @address@hidden = @address@hidden,\n\ address@hidden@tie{}unitary and @address@hidden trapezoidal,\n\ +this function returns the address@hidden of @address@hidden + @address@hidden'},\n\ +where @var{u} and @var{v} are column vectors (rank-1 update).\n\ +\n\ +If the matrix @var{Q} is not square, the matrix @var{A} gets updated by\n\ +Q*Q'*u*v' instead of u*v'.\n\ address@hidden, qrinsert, qrdelete}\n\ address@hidden deftypefn") +{ + int nargin = args.length (); + octave_value_list retval; + + octave_value argQ,argR,argu,argv; + + if (nargin == 4 && nargout == 2 + && (argQ = args (0),argQ.is_matrix_type ()) + && (argR = args (1),argR.is_matrix_type ()) + && (argu = args (2),argu.is_matrix_type ()) + && (argv = args (3),argv.is_matrix_type ())) + { + int m = argQ.rows (), n = argR.columns (), k = argQ.columns(); + if (argR.rows () == k + && argu.rows () == m && argu.columns () == 1 + && argv.rows () == n && argv.columns () == 1) + { + if ( argQ.is_real_matrix () + && argR.is_real_matrix () + && argu.is_real_matrix () + && argv.is_real_matrix () ) + { + // all real case + Matrix Q = argQ.matrix_value (); + Matrix R = argR.matrix_value (); + Matrix u = argu.matrix_value (); + Matrix v = argv.matrix_value (); + + QR fact (Q, R); + fact.update (u, v); + + retval (0) = fact.Q (); + retval (1) = fact.R (); + } + else + { + // complex case + ComplexMatrix Q = argQ.complex_matrix_value (); + ComplexMatrix R = argR.complex_matrix_value (); + ComplexMatrix u = argu.complex_matrix_value (); + ComplexMatrix v = argv.complex_matrix_value (); + + ComplexQR fact (Q, R); + fact.update (u, v); + + retval (0) = fact.Q (); + retval (1) = fact.R (); + } + } + else + { + error ("qrupdate: dimensions mismatch"); + } + } + else + { + print_usage(); + } + + return retval; +} +/* +%!test +%! A = [0.091364 0.613038 0.999083; +%! 0.594638 0.425302 0.603537; +%! 0.383594 0.291238 0.085574; +%! 0.265712 0.268003 0.238409; +%! 0.669966 0.743851 0.445057 ]; +%! +%! u = [0.85082; +%! 0.76426; +%! 0.42883; +%! 0.53010; +%! 0.80683 ]; +%! +%! v = [0.98810; +%! 0.24295; +%! 0.43167 ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrupdate(Q,R,u,v); +%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - A - u*v'),Inf) < norm(A)*1e1*eps) +%! +%!test +%! A = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i; +%! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i; +%! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i; +%! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i; +%! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ]; +%! +%! u = [0.20351 + 0.05401i; +%! 0.13141 + 0.43708i; +%! 0.29808 + 0.08789i; +%! 0.69821 + 0.38844i; +%! 0.74871 + 0.25821i ]; +%! +%! v = [0.85839 + 0.29468i; +%! 0.20820 + 0.93090i; +%! 0.86184 + 0.34689i ]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrupdate(Q,R,u,v); +%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - A - u*v'),Inf) < norm(A)*1e1*eps) +*/ diff -r 8a6965a01176 -r 205ddddc2990 src/Makefile.in --- a/src/Makefile.in Sun Feb 24 03:32:49 2008 -0500 +++ b/src/Makefile.in Sun Feb 24 13:27:33 2008 +0100 @@ -69,6 +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 \ 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 \