Mercurial > hg > octave-lyh
changeset 7553:56be6f31dd4e
implementation of QR factorization updating
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Tue, 04 Mar 2008 21:47:11 -0500 |
| parents | 6070c3bd69c4 |
| children | 40574114c514 |
| files | configure.in libcruft/ChangeLog libcruft/Makefile.in libcruft/qrupdate/Makefile.in libcruft/qrupdate/dqhqr.f libcruft/qrupdate/dqr1up.f libcruft/qrupdate/dqrdec.f libcruft/qrupdate/dqrder.f libcruft/qrupdate/dqrinc.f libcruft/qrupdate/dqrinr.f libcruft/qrupdate/dqrqhv.f libcruft/qrupdate/zqhqr.f libcruft/qrupdate/zqr1up.f libcruft/qrupdate/zqrdec.f libcruft/qrupdate/zqrder.f libcruft/qrupdate/zqrinc.f libcruft/qrupdate/zqrinr.f libcruft/qrupdate/zqrqhv.f liboctave/ChangeLog liboctave/CmplxQR.cc liboctave/CmplxQR.h liboctave/dbleQR.cc liboctave/dbleQR.h src/ChangeLog src/DLD-FUNCTIONS/qr.cc |
| diffstat | 25 files changed, 1901 insertions(+), 5 deletions(-) [+] |
line wrap: on
line diff
--- a/configure.in +++ b/configure.in @@ -1834,6 +1834,7 @@ 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])
--- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,13 @@ +2008-02-24 Jaroslav Hajek <highegg@gmail.com> + + * qrupdate/Makefile.in, 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: + New files. + * Makefile.in (CRUFT_DIRS): Add qrupdate to the list. + 2008-02-14 John W. Eaton <jwe@octave.org> * misc/f77-fcn.h (F77_XFCN): Call octave_rethrow_exception here
--- a/libcruft/Makefile.in +++ b/libcruft/Makefile.in @@ -43,7 +43,7 @@ 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)
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/Makefile.in @@ -0,0 +1,39 @@ +# 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 +# <http://www.gnu.org/licenses/>. + +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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/dqhqr.f @@ -0,0 +1,69 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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 +
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/dqr1up.f @@ -0,0 +1,52 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/dqrdec.f @@ -0,0 +1,66 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/dqrder.f @@ -0,0 +1,93 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/dqrinc.f @@ -0,0 +1,73 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/dqrinr.f @@ -0,0 +1,73 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/dqrqhv.f @@ -0,0 +1,75 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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 +
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/zqhqr.f @@ -0,0 +1,69 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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 +
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/zqr1up.f @@ -0,0 +1,53 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/zqrdec.f @@ -0,0 +1,66 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/zqrder.f @@ -0,0 +1,93 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/zqrinc.f @@ -0,0 +1,73 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/zqrinr.f @@ -0,0 +1,73 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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
new file mode 100644 --- /dev/null +++ b/libcruft/qrupdate/zqrqhv.f @@ -0,0 +1,75 @@ +c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +c +c Author: Jaroslav Hajek <highegg@gmail.com> +c +c This source is free software; you can redistribute it and/or modify +c it under the terms of the GNU General Public License as published by +c the Free Software Foundation; either version 2 of the License, or +c (at your option) any later version. +c +c This program is distributed in the hope that it will be useful, +c but WITHOUT ANY WARRANTY; without even the implied warranty of +c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +c GNU General Public License for more details. +c +c You should have received a copy of the GNU General Public License +c along with this software; see the file COPYING. If not, see +c <http://www.gnu.org/licenses/>. +c + subroutine 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 +
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,16 @@ +2008-03-04 Jaroslav Hajek <highegg@gmail.com> + + * 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. + * dbleQR.h: Provide decls. + * CmplxQR.cc (ComplexQR::update, ComplexQR::insert_col, + ComplexQR::delete_col, ComplexQR::insert_row, + ComplexQR::delete_row): New methods. + (ComplexQR::ComplexQR (const ComplexMatrix&, const ComplexMAtrix&)): + New constructor. + * CmplxQR.h: Provide decls. + 2008-03-04 Jaroslav Hajek <highegg@gmail.com> * Array-C.cc, Sparse-C.cc: Include oct-sort.cc after definitions
--- a/liboctave/CmplxQR.cc +++ b/liboctave/CmplxQR.cc @@ -21,6 +21,8 @@ */ +// updating/downdating by Jaroslav Hajek 2008 + #ifdef HAVE_CONFIG_H #include <config.h> #endif @@ -28,6 +30,8 @@ #include "CmplxQR.h" #include "f77-fcn.h" #include "lo-error.h" +#include "Range.h" +#include "idx-vector.h" extern "C" { @@ -40,6 +44,30 @@ 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,140 @@ } } + +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) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () != m || v.length () != n) + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + else + 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) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 1 || j > n+1) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + 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) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (k < m && k < n) + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + else if (j < 1 || j > n) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + 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) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 1 || j > m+1) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + 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) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 1 || j > n) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + 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 (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) + { + q.resize (q.rows (), r_nc); + r.resize (r_nc, r_nc); + } +#endif +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/CmplxQR.h +++ b/liboctave/CmplxQR.h @@ -27,8 +27,11 @@ #include <iostream> #include "CMatrix.h" +#include "CColVector.h" +#include "CRowVector.h" #include "dbleQR.h" + class OCTAVE_API ComplexQR @@ -39,6 +42,8 @@ ComplexQR (const ComplexMatrix&, QR::type = QR::std); + ComplexQR (const ComplexMatrix& q, const ComplexMatrix& r); + ComplexQR (const ComplexQR& a) : q (a.q), r (a.r) { } ComplexQR& operator = (const ComplexQR& a) @@ -59,6 +64,18 @@ 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:
--- a/liboctave/dbleQR.cc +++ b/liboctave/dbleQR.cc @@ -21,6 +21,8 @@ */ +// updating/downdating by Jaroslav Hajek 2008 + #ifdef HAVE_CONFIG_H #include <config.h> #endif @@ -28,6 +30,8 @@ #include "dbleQR.h" #include "f77-fcn.h" #include "lo-error.h" +#include "Range.h" +#include "idx-vector.h" extern "C" { @@ -38,6 +42,30 @@ 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,140 @@ } } +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) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () == m && v.length () == n) + F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), + u.data (), v.data ())); + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); +} + +void +QR::insert_col (const Matrix& u, octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 1 || j > n+1) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + 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) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (k < m && k < n) + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + else if (j < 1 || j > n) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + 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) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 1 || j > m+1) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + 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) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 1 || j > n) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + 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 (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) + { + q.resize (q.rows (), r_nc); + r.resize (r_nc, r_nc); + } +#endif +} + + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/dbleQR.h +++ b/liboctave/dbleQR.h @@ -21,12 +21,16 @@ */ +// updating/downdating by Jaroslav Hajek 2008 + #if !defined (octave_QR_h) #define octave_QR_h 1 #include <iostream> #include "dMatrix.h" +#include "dColVector.h" +#include "dRowVector.h" class OCTAVE_API @@ -45,6 +49,8 @@ QR (const Matrix&, QR::type = QR::std); + QR (const Matrix& q, const Matrix& r); + QR (const QR& a) : q (a.q), r (a.r) { } QR& operator = (const QR& a) @@ -65,6 +71,18 @@ 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 (void); + friend std::ostream& operator << (std::ostream&, const QR&); protected:
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2008-02-24 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/qr.cc (Fqrupdate, Fqrinsert, Fqrdelete): + New functions. + 2008-03-04 Ryan Rusaw <rrusaw@gmail.com> * toplev.h (octave_call_stack::element): New static function.
--- a/src/DLD-FUNCTIONS/qr.cc +++ b/src/DLD-FUNCTIONS/qr.cc @@ -20,6 +20,10 @@ */ +// The qrupdate, qrinsert, and qrdelete functions were written by +// Jaroslav Hajek <highegg@gmail.com>, Copyright (C) 2008 VZLU +// Prague, a.s., Czech Republic. + #ifdef HAVE_CONFIG_H #include <config.h> #endif @@ -179,7 +183,7 @@ int nargin = args.length (); - if (nargin < 1 || nargin > (args(0).is_sparse_type() ? 3 : 2) || nargout > 3) + if (nargin < 1 || nargin > (args(0).is_sparse_type() ? 3 : 2)) { print_usage (); return retval; @@ -342,9 +346,7 @@ } } else - { - gripe_wrong_type_arg ("qr", arg); - } + gripe_wrong_type_arg ("qr", arg); } return retval; @@ -429,6 +431,470 @@ */ +DEFUN_DLD (qrupdate, args, , + "-*- texinfo -*-\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\ +of @w{@var{A} + @var{u}*@var{v}'}, where @var{u} and @var{v} are\n\ +column vectors (rank-1 update).\n\ +\n\ +If the matrix @var{Q} is not square, the matrix @var{A} is updated by\n\ +Q*Q'*u*v' instead of u*v'.\n\ +@seealso{qr, qrinsert, qrdelete}\n\ +@end deftypefn") +{ + octave_idx_type nargin = args.length (); + octave_value_list retval; + + octave_value argQ,argR,argu,argv; + + if (nargin == 4 + && (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 ())) + { + octave_idx_type m = argQ.rows (); + octave_idx_type n = argR.columns (); + octave_idx_type 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(1) = fact.R (); + retval(0) = fact.Q (); + } + 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(1) = fact.R (); + retval(0) = fact.Q (); + } + } + 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) +*/ + +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\ +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\ +@w{[A(:,1:j-1) x A(:,j:n)]}, where @var{u} is a column vector to be\n\ +inserted into @var{A} (if @var{orient} is @code{\"col\"}), or the\n\ +QR@tie{}factorization of @w{[A(1:j-1,:);x;A(:,j:n)]}, where @var{x}\n\ +is a row vector to be inserted into @var{A} (if @var{orient} is\n\ +@code{\"row\"}).\n\ +\n\ +The default value of @var{orient} is @code{\"col\"}.\n\ +\n\ +If @var{orient} is @code{\"col\"} and the matrix @var{Q} is not square,\n\ +then what gets inserted is the projection of @var{u} onto the space\n\ +spanned by columns of @var{Q}, i.e. Q*Q'*u.\n\ +\n\ +If @var{orient} is @code{\"row\"}, @var{Q} must be square.\n\ +@seealso{qr, qrupdate, qrdelete}\n\ +@end deftypefn") +{ + octave_idx_type nargin = args.length (); + octave_value_list retval; + + octave_value argQ,argR,argj,argx,argor; + + if ((nargin == 4 || nargin == 5) + && (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 ()))) + { + octave_idx_type m = argQ.rows (); + octave_idx_type n = argR.columns (); + octave_idx_type 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.int_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(1) = fact.R (); + retval(0) = fact.Q (); + } + 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(1) = fact.R (); + retval(0) = fact.Q (); + } + + } + 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) +*/ + +DEFUN_DLD (qrdelete, args, , + "-*- texinfo -*-\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\ +@w{[A(:,1:j-1) A(:,j+1:n)]}, i.e. @var{A} with one column deleted\n\ +(if @var{orient} is \"col\"), or the QR@tie{}factorization of\n\ +@w{[A(1:j-1,:);A(:,j+1:n)]}, i.e. @var{A} with one row deleted (if\n\ +@var{orient} is \"row\").\n\ +\n\ +The default value of @var{orient} is \"col\".\n\ +\n\ +If @var{orient} is \"col\", the matrix @var{Q} is not required to\n\ +be square.\n\ +\n\ +For @sc{Matlab} compatibility, if @var{Q} is nonsquare on input, the\n\ +updated factorization is always stripped to the economical form, i.e.\n\ +@code{columns (Q) == rows (R) = columns (R)}.\n\ +\n\ +To get the less intelligent but more natural behaviour when @var{Q}\n\ +retains it shape and @var{R} loses one column, set @var{orient} to\n\ +\"col+\" instead.\n\ +\n\ +If @var{orient} is \"row\", @var{Q} must be square.\n\ +@seealso{qr, qrinsert, qrupdate}\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value_list retval; + + octave_value argQ,argR,argj,argor; + + if ((nargin == 3 || nargin == 4) + && (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 ()))) + { + octave_idx_type m = argQ.rows (); + octave_idx_type k = argQ.columns (); + octave_idx_type n = argR.columns (); + + bool row = false; + bool 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(1) = fact.R (); + retval(0) = fact.Q (); + } + 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(1) = fact.R (); + retval(0) = fact.Q (); + } + + } + 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) +*/ /* ;;; Local Variables: ***