Mercurial > octave-nkf
changeset 7700:efccca5f2ad7
more QR & Cholesky updating functions
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Mon, 07 Apr 2008 11:43:19 -0400 |
| parents | 27a5f578750c |
| children | b9f5f281b54b |
| files | libcruft/ChangeLog libcruft/qrupdate/Makefile.in libcruft/qrupdate/dchdex.f libcruft/qrupdate/dchinx.f libcruft/qrupdate/dqrqhu.f libcruft/qrupdate/dqrshc.f libcruft/qrupdate/zchdex.f libcruft/qrupdate/zchinx.f libcruft/qrupdate/zqrqhu.f libcruft/qrupdate/zqrshc.f liboctave/ChangeLog liboctave/CmplxCHOL.cc liboctave/CmplxCHOL.h liboctave/CmplxQR.cc liboctave/CmplxQR.h liboctave/dbleCHOL.cc liboctave/dbleCHOL.h liboctave/dbleQR.cc liboctave/dbleQR.h src/ChangeLog src/DLD-FUNCTIONS/chol.cc src/DLD-FUNCTIONS/qr.cc |
| diffstat | 22 files changed, 1368 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Fri Apr 04 15:57:31 2008 -0400 +++ b/libcruft/ChangeLog Mon Apr 07 11:43:19 2008 -0400 @@ -1,3 +1,10 @@ +2008-04-07 Jaroslav Hajek <highegg@gmail.com> + + * qrupdate/dqrqhu.f, qrupdate/zqrqhu.f, + * qrupdate/dqrshc.f, qrupdate/zqrshc.f, + * qrupdate/dchinx.f, qrupdate/zchinx.f, + * qrupdate/dchdex.f, qrupdate/zchdex.f: New files. + 2008-04-02 Jaroslav Hajek <highegg@gmail.com> * blas-xtra/xzdotu.f: Turn into simple wrapper for zdotu.
--- a/libcruft/qrupdate/Makefile.in Fri Apr 04 15:57:31 2008 -0400 +++ b/libcruft/qrupdate/Makefile.in Mon Apr 07 11:43:19 2008 -0400 @@ -31,8 +31,12 @@ dqrdec.f zqrdec.f \ dqrinr.f zqrinr.f \ dqrder.f zqrder.f \ + dqrshc.f zqrshc.f \ dch1up.f zch1up.f \ dch1dn.f zch1dn.f \ + dchinx.f zchinx.f \ + dchdex.f zchdex.f \ + dqrqhu.f zqrqhu.f \ dqrqhv.f zqrqhv.f \ dqhqr.f zqhqr.f
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dchdex.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,61 @@ +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 dchdex(n,R,R1,j) +c purpose: given an upper triangular matrix R that is a Cholesky +c factor of a symmetric positive definite matrix A, i.e. +c A = R'*R, this subroutine updates R -> R1 so that +c R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1]. +c (real version) +c arguments: +c n (in) the order of matrix R +c R (in) the original upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the deleted row/column + integer n,j,info + double precision R(n,n),R1(n-1,n-1) + double precision Qdum,c,s,rr + external dlacpy,dqhqr,dlartg + +c quick return if possible + if (n == 1) return + +c check arguments + info = 0 + if (n <= 0) then + info = 1 + else if (j < 1 .or. j > n) then + info = 4 + end if + if (info /= 0) then + call xerbla('DQRDEX',info) + end if + +c setup the new matrix R1 + if (j > 1) then + call dlacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1) + end if + if (j < n) then + call dlacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1) + call dqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1) +c eliminate R(n,n) + call dlartg(R1(n-1,n-1),R(n,n),c,s,rr) + R1(n-1,n-1) = rr + endif + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dchinx.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,108 @@ +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 dchinx(n,R,R1,j,u,info) +c purpose: given an upper triangular matrix R that is a Cholesky +c factor of a symmetric positive definite matrix A, i.e. +c A = R'*R, this subroutine updates R -> R1 so that +c R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u, +c jj = [1:j-1,j+1:n+1]. +c (real version) +c arguments: +c n (in) the order of matrix R +c R (in) the original upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the inserted row/column +c u (in) the vector (n+1) determining the rank-1 update +c info (out) on exit, if info = 1, the +c definiteness. + + integer n,j,info + double precision R(n,n),R1(n+1,n+1),u(n+1) + double precision rho,Qdum,w,dnrm2 + external dcopy,dlacpy,dtrsv,dnrm2 + integer jj + +c quick return if possible + if (n == 0) then + if (u(1) <= 0) then + info = 1 + return + else + R(1,1) = sqrt(u(1)) + end if + end if + +c check arguments + info = 0 + if (n < 0) then + info = 1 + else if (j < 1 .or. j > n+1) then + info = 4 + end if + if (info /= 0) then + call xerbla('DQRINX',info) + end if + +c copy shifted vector + if (j > 1) then + call dcopy(j-1,u,1,R1(1,j),1) + end if + w = u(j) + if (j < n+1) then + call dcopy(n-j+1,u(j+1),1,R1(j,j),1) + end if + +c check for singularity of R + do i = 1,n + if (R(i,i) == 0d0) then + info = 2 + return + end if + end do +c form R' \ u + call dtrsv('U','T','N',n,R,n,R1(1,j),1) + rho = dnrm2(n,R1(1,j),1) +c check positive definiteness + rho = u(j) - rho**2 + if (rho <= 0d0) then + info = 1 + return + end if + R1(n+1,n+1) = sqrt(rho) + +c setup the new matrix R1 + do i = 1,n+1 + R1(n+1,i) = 0d0 + end do + if (j > 1) then + call dlacpy('0',n,j-1,R(1,1),n,R1(1,1),n+1) + end if + if (j <= n) then + call dlacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1) +c retriangularize + jj = min(j+1,n) + call dqrqhu(0,n+1-j,n-j,Qdum,1,R1(j,jj),n+1,R1(j,j),w) + R1(j,j) = w + do jj = j+1,n + R1(jj,j) = 0d0 + end do + end if + + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqrqhu.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,78 @@ +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 dqrqhu(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 a k-vector u, +c this subroutine updates the matrices Q -> Q1 and +c R -> R1 so that Q1 = Q*G', R1 = G*R, u1(2:k) = 0 +c with G orthogonal, R1 upper Hessenberg, and u1 = G*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. +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 k-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 + external dlartg,drot + integer i,info +c quick return if possible. + if (k <= 0) 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('DQRQHU',info) + end if + rr = u(k) + do i = k-1,1,-1 + w = rr + if (w /= 0d0) then + call dlartg(u(i),w,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 if necessary + if (m > 0) then + call drot(m,Q(1,i),1,Q(1,i+1),1,c,s) + end if + else +c no rotation necessary + rr = u(i) + end if + end do + end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/dqrshc.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,97 @@ +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 dqrshc(m,n,k,Q,R,i,j) +c purpose: updates a QR factorization after circular shift of +c columns. +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 R -> R1 so that Q1 is again unitary, R1 upper trapezoidal, +c and +c Q1*R1 = A(:,p), where A = Q*R and p is the permutation +c [1:i-1,shift(i:j,-1),j+1:n] if i < j or +c [1:j-1,shift(j:i,+1),i+1:n] if j > i. +c if m == 0, the matrix Q is ignored. +c (real version) +c arguments: +c m (in) number of rows of the matrix Q, or 0 if Q is not needed. +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 R (io) on entry, the upper trapezoidal m-by-n matrix R. +c on exit, the updated matrix R1. +c i (in) the first index determining the range (see above) +c j (in) the second index determining the range (see above) +c + integer m,n,k,i,j + double precision Q(m,k),R(k,n) + external dswap,dqhqr + double precision w + integer l,jj,kk,info + +c quick return if possible + if (k <= 0 .or. n <= 1) return + info = 0 + if (m /= 0 .and. k > m) then + info = 3 + else if (i < 1 .or. i > n) then + info = 6 + else if (j < 1 .or. j > n) then + info = 7 + end if + if (info /= 0) then + call xerbla('DQRSHC',info) + end if + + if (i < j) then +c shift columns + do l = i,j-1 + call dswap(min(k,l+1),R(1,l),1,R(1,l+1),1) + end do +c retriangularize + if (i < k) then + kk = min(k,j) + if (m > 0) then + call dqhqr(m,n+1-i,kk+1-i,Q(1,i),m,R(i,i),k) + else + call dqhqr(0,n+1-i,kk+1-i,Q,1,R(i,i),k) + endif + end if + else if (j < i) then +c shift columns + do l = i,j+1,-1 + call dswap(min(k,i),R(1,l),1,R(1,l-1),1) + end do +c retriangularize + if (j < k) then + jj = min(j+1,n) + kk = min(k,i) + if (m > 0) then + call dqrqhu(m,n-j,kk+1-j,Q(1,j),m,R(j,jj),k,R(j,j),w) + else + call dqrqhu(0,n-j,kk+1-j,Q,1,R(j,jj),k,R(j,j),w) + end if + R(j,j) = w + do jj = j+1,kk + R(jj,j) = 0 + end do + end if + end if + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zchdex.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,62 @@ +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 zchdex(n,R,R1,j) +c purpose: given an upper triangular matrix R that is a Cholesky +c factor of a symmetric positive definite matrix A, i.e. +c A = R'*R, this subroutine updates R -> R1 so that +c R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1]. +c (real version) +c arguments: +c n (in) the order of matrix R +c R (in) the original upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the deleted row/column + integer n,j,info + double complex R(n,n),R1(n-1,n-1) + double precision c + double complex Qdum,s,rr + external zlacpy,zqhqr,zlartg + +c quick return if possible + if (n == 1) return + +c check arguments + info = 0 + if (n <= 0) then + info = 1 + else if (j < 1 .or. j > n) then + info = 4 + end if + if (info /= 0) then + call xerbla('ZQRDEX',info) + end if + +c setup the new matrix R1 + if (j > 1) then + call zlacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1) + end if + if (j < n) then + call zlacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1) + call zqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1) +c eliminate R(n,n) + call zlartg(R1(n-1,n-1),R(n,n),c,s,rr) + R1(n-1,n-1) = rr + endif + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zchinx.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,109 @@ +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 zchinx(n,R,R1,j,u,info) +c purpose: given an upper triangular matrix R that is a Cholesky +c factor of a symmetric positive definite matrix A, i.e. +c A = R'*R, this subroutine updates R -> R1 so that +c R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u, +c jj = [1:j-1,j+1:n+1]. +c (complex version) +c arguments: +c n (in) the order of matrix R +c R (in) the original upper trapezoidal matrix R +c R1 (out) the updated matrix R1 +c j (in) the position of the inserted row/column +c u (in) the vector (n+1) determining the rank-1 update +c info (out) on exit, if info = 1, the +c definiteness. + + integer n,j,info + double complex R(n,n),R1(n+1,n+1),u(n+1) + double precision rho,dznrm2 + double complex Qdum,w + external zcopy,zlacpy,ztrsv,dznrm2 + integer jj + +c quick return if possible + if (n == 0) then + if (real(u(1)) <= 0) then + info = 1 + return + else + R(1,1) = sqrt(real(u(1))) + end if + end if + +c check arguments + info = 0 + if (n < 0) then + info = 1 + else if (j < 1 .or. j > n+1) then + info = 4 + end if + if (info /= 0) then + call xerbla('ZQRINX',info) + end if + +c copy shifted vector + if (j > 1) then + call zcopy(j-1,u,1,R1(1,j),1) + end if + w = u(j) + if (j < n+1) then + call zcopy(n-j+1,u(j+1),1,R1(j,j),1) + end if + +c check for singularity of R + do i = 1,n + if (R(i,i) == 0d0) then + info = 2 + return + end if + end do +c form R' \ u + call ztrsv('U','T','N',n,R,n,R1(1,j),1) + rho = dznrm2(n,R1(1,j),1) +c check positive definiteness + rho = u(j) - rho**2 + if (rho <= 0d0) then + info = 1 + return + end if + R1(n+1,n+1) = sqrt(rho) + +c setup the new matrix R1 + do i = 1,n+1 + R1(n+1,i) = 0d0 + end do + if (j > 1) then + call zlacpy('0',j-1,n,R(1,1),n,R1(1,1),n+1) + end if + if (j <= n) then + call zlacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1) +c retriangularize + jj = min(j+1,n) + call zqrqhu(0,n+1-j,n-j,Qdum,1,R1(j,jj),n+1,R1(j,j),w) + R1(j,j) = w + do jj = j+1,n + R1(jj,j) = 0d0 + end do + end if + + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqrqhu.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,78 @@ +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 zqrqhu(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 a k-vector u, +c this subroutine updates the matrices Q -> Q1 and +c R -> R1 so that Q1 = Q*G', R1 = G*R, u1(2:k) = 0 +c with G unitary, R1 upper Hessenberg, and u1 = G*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. +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 k-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 + external zlartg,zrot + integer i,info +c quick return if possible. + if (k <= 0) 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('ZQRQHU',info) + end if + rr = u(k) + do i = k-1,1,-1 + w = rr + if (w /= dcmplx(0d0,0d0)) then + call zlartg(u(i),w,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 if necessary + if (m > 0) then + call zrot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s)) + end if + else +c no rotation necessary + rr = u(i) + end if + end do + end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/qrupdate/zqrshc.f Mon Apr 07 11:43:19 2008 -0400 @@ -0,0 +1,97 @@ +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 zqrshc(m,n,k,Q,R,i,j) +c purpose: updates a QR factorization after circular shift of +c columns. +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 R -> R1 so that Q1 is again unitary, R1 upper trapezoidal, +c and +c Q1*R1 = A(:,p), where A = Q*R and p is the permutation +c [1:i-1,shift(i:j,-1),j+1:n] if i < j or +c [1:j-1,shift(j:i,+1),i+1:n] if j > i. +c if m == 0, the matrix Q is ignored. +c (real version) +c arguments: +c m (in) number of rows of the matrix Q, or 0 if Q is not needed. +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 R (io) on entry, the upper trapezoidal m-by-n matrix R. +c on exit, the updated matrix R1. +c i (in) the first index determining the range (see above) +c j (in) the second index determining the range (see above) +c + integer m,n,k,i,j + double complex Q(m,k),R(k,n) + external zswap,zqhqr + double complex w + integer l,jj,kk,info + +c quick return if possible + if (k <= 0 .or. n <= 1) return + info = 0 + if (m /= 0 .and. k > m) then + info = 3 + else if (i < 1 .or. i > n) then + info = 6 + else if (j < 1 .or. j > n) then + info = 7 + end if + if (info /= 0) then + call xerbla('ZQRSHC',info) + end if + + if (i < j) then +c shift columns + do l = i,j-1 + call zswap(min(k,l+1),R(1,l),1,R(1,l+1),1) + end do +c retriangularize + if (i < k) then + kk = min(k,j) + if (m > 0) then + call zqhqr(m,n+1-i,kk+1-i,Q(1,i),m,R(i,i),k) + else + call zqhqr(0,n+1-i,kk+1-i,Q,1,R(i,i),k) + endif + end if + else if (j < i) then +c shift columns + do l = i,j+1,-1 + call zswap(min(k,i),R(1,l),1,R(1,l-1),1) + end do +c retriangularize + if (j < k) then + jj = min(j+1,n) + kk = min(k,i) + if (m > 0) then + call zqrqhu(m,n-j,kk+1-j,Q(1,j),m,R(j,jj),k,R(j,j),w) + else + call zqrqhu(0,n-j,kk+1-j,Q,1,R(j,jj),k,R(j,j),w) + end if + R(j,j) = w + do jj = j+1,kk + R(jj,j) = 0 + end do + end if + end if + end
--- a/liboctave/ChangeLog Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/ChangeLog Mon Apr 07 11:43:19 2008 -0400 @@ -1,3 +1,12 @@ +2008-04-07 Jaroslav Hajek <highegg@gmail.com> + + * dbleQR.h, dbleQR.cc (QR::shift_cols): New method. + * CmplxQR.h, CmplxQR.cc (ComplexQR::shift_cols): New method. + * dbleCHOL.h, dbleCHOL.cc (CHOL::insert_sym, CHOL::delete_sym, + CHOL::shift_sym): New methods. + * CmplxCHOL.h, CmplxCHOL.cc (ComplexCHOL::insert_sym, + ComplexCHOL::delete_sym, ComplexCHOL::shift_sym): New methods. + 2008-04-03 John W. Eaton <jwe@octave.org> * lo-sysdep.cc [__WIN32__ && ! __CYGWIN__]: Include windows.h.
--- a/liboctave/CmplxCHOL.cc Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/CmplxCHOL.cc Mon Apr 07 11:43:19 2008 -0400 @@ -57,6 +57,17 @@ F77_RET_T F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, Complex*, double*, octave_idx_type&); + + F77_RET_T + F77_FUNC (zqrshc, ZQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + Complex*, Complex*, const octave_idx_type&, const octave_idx_type&); + + F77_RET_T + F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, const Complex*, Complex*, const octave_idx_type&, + const Complex*, octave_idx_type&); + + F77_RET_T + F77_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, const Complex*, Complex*, const octave_idx_type&); } octave_idx_type @@ -210,6 +221,59 @@ return info; } +octave_idx_type +ComplexCHOL::insert_sym (const ComplexMatrix& u, octave_idx_type j) +{ + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () != n+1) + (*current_liboctave_error_handler) ("CHOL insert dimension mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("CHOL insert index out of range"); + else + { + ComplexMatrix chol_mat1 (n+1, n+1); + + F77_XFCN (zchinx, ZCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), + j+1, u.data (), info)); + + chol_mat = chol_mat1; + } + + return info; +} + +void +ComplexCHOL::delete_sym (octave_idx_type j) +{ + octave_idx_type n = chol_mat.rows (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("CHOL insert index out of range"); + else + { + ComplexMatrix chol_mat1 (n-1, n-1); + + F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + + chol_mat = chol_mat1; + } +} + +void +ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j) +{ + octave_idx_type n = chol_mat.rows (); + Complex dummy; + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("CHOL shift index out of range"); + else + F77_XFCN (zqrshc, ZQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); +} + ComplexMatrix chol2inv (const ComplexMatrix& r) {
--- a/liboctave/CmplxCHOL.h Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/CmplxCHOL.h Mon Apr 07 11:43:19 2008 -0400 @@ -71,6 +71,12 @@ octave_idx_type downdate (const ComplexMatrix& u); + octave_idx_type insert_sym (const ComplexMatrix& u, octave_idx_type j); + + void delete_sym (octave_idx_type j); + + void shift_sym (octave_idx_type i, octave_idx_type j); + friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexCHOL& a); private:
--- a/liboctave/CmplxQR.cc Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/CmplxQR.cc Mon Apr 07 11:43:19 2008 -0400 @@ -68,6 +68,10 @@ F77_FUNC (zqrder, ZQRDER) (const octave_idx_type&, const octave_idx_type&, const Complex*, Complex*, const Complex*, Complex *, const octave_idx_type&); + + F77_RET_T + F77_FUNC (zqrshc, ZQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + Complex*, Complex*, const octave_idx_type&, const octave_idx_type&); } ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type) @@ -272,6 +276,19 @@ } void +ComplexQR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("QR shift index out of range"); + else + F77_XFCN (zqrshc, ZQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); +} + +void ComplexQR::economize (void) { octave_idx_type r_nc = r.columns ();
--- a/liboctave/CmplxQR.h Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/CmplxQR.h Mon Apr 07 11:43:19 2008 -0400 @@ -76,6 +76,8 @@ void delete_row (octave_idx_type j); + void shift_cols (octave_idx_type i, octave_idx_type j); + void economize(); friend std::ostream& operator << (std::ostream&, const ComplexQR&);
--- a/liboctave/dbleCHOL.cc Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/dbleCHOL.cc Mon Apr 07 11:43:19 2008 -0400 @@ -57,6 +57,17 @@ F77_RET_T F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, octave_idx_type&); + + F77_RET_T + F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + double*, double*, const octave_idx_type&, const octave_idx_type&); + + F77_RET_T + F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, const double*, double*, const octave_idx_type&, + const double*, octave_idx_type&); + + F77_RET_T + F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, const double*, double*, const octave_idx_type&); } octave_idx_type @@ -214,6 +225,59 @@ return info; } +octave_idx_type +CHOL::insert_sym (const Matrix& u, octave_idx_type j) +{ + octave_idx_type info = -1; + + octave_idx_type n = chol_mat.rows (); + + if (u.length () != n+1) + (*current_liboctave_error_handler) ("CHOL insert dimension mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("CHOL insert index out of range"); + else + { + Matrix chol_mat1 (n+1, n+1); + + F77_XFCN (dchinx, DCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), + j+1, u.data (), info)); + + chol_mat = chol_mat1; + } + + return info; +} + +void +CHOL::delete_sym (octave_idx_type j) +{ + octave_idx_type n = chol_mat.rows (); + + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("CHOL insert index out of range"); + else + { + Matrix chol_mat1 (n-1, n-1); + + F77_XFCN (dchdex, DCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + + chol_mat = chol_mat1; + } +} + +void +CHOL::shift_sym (octave_idx_type i, octave_idx_type j) +{ + octave_idx_type n = chol_mat.rows (); + double dummy; + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("CHOL shift index out of range"); + else + F77_XFCN (dqrshc, DQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); +} + Matrix chol2inv (const Matrix& r) {
--- a/liboctave/dbleCHOL.h Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/dbleCHOL.h Mon Apr 07 11:43:19 2008 -0400 @@ -68,6 +68,12 @@ octave_idx_type downdate (const Matrix& u); + octave_idx_type insert_sym (const Matrix& u, octave_idx_type j); + + void delete_sym (octave_idx_type j); + + void shift_sym (octave_idx_type i, octave_idx_type j); + friend OCTAVE_API std::ostream& operator << (std::ostream& os, const CHOL& a); private:
--- a/liboctave/dbleQR.cc Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/dbleQR.cc Mon Apr 07 11:43:19 2008 -0400 @@ -64,6 +64,10 @@ F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&, const double*, double*, const double*, double *, const octave_idx_type&); + + F77_RET_T + F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + double*, double*, const octave_idx_type&, const octave_idx_type&); } QR::QR (const Matrix& a, QR::type qr_type) @@ -261,6 +265,19 @@ } void +QR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("QR shift index out of range"); + else + F77_XFCN (dqrshc, DQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); +} + +void QR::economize (void) { octave_idx_type r_nc = r.columns ();
--- a/liboctave/dbleQR.h Fri Apr 04 15:57:31 2008 -0400 +++ b/liboctave/dbleQR.h Mon Apr 07 11:43:19 2008 -0400 @@ -81,6 +81,8 @@ void delete_row (octave_idx_type j); + void shift_cols (octave_idx_type i, octave_idx_type j); + void economize (void); friend std::ostream& operator << (std::ostream&, const QR&);
--- a/src/ChangeLog Fri Apr 04 15:57:31 2008 -0400 +++ b/src/ChangeLog Mon Apr 07 11:43:19 2008 -0400 @@ -1,3 +1,9 @@ +2008-04-07 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/qr.cc (Fqrshift): New function. + * DLD-FUNCTIONS/chol.cc (Fcholinsert, Fcholdelete, Fcholshift): + New functions. + 2008-04-04 John W. Eaton <jwe@octave.org> * parse.y (make_constant): Handle escape sequences in dq-strings.
--- a/src/DLD-FUNCTIONS/chol.cc Fri Apr 04 15:57:31 2008 -0400 +++ b/src/DLD-FUNCTIONS/chol.cc Mon Apr 07 11:43:19 2008 -0400 @@ -21,9 +21,9 @@ */ -// The cholupdate function was written by Jaroslav Hajek -// <highegg@gmail.com>, Copyright (C) 2008 VZLU Prague, a.s., Czech -// Republic. +// The cholupdate, cholinsert, choldelete and cholshift 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> @@ -600,6 +600,350 @@ %! assert(norm(R1 - R,Inf) < 1e1*eps) */ +DEFUN_DLD (cholinsert, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})\n\ +Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\ +positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\ +return the QR@tie{}factorization of\n\ +@var{A1}, where @w{A1(p,p) = A}, @w{A1(:,j) = A1(j,:)' = u} and\n\ +@w{p = [1:j-1,j+1:n+1]}. @w{u(j)} should be positive.\n\ +On return, @var{info} is set to\n\ +@itemize\n\ +@item 0 if the insertion was successful,\n\ +@item 1 if @var{A1} is not positive definite,\n\ +@item 2 if @var{R} is singular.\n\ +@end itemize\n\ +\n\ +If @var{info} is not present, an error message is printed in cases 1 and 2.\n\ +@seealso{chol, cholupdate, choldelete}\n\ +@end deftypefn") +{ + octave_idx_type nargin = args.length (); + + octave_value_list retval; + + if (nargin != 3) + { + print_usage (); + return retval; + } + + octave_value argr = args(0); + octave_value argj = args(1); + octave_value argu = args(2); + + if (argr.is_matrix_type () && argu.is_matrix_type () + && argj.is_real_scalar ()) + { + octave_idx_type n = argr.rows (); + octave_idx_type j = argj.scalar_value (); + + if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1) + { + if (j > 0 && j <= n+1) + { + if (argr.is_real_matrix () && argu.is_real_matrix ()) + { + // real case + Matrix R = argr.matrix_value (); + Matrix u = argu.matrix_value (); + + CHOL fact; + fact.set (R); + int err = fact.insert_sym (u, j-1); + + if (nargout > 1) + retval(1) = err; + else if (err) + error ("cholinsert: insertion violates positiveness"); + + retval(0) = fact.chol_matrix (); + } + else + { + // complex case + ComplexMatrix R = argr.complex_matrix_value (); + ComplexMatrix u = argu.complex_matrix_value (); + + ComplexCHOL fact; + fact.set (R); + int err = fact.insert_sym (u, j-1); + + if (nargout > 1) + retval(1) = err; + else if (err) + error ("cholinsert: insertion violates positiveness"); + + retval(0) = fact.chol_matrix (); + } + } + else + error ("cholinsert: index out of range"); + } + else + error ("cholinsert: dimension mismatch"); + } + else + print_usage (); + + return retval; +} + +/* +%!test +%! A = [ 0.436997 -0.131721 0.124120 -0.061673 ; +%! -0.131721 0.738529 0.019851 -0.140295 ; +%! 0.124120 0.019851 0.354879 -0.059472 ; +%! -0.061673 -0.140295 -0.059472 0.600939 ]; +%! +%! u = [ 0.35080 ; +%! 0.63930 ; +%! 3.31057 ; +%! -0.13825 ; +%! 0.45266 ]; +%! +%! R = chol(A); +%! +%! j = 3; p = [1:j-1, j+1:5]; +%! R1 = cholinsert(R,j,u); A1 = R1'*R1; +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(A1(p,p) - A,Inf) < 1e1*eps) +%! +%!test +%! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; +%! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ; +%! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ; +%! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ]; +%! +%! u = [ 0.35080 + 0.04298i; +%! 0.63930 + 0.23778i; +%! 3.31057 + 0.00000i; +%! -0.13825 + 0.19879i; +%! 0.45266 + 0.50020i]; +%! +%! R = chol(A); +%! +%! j = 3; p = [1:j-1, j+1:5]; +%! R1 = cholinsert(R,j,u); A1 = R1'*R1; +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(A1(p,p) - A,Inf) < 1e1*eps) +%! +*/ + +DEFUN_DLD (choldelete, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{R1} =} choldelete (@var{R}, @var{j})\n\ +Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\ +positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\ +return the QR@tie{}factorization of @w{A(p,p)}, where @w{p = [1:j-1,j+1:n+1]}.\n\ +@seealso{chol, cholupdate, cholinsert}\n\ +@end deftypefn") +{ + octave_idx_type nargin = args.length (); + + octave_value_list retval; + + if (nargin != 2) + { + print_usage (); + return retval; + } + + octave_value argr = args(0); + octave_value argj = args(1); + + if (argr.is_matrix_type () && argj.is_real_scalar ()) + { + octave_idx_type n = argr.rows (); + octave_idx_type j = argj.scalar_value (); + + if (argr.columns () == n) + { + if (j > 0 && j <= n) + { + if (argr.is_real_matrix ()) + { + // real case + Matrix R = argr.matrix_value (); + + CHOL fact; + fact.set (R); + fact.delete_sym (j-1); + + retval(0) = fact.chol_matrix (); + } + else + { + // complex case + ComplexMatrix R = argr.complex_matrix_value (); + + ComplexCHOL fact; + fact.set (R); + fact.delete_sym (j-1); + + retval(0) = fact.chol_matrix (); + } + } + else + error ("choldelete: index out of range"); + } + else + error ("choldelete: dimension mismatch"); + } + else + print_usage (); + + return retval; +} + +/* +%!test +%! A = [ 0.436997 -0.131721 0.124120 -0.061673 ; +%! -0.131721 0.738529 0.019851 -0.140295 ; +%! 0.124120 0.019851 0.354879 -0.059472 ; +%! -0.061673 -0.140295 -0.059472 0.600939 ]; +%! +%! R = chol(A); +%! +%! j = 3; p = [1:j-1,j+1:4]; +%! R1 = choldelete(R,j); +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) +%! +%!test +%! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; +%! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ; +%! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ; +%! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ]; +%! +%! R = chol(A); +%! +%! j = 3; p = [1:j-1,j+1:4]; +%! R1 = choldelete(R,j); +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) +*/ + +DEFUN_DLD (cholshift, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})\n\ +Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\ +positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\ +return the QR@tie{}factorization of\n\ +@w{@var{A}(p,p)}, where @w{p} is the permutation @*\n\ +@code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*\n\ + or @*\n\ +@code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @*\n\ +\n\ +@seealso{chol, cholinsert, choldelete}\n\ +@end deftypefn") +{ + octave_idx_type nargin = args.length (); + + octave_value_list retval; + + if (nargin != 3) + { + print_usage (); + return retval; + } + + octave_value argr = args(0); + octave_value argi = args(1); + octave_value argj = args(2); + + if (argr.is_matrix_type () && argi.is_real_scalar () && argj.is_real_scalar ()) + { + octave_idx_type n = argr.rows (); + octave_idx_type i = argi.scalar_value (); + octave_idx_type j = argj.scalar_value (); + + if (argr.columns () == n) + { + if (j > 0 && j <= n+1 && i > 0 && i <= n+1) + { + if (argr.is_real_matrix ()) + { + // real case + Matrix R = argr.matrix_value (); + + CHOL fact; + fact.set (R); + fact.shift_sym (i-1, j-1); + + retval(0) = fact.chol_matrix (); + } + else + { + // complex case + ComplexMatrix R = argr.complex_matrix_value (); + + ComplexCHOL fact; + fact.set (R); + fact.shift_sym (i-1, j-1); + + retval(0) = fact.chol_matrix (); + } + } + else + error ("cholshift: index out of range"); + } + else + error ("cholshift: dimension mismatch"); + } + else + print_usage (); + + return retval; +} + +/* +%!test +%! A = [ 0.436997 -0.131721 0.124120 -0.061673 ; +%! -0.131721 0.738529 0.019851 -0.140295 ; +%! 0.124120 0.019851 0.354879 -0.059472 ; +%! -0.061673 -0.140295 -0.059472 0.600939 ]; +%! +%! R = chol(A); +%! +%! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; +%! R1 = cholshift(R,i,j); +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) +%! +%! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4]; +%! R1 = cholshift(R,i,j); +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) +%! +%!test +%! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; +%! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ; +%! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ; +%! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ]; +%! +%! R = chol(A); +%! +%! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; +%! R1 = cholshift(R,i,j); +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) +%! +%! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4]; +%! R1 = cholshift(R,i,j); +%! +%! assert(norm(triu(R1)-R1,Inf) == 0) +%! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) +*/ + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/src/DLD-FUNCTIONS/qr.cc Fri Apr 04 15:57:31 2008 -0400 +++ b/src/DLD-FUNCTIONS/qr.cc Mon Apr 07 11:43:19 2008 -0400 @@ -20,7 +20,7 @@ */ -// The qrupdate, qrinsert, and qrdelete functions were written by +// The qrupdate, qrinsert, qrdelete and qrshift functions were written by // Jaroslav Hajek <highegg@gmail.com>, Copyright (C) 2008 VZLU // Prague, a.s., Czech Republic. @@ -913,6 +913,132 @@ %! assert(norm(vec(Q*R - [A(1:2,:);A(4:5,:)]),Inf) < norm(A)*1e1*eps) */ +DEFUN_DLD (qrshift, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}] =} qrshift (@var{Q}, @var{R}, @var{i}, @var{j})\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}(:,p)}, where @w{p} is the permutation @*\n\ +@code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*\n\ + or @*\n\ +@code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @*\n\ +\n\ +@seealso{qr, qrinsert, qrdelete}\n\ +@end deftypefn") +{ + octave_idx_type nargin = args.length (); + octave_value_list retval; + + if (nargin != 4) + { + print_usage (); + return retval; + } + + octave_value argq = args(0); + octave_value argr = args(1); + octave_value argi = args(2); + octave_value argj = args(3); + + if (argq.is_matrix_type () && argr.is_matrix_type () + && argi.is_real_scalar () && argj.is_real_scalar ()) + { + octave_idx_type m = argq.rows (); + octave_idx_type n = argr.columns (); + octave_idx_type k = argq.columns (); + + if (argr.rows () == k) + { + octave_idx_type i = argi.scalar_value (); + octave_idx_type j = argj.scalar_value (); + if (i > 1 && i <= n && j > 1 && j <= n) + { + if (argq.is_real_matrix () + && argr.is_real_matrix ()) + { + // all real case + Matrix Q = argq.matrix_value (); + Matrix R = argr.matrix_value (); + + QR fact (Q, R); + fact.shift_cols (i-1, j-1); + + 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); + fact.shift_cols (i-1, j-1); + + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + } + else + error ("qrshift: index out of range"); + } + else + error ("qrshift: 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 ].'; +%! +%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrshift(Q,R,i,j); +%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - A(:,p)),Inf) < norm(A)*1e1*eps) +%! +%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrshift(Q,R,i,j); +%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - A(:,p)),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 ].'; +%! +%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrshift(Q,R,i,j); +%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - A(:,p)),Inf) < norm(A)*1e1*eps) +%! +%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5]; +%! +%! [Q,R] = qr(A); +%! [Q,R] = qrshift(Q,R,i,j); +%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps) +%! assert(norm(vec(triu(R)-R),Inf) == 0) +%! assert(norm(vec(Q*R - A(:,p)),Inf) < norm(A)*1e1*eps) +*/ + /* ;;; Local Variables: *** ;;; mode: C++ ***