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