# HG changeset patch # User Jaroslav Hajek # Date 1232482602 -3600 # Node ID d66c9b6e506ab48a61da39b04218d01433816364 # Parent 3d8a914c580e99d0294568ce00e0bbdb8f52fc88 imported patch qrupdate.diff diff -r 3d8a914c580e -r d66c9b6e506a Makeconf.in --- a/Makeconf.in Tue Jan 20 21:15:17 2009 +0100 +++ b/Makeconf.in Tue Jan 20 21:16:42 2009 +0100 @@ -239,6 +239,7 @@ CHOLMOD_LIBS = @CHOLMOD_LIBS@ CXSPARSE_LIBS = @CXSPARSE_LIBS@ OPENGL_LIBS = @OPENGL_LIBS@ +QRUPDATE_LIBS = @QRUPDATE_LIBS@ ARPACK_LIBS = @ARPACK_LIBS@ LIBS = @LIBS@ diff -r 3d8a914c580e -r d66c9b6e506a configure.in --- a/configure.in Tue Jan 20 21:15:17 2009 +0100 +++ b/configure.in Tue Jan 20 21:16:42 2009 +0100 @@ -865,6 +865,34 @@ AC_MSG_WARN($warn_blas_f77_incompatible) fi +QRUPDATE_LIBS= +AC_SUBST(QRUPDATE_LIBS) + +# Check for the qrupdate library +AC_ARG_WITH(qrupdate, + [AS_HELP_STRING([--without-qrupdate], + [don't use qrupdate, disable QR & Cholesky updating functions])], + with_qrupdate=$withval, with_qrupdate=yes) + +warn_qrupdate="qrupdate not found. The QR & Cholesky updating function will not be available." +if test "$with_qrupdate" = yes; then + with_qrupdate=no + if $have_fortran_compiler; then + AC_F77_FUNC(sqr1up) + elif $have_f2c; then + sqr1up=sqr1up_ + fi + AC_CHECK_LIB(qrupdate, $sqr1up, + [QRUPDATE_LIBS="-lqrupdate"; with_qrupdate=yes], [], [$BLAS_LIBS $FLIBS]) + if test "$with_qrupdate" = yes; then + AC_DEFINE(HAVE_QRUPDATE, 1, [Define if the qrupdate library is used.]) + warn_qrupdate= + fi +fi +if test -n "$warn_qrupdate"; then + AC_MSG_WARN($warn_qrupdate) +fi + # Check for AMD library AMD_LIBS= AC_SUBST(AMD_LIBS) @@ -1997,7 +2025,6 @@ libcruft/lapack/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]) @@ -2034,6 +2061,7 @@ CHOLMOD libraries: $CHOLMOD_LIBS CXSPARSE libraries: $CXSPARSE_LIBS ARPACK libraries: $ARPACK_LIBS + QRUPDATE libraries: $QRUPDATE_LIBS HDF5 libraries: $HDF5_LIBS CURL libraries: $CURL_LIBS REGEX libraries: $REGEX_LIBS @@ -2120,6 +2148,11 @@ warn_msg_printed=true fi +if test -n "$warn_qrupdate"; then + AC_MSG_WARN($warn_qrupdate) + warn_msg_printed=true +fi + if test -n "$warn_amd"; then AC_MSG_WARN($warn_amd) warn_msg_printed=true diff -r 3d8a914c580e -r d66c9b6e506a libcruft/Makefile.in --- a/libcruft/Makefile.in Tue Jan 20 21:15:17 2009 +0100 +++ b/libcruft/Makefile.in Tue Jan 20 21:16:42 2009 +0100 @@ -43,7 +43,7 @@ CRUFT_DIRS = amos @BLAS_DIR@ blas-xtra daspk dasrt dassl \ @FFT_DIR@ @LAPACK_DIR@ lapack-xtra \ - misc odepack ordered-qz quadpack qrupdate ranlib \ + misc odepack ordered-qz quadpack ranlib \ slatec-err slatec-fn villad SUBDIRS = $(CRUFT_DIRS) diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/Makefile.in --- a/libcruft/qrupdate/Makefile.in Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,58 +0,0 @@ -# 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 -# . - -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 \ - 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 \ - sch1up.f cch1up.f \ - sqrinc.f cqrinc.f \ - sqrdec.f cqrdec.f \ - sqrinr.f cqrinr.f \ - sqrder.f cqrder.f \ - sqrshc.f cqrshc.f \ - sqr1up.f cqr1up.f \ - sch1dn.f cch1dn.f \ - schinx.f cchinx.f \ - schdex.f cchdex.f \ - sqrqhu.f cqrqhu.f \ - sqrqhv.f cqrqhv.f \ - sqhqr.f cqhqr.f - -include $(TOPDIR)/Makeconf - -include ../Makerules diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cch1dn.f --- a/libcruft/qrupdate/cch1dn.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -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 cch1dn(n,R,u,w,info) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a hermitian positive definite matrix A, i.e. -c A = R'*R, this subroutine downdates R -> R1 so that -c R1'*R1 = A - u*u' -c (complex version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n,info - complex R(n,n),u(n) - real w(n) - external ctrsv,clartg,scnrm2 - real rho,scnrm2 - complex crho,rr,ui,t - integer i,j - -c quick return if possible - if (n <= 0) return -c check for singularity of R - do i = 1,n - if (R(i,i) == 0e0) then - info = 2 - return - end if - end do -c form R' \ u - call ctrsv('U','C','N',n,R,n,u,1) - rho = scnrm2(n,u,1) -c check positive definiteness - rho = 1 - rho**2 - if (rho <= 0e0) then - info = 1 - return - end if - crho = sqrt(rho) -c eliminate R' \ u - do i = n,1,-1 - ui = u(i) -c generate next rotation - call clartg(crho,ui,w(i),u(i),rr) - crho = rr - end do -c apply rotations - do i = n,1,-1 - ui = 0e0 - do j = i,1,-1 - t = w(j)*ui + u(j)*R(j,i) - R(j,i) = w(j)*R(j,i) - conjg(u(j))*ui - ui = t - end do - end do - - info = 0 - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cch1up.f --- a/libcruft/qrupdate/cch1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,56 +0,0 @@ -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 cch1up(n,R,u,w) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a hermitian positive definite matrix A, i.e. -c A = R'*R, this subroutine updates R -> R1 so that -c R1'*R1 = A + u*u' or A - u*u' -c (complex version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a real workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n - complex R(n,n),u(n) - real w(n) - external clartg - complex rr,ui,t - integer i,j - - do i = 1,n -c apply stored rotations, column-wise - ui = conjg(u(i)) - do j = 1,i-1 - t = w(j)*R(j,i) + u(j)*ui - ui = w(j)*ui - conjg(u(j))*R(j,i) - R(j,i) = t - end do -c generate next rotation - call clartg(R(i,i),ui,w(i),u(i),rr) - R(i,i) = rr - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cchdex.f --- a/libcruft/qrupdate/cchdex.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,62 +0,0 @@ -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 cchdex(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 (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 deleted row/column - integer n,j,info - complex R(n,n),R1(n-1,n-1) - real c - complex Qdum,s,rr - external xerbla,clacpy,cqhqr,clartg - -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('CCHDEX',info) - end if - -c setup the new matrix R1 - if (j > 1) then - call clacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1) - end if - if (j < n) then - call clacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1) - call cqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1) -c eliminate R(n,n) - call clartg(R1(n-1,n-1),R(n,n),c,s,rr) - R1(n-1,n-1) = rr - endif - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cchinx.f --- a/libcruft/qrupdate/cchinx.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,109 +0,0 @@ -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 cchinx(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 - complex R(n,n),R1(n+1,n+1),u(n+1) - real rho,scnrm2 - complex Qdum,w - external xerbla,ccopy,clacpy,ctrsv,scnrm2,cqrqhu - 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('CCHINX',info) - end if - -c copy shifted vector - if (j > 1) then - call ccopy(j-1,u,1,R1(1,j),1) - end if - w = u(j) - if (j < n+1) then - call ccopy(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) == 0e0) then - info = 2 - return - end if - end do -c form R' \ u - call ctrsv('U','T','N',n,R,n,R1(1,j),1) - rho = scnrm2(n,R1(1,j),1) -c check positive definiteness - rho = u(j) - rho**2 - if (rho <= 0e0) 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) = 0e0 - end do - if (j > 1) then - call clacpy('0',j-1,n,R(1,1),n,R1(1,1),n+1) - end if - if (j <= n) then - call clacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1) -c retriangularize - jj = min(j+1,n) - call cqrqhu(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) = 0e0 - end do - end if - - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqhqr.f --- a/libcruft/qrupdate/cqhqr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,69 +0,0 @@ -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 cqhqr(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 - complex Q(ldq,*),R(ldr,*) - real c - complex s,rr - external xerbla,clartg,crot - 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('CQHQR',info) - end if -c triangularize - do i = 1,min(k-1,n) - call clartg(R(i,i),R(i+1,i),c,s,rr) - R(i,i) = rr - R(i+1,i) = 0e0 - if (i < n) then - call crot(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 crot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s)) - end if - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqr1up.f --- a/libcruft/qrupdate/cqr1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -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 cqr1up(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 - complex Q(m,k),R(k,n),u(m),v(n) - complex w - external cqrqhv,cqhqr - integer i -c quick return if possible - if (m <= 0 .or. n <= 0) return -c eliminate tail of Q'*u - call cqrqhv(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 cqhqr(m,n,k,Q,m,R,k) - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqrdec.f --- a/libcruft/qrupdate/cqrdec.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -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 cqrdec(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 - complex Q(m,k),R(k,n),R1(k,n-1) - external xerbla,ccopy,cqhqr - 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('CQRDEC',info) - end if -c copy leading portion - call ccopy(k*(j-1),R,1,R1,1) - if (j < n) then -c copy trailing portion of R - call ccopy(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 cqhqr(m,n-j,k-j+1,Q(1,j),m,R1(j,j),k) - end if - end if - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqrder.f --- a/libcruft/qrupdate/cqrder.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -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 cqrder(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 - complex Q(m,m),Q1(m-1,m-1),R(m,n),R1(m-1,n) - real c - complex s,rr,w - external xerbla,clacpy,clartg,crot,csscal,caxpy - 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('CQRDER',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 clacpy('0',j-1,m-1,Q(1,2),m,Q1(1,1),m-1) - end if - if (j < m) then - call clacpy('0',m-j,m-1,Q(j+1,2),m,Q1(j,1),m-1) - end if -c setup the new matrix R1 - call clacpy('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 clartg(Q(j,i),w,c,s,rr) - w = rr -c apply rotation to rows of R1 - if (i <= n) then - call crot(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 crot(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 clartg(Q(j,1),w,c,s,rr) - w = rr - call csscal(n,c,R1(1,1),m-1) - call caxpy(n,-s,R(1,1),m,R1(1,1),m-1) -c apply rotation to columns of Q1 - call csscal(m-1,c,Q1(1,1),1) - if (j > 1) then - call caxpy(j-1,-conjg(s),Q(1,1),1,Q1(1,1),1) - end if - if (j < m) then - call caxpy(m-j,-conjg(s),Q(j+1,1),1,Q1(j,1),1) - end if - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqrinc.f --- a/libcruft/qrupdate/cqrinc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,74 +0,0 @@ -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 cqrinc(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 - complex Q(m,k),R(k,n),R1(k,n+1),x(m) - complex w - external xerbla,ccopy,cqrqhv,cgemv - 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('CQRINC',info) - end if -c copy leading portion of R - call ccopy(k*(j-1),R,1,R1,1) - if (j <= n) then - call ccopy(k*(n+1-j),R(1,j),1,R1(1,j+1),1) - end if - call cgemv('C',m,min(k,j-1),cmplx(1e0),Q,m,x,1, - + cmplx(0e0),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 cqrqhv(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) = 0e0 - end do - end if - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqrinr.f --- a/libcruft/qrupdate/cqrinr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -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 cqrinr(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 - complex Q(m,m),Q1(m+1,m+1),R(m,n),R1(m+1,n),x(n) - external xerbla,clacpy,ccopy,cqhqr - 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('CQRINR',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 clacpy('0',j-1,m,Q(1,1),m,Q1(1,2),m+1) - end if - if (j <= m) then - call clacpy('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) = 0e0 - Q1(j,i) = 0e0 - end do - Q1(j,1) = 1e0 -c setup the new matrix R1 - call ccopy(n,x,1,R1(1,1),m+1) - call clacpy('0',m,n,R(1,1),m,R1(2,1),m+1) -c rotate to form proper QR - call cqhqr(m+1,n,m+1,Q1,m+1,R1,m+1) - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqrqhu.f --- a/libcruft/qrupdate/cqrqhu.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -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 cqrqhu(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 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 u (in) the k-vector u. -c rr (out) the first element of Q1'*u on exit. -c -c if Q is unitary, so is Q1. It is not strictly -c necessary, however. - integer m,n,k,ldq,ldr - complex Q(ldq,*),R(ldr,*),u(*),rr - real c - complex s,w - external xerbla,clartg,crot - 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('CQRQHU',info) - end if - rr = u(k) - do i = k-1,1,-1 - w = rr - if (w /= cmplx(0e0,0e0)) then - call clartg(u(i),w,c,s,rr) -c apply rotation to rows of R if necessary - if (i <= n) then - call crot(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 crot(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 - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqrqhv.f --- a/libcruft/qrupdate/cqrqhv.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -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 cqrqhv(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 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 u (in) the m-vector u. -c rr (out) the first element of Q1'*u on exit. -c -c if Q is unitary, so is Q1. It is not strictly -c necessary, however. - integer m,n,k,ldq,ldr - complex Q(ldq,*),R(ldr,*),u(*),rr - real c - complex s,w,w1,cdotc - external xerbla,cdotc,clartg,crot - 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('CQRQHV',info) - end if -c form each element of w = Q'*u when necessary. - rr = cdotc(m,Q(1,k),1,u,1) - do i = k-1,1,-1 - w1 = rr - w = cdotc(m,Q(1,i),1,u,1) - call clartg(w,w1,c,s,rr) -c apply rotation to rows of R if necessary - if (i <= n) then - call crot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s) - end if -c apply rotation to columns of Q - call crot(m,Q(1,i),1,Q(1,i+1),1,c,conjg(s)) - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/cqrshc.f --- a/libcruft/qrupdate/cqrshc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -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 cqrshc(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 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 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 (complex 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 - complex Q(m,k),R(k,n) - external xerbla,cswap,cqhqr,cqrqhu - 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('CQRSHC',info) - end if - - if (i < j) then -c shift columns - do l = i,j-1 - call cswap(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 cqhqr(m,n+1-i,kk+1-i,Q(1,i),m,R(i,i),k) - else - call cqhqr(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 cswap(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 cqrqhu(m,n-j,kk+1-j,Q(1,j),m,R(j,jj),k,R(j,j),w) - else - call cqrqhu(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 diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dch1dn.f --- a/libcruft/qrupdate/dch1dn.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic -c -c Author: Jaroslav Hajek -c -c This source is free software; you can redistribute it and/or modify -c it under the terms of the GNU General Public License as published by -c the Free Software Foundation; either version 2 of the License, or -c (at your option) any later version. -c -c This program is distributed in the hope that it will be useful, -c but WITHOUT ANY WARRANTY; without even the implied warranty of -c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -c GNU General Public License for more details. -c -c You should have received a copy of the GNU General Public License -c along with this software; see the file COPYING. If not, see -c . -c - subroutine dch1dn(n,R,u,w,info) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a symmetric positive definite matrix A, i.e. -c A = R'*R, this subroutine downdates R -> R1 so that -c R1'*R1 = A - u*u' -c (real version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n,info - double precision R(n,n),u(n) - double precision w(n) - external dtrsv,dlartg,dnrm2 - double precision rho,dnrm2 - double precision rr,ui,t - integer i,j - -c quick return if possible - if (n <= 0) return -c check for singularity of R - do i = 1,n - if (R(i,i) == 0d0) then - info = 2 - return - end if - end do -c form R' \ u - call dtrsv('U','T','N',n,R,n,u,1) - rho = dnrm2(n,u,1) -c check positive definiteness - rho = 1 - rho**2 - if (rho <= 0d0) then - info = 1 - return - end if - rho = sqrt(rho) -c eliminate R' \ u - do i = n,1,-1 - ui = u(i) -c generate next rotation - call dlartg(rho,ui,w(i),u(i),rr) - rho = rr - end do -c apply rotations - do i = n,1,-1 - ui = 0d0 - do j = i,1,-1 - t = w(j)*ui + u(j)*R(j,i) - R(j,i) = w(j)*R(j,i) - u(j)*ui - ui = t - end do - end do - - info = 0 - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dch1up.f --- a/libcruft/qrupdate/dch1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,57 +0,0 @@ -c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic -c -c Author: Jaroslav Hajek -c -c This source is free software; you can redistribute it and/or modify -c it under the terms of the GNU General Public License as published by -c the Free Software Foundation; either version 2 of the License, or -c (at your option) any later version. -c -c This program is distributed in the hope that it will be useful, -c but WITHOUT ANY WARRANTY; without even the implied warranty of -c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -c GNU General Public License for more details. -c -c You should have received a copy of the GNU General Public License -c along with this software; see the file COPYING. If not, see -c . -c - - subroutine dch1up(n,R,u,w) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a symmetric positive definite matrix A, i.e. -c A = R'*R, this subroutine updates R -> R1 so that -c R1'*R1 = A + u*u' -c (real version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n - double precision R(n,n),u(n) - double precision w(n) - external dlartg - double precision rr,ui,t - integer i,j - - do i = 1,n -c apply stored rotations, column-wise - ui = u(i) - do j = 1,i-1 - t = w(j)*R(j,i) + u(j)*ui - ui = w(j)*ui - u(j)*R(j,i) - R(j,i) = t - end do -c generate next rotation - call dlartg(R(i,i),ui,w(i),u(i),rr) - R(i,i) = rr - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dchdex.f --- a/libcruft/qrupdate/dchdex.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,61 +0,0 @@ -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 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 xerbla,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('DCHDEX',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 diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dchinx.f --- a/libcruft/qrupdate/dchinx.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,108 +0,0 @@ -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 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 xerbla,dcopy,dlacpy,dtrsv,dnrm2,dqrqhu - 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('DCHINX',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 diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqhqr.f --- a/libcruft/qrupdate/dqhqr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,69 +0,0 @@ -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 xerbla,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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqr1up.f --- a/libcruft/qrupdate/dqr1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,52 +0,0 @@ -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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqrdec.f --- a/libcruft/qrupdate/dqrdec.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -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 xerbla,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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqrder.f --- a/libcruft/qrupdate/dqrder.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqrinc.f --- a/libcruft/qrupdate/dqrinc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -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 -c orthogonal, 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,dgemv - 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('DQRINC',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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqrinr.f --- a/libcruft/qrupdate/dqrinr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqrqhu.f --- a/libcruft/qrupdate/dqrqhu.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -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 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 xerbla,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 - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqrqhv.f --- a/libcruft/qrupdate/dqrqhv.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -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 xerbla,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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/dqrshc.f --- a/libcruft/qrupdate/dqrshc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -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 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 orthogonal, R1 upper -c trapezoidal, 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 xerbla,dswap,dqhqr,dqrqhu - 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 diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sch1dn.f --- a/libcruft/qrupdate/sch1dn.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -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 sch1dn(n,R,u,w,info) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a symmetric positive definite matrix A, i.e. -c A = R'*R, this subroutine downdates R -> R1 so that -c R1'*R1 = A - u*u' -c (real version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n,info - real R(n,n),u(n) - real w(n) - external strsv,slartg,snrm2 - real rho,snrm2 - real rr,ui,t - integer i,j - -c quick return if possible - if (n <= 0) return -c check for singularity of R - do i = 1,n - if (R(i,i) == 0e0) then - info = 2 - return - end if - end do -c form R' \ u - call strsv('U','T','N',n,R,n,u,1) - rho = snrm2(n,u,1) -c check positive definiteness - rho = 1 - rho**2 - if (rho <= 0e0) then - info = 1 - return - end if - rho = sqrt(rho) -c eliminate R' \ u - do i = n,1,-1 - ui = u(i) -c generate next rotation - call slartg(rho,ui,w(i),u(i),rr) - rho = rr - end do -c apply rotations - do i = n,1,-1 - ui = 0e0 - do j = i,1,-1 - t = w(j)*ui + u(j)*R(j,i) - R(j,i) = w(j)*R(j,i) - u(j)*ui - ui = t - end do - end do - - info = 0 - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sch1up.f --- a/libcruft/qrupdate/sch1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,57 +0,0 @@ -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 sch1up(n,R,u,w) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a symmetric positive definite matrix A, i.e. -c A = R'*R, this subroutine updates R -> R1 so that -c R1'*R1 = A + u*u' -c (real version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n - real R(n,n),u(n) - real w(n) - external slartg - real rr,ui,t - integer i,j - - do i = 1,n -c apply stored rotations, column-wise - ui = u(i) - do j = 1,i-1 - t = w(j)*R(j,i) + u(j)*ui - ui = w(j)*ui - u(j)*R(j,i) - R(j,i) = t - end do -c generate next rotation - call slartg(R(i,i),ui,w(i),u(i),rr) - R(i,i) = rr - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/schdex.f --- a/libcruft/qrupdate/schdex.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,61 +0,0 @@ -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 schdex(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 - real R(n,n),R1(n-1,n-1) - real Qdum,c,s,rr - external xerbla,slacpy,sqhqr,slartg - -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('SQRDEX',info) - end if - -c setup the new matrix R1 - if (j > 1) then - call slacpy('0',n-1,j-1,R(1,1),n,R1(1,1),n-1) - end if - if (j < n) then - call slacpy('0',n-1,n-j,R(1,j+1),n,R1(1,j),n-1) - call sqhqr(0,n-j,n-j,Qdum,1,R1(j,j),n-1) -c eliminate R(n,n) - call slartg(R1(n-1,n-1),R(n,n),c,s,rr) - R1(n-1,n-1) = rr - endif - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/schinx.f --- a/libcruft/qrupdate/schinx.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,108 +0,0 @@ -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 schinx(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 - real R(n,n),R1(n+1,n+1),u(n+1) - real rho,Qdum,w,snrm2 - external xerbla,scopy,slacpy,strsv,snrm2,sqrqhu - 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('SCHINX',info) - end if - -c copy shifted vector - if (j > 1) then - call scopy(j-1,u,1,R1(1,j),1) - end if - w = u(j) - if (j < n+1) then - call scopy(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) == 0e0) then - info = 2 - return - end if - end do -c form R' \ u - call strsv('U','T','N',n,R,n,R1(1,j),1) - rho = snrm2(n,R1(1,j),1) -c check positive definiteness - rho = u(j) - rho**2 - if (rho <= 0e0) 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) = 0e0 - end do - if (j > 1) then - call slacpy('0',n,j-1,R(1,1),n,R1(1,1),n+1) - end if - if (j <= n) then - call slacpy('0',n,n-j+1,R(1,j),n,R1(1,j+1),n+1) -c retriangularize - jj = min(j+1,n) - call sqrqhu(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) = 0e0 - end do - end if - - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqhqr.f --- a/libcruft/qrupdate/sqhqr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,69 +0,0 @@ -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 sqhqr(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 - real Q(ldq,*),R(ldr,*) - real c - real s,rr - external xerbla,slartg,srot - 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('SQHQR',info) - end if -c triangularize - do i = 1,min(k-1,n) - call slartg(R(i,i),R(i+1,i),c,s,rr) - R(i,i) = rr - R(i+1,i) = 0e0 - if (i < n) then - call srot(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 srot(m,Q(1,i),1,Q(1,i+1),1,c,s) - end if - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqr1up.f --- a/libcruft/qrupdate/sqr1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,52 +0,0 @@ -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 sqr1up(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 - real Q(m,k),R(k,n),u(m),v(n) - real w - external sqrqhv,sqhqr,saxpy -c quick return if possible - if (m <= 0 .or. n <= 0) return -c eliminate tail of Q'*u - call sqrqhv(m,n,k,Q,m,R,m,u,w) -c update R - - call saxpy(n,w,v,1,R(1,1),m) - -c retriangularize R - call sqhqr(m,n,k,Q,m,R,k) - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqrdec.f --- a/libcruft/qrupdate/sqrdec.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -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 sqrdec(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 - real Q(m,k),R(k,n),R1(k,n-1) - external xerbla,scopy,sqhqr - 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('SQRDEC',info) - end if -c copy leading portion - call scopy(k*(j-1),R,1,R1,1) - if (j < n) then -c copy trailing portion of R - call scopy(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 sqhqr(m,n-j,k-j+1,Q(1,j),m,R1(j,j),k) - end if - end if - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqrder.f --- a/libcruft/qrupdate/sqrder.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -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 sqrder(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 - real Q(m,m),Q1(m-1,m-1),R(m,n),R1(m-1,n) - real c - real s,rr,w - external xerbla,slacpy,slartg,srot,sscal,saxpy - 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('SQRDER',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 slacpy('0',j-1,m-1,Q(1,2),m,Q1(1,1),m-1) - end if - if (j < m) then - call slacpy('0',m-j,m-1,Q(j+1,2),m,Q1(j,1),m-1) - end if -c setup the new matrix R1 - call slacpy('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 slartg(Q(j,i),w,c,s,rr) - w = rr -c apply rotation to rows of R1 - if (i <= n) then - call srot(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 srot(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 slartg(Q(j,1),w,c,s,rr) - w = rr - call sscal(n,c,R1(1,1),m-1) - call saxpy(n,-s,R(1,1),m,R1(1,1),m-1) -c apply rotation to columns of Q1 - call sscal(m-1,c,Q1(1,1),1) - if (j > 1) then - call saxpy(j-1,-s,Q(1,1),1,Q1(1,1),1) - end if - if (j < m) then - call saxpy(m-j,-s,Q(j+1,1),1,Q1(j,1),1) - end if - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqrinc.f --- a/libcruft/qrupdate/sqrinc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -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 sqrinc(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 -c orthogonal, 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 - real Q(m,k),R(k,n),R1(k,n+1),x(m) - - - real w - external xerbla,scopy,sqrqhv,sgemv - 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('SQRINC',info) - end if -c copy leading portion of R - call scopy(k*(j-1),R,1,R1,1) - if (j <= n) then - call scopy(k*(n+1-j),R(1,j),1,R1(1,j+1),1) - end if - call sgemv('T',m,min(k,j-1),1e0,Q,m,x,1,0e0,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 sqrqhv(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) = 0e0 - end do - end if - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqrinr.f --- a/libcruft/qrupdate/sqrinr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -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 sqrinr(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 - real Q(m,m),Q1(m+1,m+1),R(m,n),R1(m+1,n),x(n) - external xerbla,slacpy,scopy,sqhqr - 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('SQRINR',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 slacpy('0',j-1,m,Q(1,1),m,Q1(1,2),m+1) - end if - if (j <= m) then - call slacpy('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) = 0e0 - Q1(j,i) = 0e0 - end do - Q1(j,1) = 1e0 -c setup the new matrix R1 - call scopy(n,x,1,R1(1,1),m+1) - call slacpy('0',m,n,R(1,1),m,R1(2,1),m+1) -c rotate to form proper QR - call sqhqr(m+1,n,m+1,Q1,m+1,R1,m+1) - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqrqhu.f --- a/libcruft/qrupdate/sqrqhu.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -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 sqrqhu(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 - real Q(ldq,*),R(ldr,*),u(*),rr - real c - real s,w - external xerbla,slartg,srot - 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('SQRQHU',info) - end if - rr = u(k) - do i = k-1,1,-1 - w = rr - if (w /= 0e0) then - call slartg(u(i),w,c,s,rr) -c apply rotation to rows of R if necessary - if (i <= n) then - call srot(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 srot(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 - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqrqhv.f --- a/libcruft/qrupdate/sqrqhv.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -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 sqrqhv(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 - real Q(ldq,*),R(ldr,*),u(*),rr - real c - real s,w,w1,sdot - external xerbla,sdot,slartg,srot - 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('SQRQHV',info) - end if -c form each element of w = Q'*u when necessary. - rr = sdot(m,Q(1,k),1,u,1) - do i = k-1,1,-1 - w1 = rr - w = sdot(m,Q(1,i),1,u,1) - call slartg(w,w1,c,s,rr) -c apply rotation to rows of R if necessary - if (i <= n) then - call srot(n+1-i,R(i,i),ldr,R(i+1,i),ldr,c,s) - end if -c apply rotation to columns of Q - call srot(m,Q(1,i),1,Q(1,i+1),1,c,s) - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/sqrshc.f --- a/libcruft/qrupdate/sqrshc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -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 sqrshc(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 orthogonal, R1 upper -c trapezoidal, 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 - real Q(m,k),R(k,n) - external xerbla,sswap,sqhqr,sqrqhu - real 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('SQRSHC',info) - end if - - if (i < j) then -c shift columns - do l = i,j-1 - call sswap(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 sqhqr(m,n+1-i,kk+1-i,Q(1,i),m,R(i,i),k) - else - call sqhqr(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 sswap(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 sqrqhu(m,n-j,kk+1-j,Q(1,j),m,R(j,jj),k,R(j,j),w) - else - call sqrqhu(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 diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zch1dn.f --- a/libcruft/qrupdate/zch1dn.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic -c -c Author: Jaroslav Hajek -c -c This source is free software; you can redistribute it and/or modify -c it under the terms of the GNU General Public License as published by -c the Free Software Foundation; either version 2 of the License, or -c (at your option) any later version. -c -c This program is distributed in the hope that it will be useful, -c but WITHOUT ANY WARRANTY; without even the implied warranty of -c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -c GNU General Public License for more details. -c -c You should have received a copy of the GNU General Public License -c along with this software; see the file COPYING. If not, see -c . -c - subroutine zch1dn(n,R,u,w,info) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a hermitian positive definite matrix A, i.e. -c A = R'*R, this subroutine downdates R -> R1 so that -c R1'*R1 = A - u*u' -c (complex version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n,info - double complex R(n,n),u(n) - double precision w(n) - external ztrsv,zlartg,dznrm2 - double precision rho,dznrm2 - double complex crho,rr,ui,t - integer i,j - -c quick return if possible - if (n <= 0) return -c check for singularity of R - do i = 1,n - if (R(i,i) == 0d0) then - info = 2 - return - end if - end do -c form R' \ u - call ztrsv('U','C','N',n,R,n,u,1) - rho = dznrm2(n,u,1) -c check positive definiteness - rho = 1 - rho**2 - if (rho <= 0d0) then - info = 1 - return - end if - crho = sqrt(rho) -c eliminate R' \ u - do i = n,1,-1 - ui = u(i) -c generate next rotation - call zlartg(crho,ui,w(i),u(i),rr) - crho = rr - end do -c apply rotations - do i = n,1,-1 - ui = 0d0 - do j = i,1,-1 - t = w(j)*ui + u(j)*R(j,i) - R(j,i) = w(j)*R(j,i) - conjg(u(j))*ui - ui = t - end do - end do - - info = 0 - end diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zch1up.f --- a/libcruft/qrupdate/zch1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,56 +0,0 @@ -c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic -c -c Author: Jaroslav Hajek -c -c This source is free software; you can redistribute it and/or modify -c it under the terms of the GNU General Public License as published by -c the Free Software Foundation; either version 2 of the License, or -c (at your option) any later version. -c -c This program is distributed in the hope that it will be useful, -c but WITHOUT ANY WARRANTY; without even the implied warranty of -c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -c GNU General Public License for more details. -c -c You should have received a copy of the GNU General Public License -c along with this software; see the file COPYING. If not, see -c . -c - subroutine zch1up(n,R,u,w) -c purpose: given an upper triangular matrix R that is a Cholesky -c factor of a hermitian positive definite matrix A, i.e. -c A = R'*R, this subroutine updates R -> R1 so that -c R1'*R1 = A + u*u' or A - u*u' -c (complex version) -c arguments: -c n (in) the order of matrix R -c R (io) on entry, the upper triangular matrix R -c on exit, the updated matrix R1 -c u (io) the vector determining the rank-1 update -c on exit, u is destroyed. -c w (w) a real workspace vector of size n -c -c NOTE: the workspace vector is used to store the rotations -c so that R does not need to be traversed by rows. -c - integer n - double complex R(n,n),u(n) - double precision w(n) - external zlartg - double complex rr,ui,t - integer i,j - - do i = 1,n -c apply stored rotations, column-wise - ui = conjg(u(i)) - do j = 1,i-1 - t = w(j)*R(j,i) + u(j)*ui - ui = w(j)*ui - conjg(u(j))*R(j,i) - R(j,i) = t - end do -c generate next rotation - call zlartg(R(i,i),ui,w(i),u(i),rr) - R(i,i) = rr - end do - end - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zchdex.f --- a/libcruft/qrupdate/zchdex.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,62 +0,0 @@ -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 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 (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 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 xerbla,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('ZCHDEX',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 diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zchinx.f --- a/libcruft/qrupdate/zchinx.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,109 +0,0 @@ -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 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 xerbla,zcopy,zlacpy,ztrsv,dznrm2,zqrqhu - integer jj - -c quick return if possible - if (n == 0) then - if (dble(u(1)) <= 0) then - info = 1 - return - else - R(1,1) = sqrt(dble(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('ZCHINX',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 diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqhqr.f --- a/libcruft/qrupdate/zqhqr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,69 +0,0 @@ -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 xerbla,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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqr1up.f --- a/libcruft/qrupdate/zqr1up.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqrdec.f --- a/libcruft/qrupdate/zqrdec.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -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 xerbla,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('ZQRDEC',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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqrder.f --- a/libcruft/qrupdate/zqrder.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -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,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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqrinc.f --- a/libcruft/qrupdate/zqrinc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,74 +0,0 @@ -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,zgemv - 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('ZQRINC',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),dcmplx(1d0),Q,m,x,1, - + dcmplx(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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqrinr.f --- a/libcruft/qrupdate/zqrinr.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -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 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 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,zcopy,zqhqr - 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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqrqhu.f --- a/libcruft/qrupdate/zqrqhu.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -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 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 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 u (in) the k-vector u. -c rr (out) the first element of Q1'*u on exit. -c -c if Q is unitary, 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 xerbla,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 - diff -r 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqrqhv.f --- a/libcruft/qrupdate/zqrqhv.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -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 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 u (in) the m-vector u. -c rr (out) the first element of Q1'*u on exit. -c -c if Q is unitary, 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 xerbla,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 3d8a914c580e -r d66c9b6e506a libcruft/qrupdate/zqrshc.f --- a/libcruft/qrupdate/zqrshc.f Tue Jan 20 21:15:17 2009 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -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 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 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 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 (complex 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 xerbla,zswap,zqhqr,zqrqhu - 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 diff -r 3d8a914c580e -r d66c9b6e506a liboctave/ChangeLog --- a/liboctave/ChangeLog Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/ChangeLog Tue Jan 20 21:16:42 2009 +0100 @@ -1,3 +1,71 @@ +2009-01-17 Jaroslav Hajek + + * floatQR.h (FloatQR::update, FloatQR::insert_col, + FloatQR::insert_row, FloatQR::delete_col, FloatQR::delete_row, + FloatQR::shift_col): Update interfaces. + + * floatQR.cc: Update external decls for qrupdate routines. + (FloatQR::update, FloatQR::insert_col, FloatQR::insert_row, + FloatQR::delete_col, FloatQR::delete_row, FloatQR::shift_col): Reflect + changes in qrupdate interfaces, implement batch updates. + + * dbleQR.h (QR::update, QR::insert_col, QR::insert_row, + QR::delete_col, QR::delete_row, QR::shift_col): Update interfaces. + + * dbleQR.cc: Update external decls for qrupdate routines. + (QR::update, QR::insert_col, QR::insert_row, QR::delete_col, + QR::delete_row, QR::shift_col): Reflect changes in qrupdate + interfaces, implement batch updates. + + * fCmplxQR.h (FloatComplexQR::update, FloatComplexQR::insert_col, + FloatComplexQR::insert_row, FloatComplexQR::delete_col, + FloatComplexQR::delete_row, FloatComplexQR::shift_col): Update + interfaces. + + * fCmplxQR.cc: Update external decls for qrupdate routines. + (FloatComplexQR::update, FloatComplexQR::insert_col, + FloatComplexQR::insert_row, FloatComplexQR::delete_col, + FloatComplexQR::delete_row, FloatComplexQR::shift_col): Reflect + changes in qrupdate interfaces, + implement batch updates. + + * CmplxQR.h (ComplexQR::update, ComplexQR::insert_col, + ComplexQR::insert_row, ComplexQR::delete_col, ComplexQR::delete_row, + ComplexQR::shift_col): Update interfaces. + + * CmplxQR.cc: Update external decls for qrupdate routines. + (ComplexQR::update, ComplexQR::insert_col, + ComplexQR::insert_row, ComplexQR::delete_col, ComplexQR::delete_row, + ComplexQR::shift_col): Reflect changes in qrupdate interfaces, + implement batch updates. + + * floatCHOL.h (FloatCHOL::update, FloatCHOL::downdate, + FloatCHOL::insert_sym): Update interfaces. + * floatCHOL.cc: Update external decls for qrupdate routines. + (FloatCHOL::update, FloatCHOL::downdate, FloatCHOL::insert_sym, + FloatCHOL::delete_sym, FloatCHOL::shift_sym): Reflect changes in + qrupdate interfaces, + + * CHOL.h (CHOL::update, CHOL::downdate, CHOL::insert_sym): Update + interfaces. + * CHOL.cc: Update external decls for qrupdate routines. + (CHOL::update, CHOL::downdate, CHOL::insert_sym, CHOL::delete_sym, + CHOL::shift_sym): Reflect changes in qrupdate interfaces, + + * fCmplxCHOL.h (FloatComplexCHOL::update, FloatComplexCHOL::downdate, + FloatComplexCHOL::insert_sym): Update interfaces. + * fCmplxCHOL.cc: Update external decls for qrupdate routines. + (FloatComplexCHOL::update, FloatComplexCHOL::downdate, + FloatComplexCHOL::insert_sym, FloatComplexCHOL::delete_sym, + FloatComplexCHOL::shift_sym): Reflect changes in qrupdate interfaces, + + * CmplxCHOL.h (ComplexCHOL::update, ComplexCHOL::downdate, + ComplexCHOL::insert_sym): Update interfaces. + * CmplxCHOL.cc: Update external decls for qrupdate routines. + (ComplexCHOL::update, ComplexCHOL::downdate, ComplexCHOL::insert_sym, + ComplexCHOL::delete_sym, ComplexCHOL::shift_sym): Reflect changes in + qrupdate interfaces, + 2009-01-17 Jaroslav Hajek * Array.h (Array): Document internal use of slice_data and diff -r 3d8a914c580e -r d66c9b6e506a liboctave/CmplxCHOL.cc --- a/liboctave/CmplxCHOL.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/CmplxCHOL.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #ifdef HAVE_CONFIG_H #include #endif @@ -52,23 +51,30 @@ Complex*, const octave_idx_type&, const double&, double&, Complex*, double*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); +#ifdef HAVE_QRUPDATE + F77_RET_T - F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, Complex*, double*); + F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, const octave_idx_type&, + Complex*, double*); F77_RET_T - F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, Complex*, double*, + F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, const octave_idx_type&, + Complex*, double*, octave_idx_type&); + + F77_RET_T + F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, 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_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, double*); 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&); + F77_FUNC (zchshx, ZCHSHX) (const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + Complex*, double*); +#endif } octave_idx_type @@ -182,26 +188,28 @@ (*current_liboctave_error_handler) ("CHOL requires square matrix"); } +#ifdef HAVE_QRUPDATE + void -ComplexCHOL::update (const ComplexMatrix& u) +ComplexCHOL::update (const ComplexColumnVector& u) { octave_idx_type n = chol_mat.rows (); if (u.length () == n) { - ComplexMatrix tmp = u; + ComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w)); + F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), rw)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); } octave_idx_type -ComplexCHOL::downdate (const ComplexMatrix& u) +ComplexCHOL::downdate (const ComplexColumnVector& u) { octave_idx_type info = -1; @@ -209,38 +217,40 @@ if (u.length () == n) { - ComplexMatrix tmp = u; + ComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (double, w, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w, info)); + F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), rw, info)); } else - (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); return info; } octave_idx_type -ComplexCHOL::insert_sym (const ComplexMatrix& u, octave_idx_type j) +ComplexCHOL::insert_sym (const ComplexColumnVector& 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"); + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("CHOL insert index out of range"); + (*current_liboctave_error_handler) ("cholinsert: index out of range"); else { - ComplexMatrix chol_mat1 (n+1, n+1); + ComplexColumnVector utmp = u; + + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zchinx, ZCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), - j+1, u.data (), info)); + chol_mat.resize (n+1, n+1); - chol_mat = chol_mat1; + F77_XFCN (zchinx, ZCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), rw, info)); } return info; @@ -252,14 +262,15 @@ octave_idx_type n = chol_mat.rows (); if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("CHOL delete index out of range"); + (*current_liboctave_error_handler) ("choldelete: index out of range"); else { - ComplexMatrix chol_mat1 (n-1, n-1); + OCTAVE_LOCAL_BUFFER (double, rw, n); - F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, rw)); - chol_mat = chol_mat1; + chol_mat.resize (n-1, n-1); } } @@ -267,14 +278,21 @@ 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"); + (*current_liboctave_error_handler) ("cholshift: index out of range"); else - F77_XFCN (zqrshc, ZQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (Complex, w, n); + OCTAVE_LOCAL_BUFFER (double, rw, n); + + F77_XFCN (zchshx, ZCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w, rw)); + } } +#endif + ComplexMatrix chol2inv (const ComplexMatrix& r) { diff -r 3d8a914c580e -r d66c9b6e506a liboctave/CmplxCHOL.h --- a/liboctave/CmplxCHOL.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/CmplxCHOL.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,14 +22,13 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #if !defined (octave_ComplexCHOL_h) #define octave_ComplexCHOL_h 1 #include #include "CMatrix.h" +#include "CColVector.h" class OCTAVE_API @@ -67,16 +67,20 @@ void set (const ComplexMatrix& R); - void update (const ComplexMatrix& u); +#ifdef HAVE_QRUPDATE + + void update (const ComplexColumnVector& u); - octave_idx_type downdate (const ComplexMatrix& u); + octave_idx_type downdate (const ComplexColumnVector& u); - octave_idx_type insert_sym (const ComplexMatrix& u, octave_idx_type j); + octave_idx_type insert_sym (const ComplexColumnVector& u, octave_idx_type j); void delete_sym (octave_idx_type j); void shift_sym (octave_idx_type i, octave_idx_type j); +#endif + friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexCHOL& a); private: diff -r 3d8a914c580e -r d66c9b6e506a liboctave/CmplxQR.cc --- a/liboctave/CmplxQR.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/CmplxQR.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #ifdef HAVE_CONFIG_H #include #endif @@ -32,6 +31,7 @@ #include "lo-error.h" #include "Range.h" #include "idx-vector.h" +#include "oct-locbuf.h" extern "C" { @@ -45,33 +45,40 @@ Complex*, const octave_idx_type&, Complex*, Complex*, const octave_idx_type&, octave_idx_type&); - // these come from qrupdate +#ifdef HAVE_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*); + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + Complex*, Complex*, Complex*, double*); 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*); + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, const Complex*, double*); 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&); + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, double*); 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*); + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, const Complex*, double*); 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&); + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, Complex*, double*); 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&); + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + Complex*, double*); + +#endif } ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type) @@ -171,6 +178,27 @@ this->r = r_arg; } +#ifdef HAVE_QRUPDATE + +void +ComplexQR::update (const ComplexColumnVector& u, const ComplexColumnVector& v) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () == m && v.length () == n) + { + ComplexColumnVector utmp = u, vtmp = v; + OCTAVE_LOCAL_BUFFER (Complex, w, k); + OCTAVE_LOCAL_BUFFER (double, rw, k); + F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w, rw)); + } + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); +} + void ComplexQR::update (const ComplexMatrix& u, const ComplexMatrix& v) { @@ -178,32 +206,94 @@ octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); - if (u.length () == m && v.length () == n) - F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), - u.data (), v.data ())); + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + OCTAVE_LOCAL_BUFFER (Complex, w, k); + OCTAVE_LOCAL_BUFFER (double, rw, k); + for (octave_idx_type i = 0; i < u.cols (); i++) + { + ComplexColumnVector utmp = u.column (i), vtmp = v.column (i); + F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w, rw)); + } + } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void -ComplexQR::insert_col (const ComplexMatrix& u, octave_idx_type j) +ComplexQR::insert_col (const ComplexColumnVector& u, octave_idx_type j) { octave_idx_type m = q.rows (); octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); if (u.length () != m) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - ComplexMatrix r1 (m,n+1); + if (k < m) + { + q.resize (m, k+1); + r.resize (k+1, n+1); + } + else + { + r.resize (k, n+1); + } + + ComplexColumnVector utmp = u; + OCTAVE_LOCAL_BUFFER (double, rw, k); + F77_XFCN (zqrinc, ZQRINC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + utmp.data (), rw)); + } +} + +void +ComplexQR::insert_col (const ComplexMatrix& u, const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); - F77_XFCN (zqrinc, ZQRINC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1, u.data ())); + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); - r = r1; + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + octave_idx_type kmax = std::min (k + nj, m); + if (k < m) + { + q.resize (m, kmax); + r.resize (kmax, n + nj); + } + else + { + r.resize (k, n + nj); + } + + OCTAVE_LOCAL_BUFFER (double, rw, kmax); + for (octave_idx_type i = 0; i < js.length (); i++) + { + ComplexColumnVector utmp = u.column (jsi(i)); + F77_XFCN (zqrinc, ZQRINC, (m, n + i, std::min (kmax, k + i), + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, + utmp.data (), rw)); + } } } @@ -214,41 +304,87 @@ octave_idx_type k = r.rows (); octave_idx_type n = r.columns (); - if (k < m && k < n) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); - else if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - ComplexMatrix r1 (k, n-1); + OCTAVE_LOCAL_BUFFER (double, rw, k); + F77_XFCN (zqrdec, ZQRDEC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, rw)); - F77_XFCN (zqrdec, ZQRDEC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1)); - - r = r1; + if (k < m) + { + q.resize (m, k-1); + r.resize (k-1, n-1); + } + else + { + r.resize (k, n-1); + } } } void -ComplexQR::insert_row (const ComplexMatrix& u, octave_idx_type j) +ComplexQR::delete_col (const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + OCTAVE_LOCAL_BUFFER (double, rw, k); + for (octave_idx_type i = 0; i < js.length (); i++) + { + F77_XFCN (zqrdec, ZQRDEC, (m, n - i, k == m ? k : k - i, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, rw)); + } + if (k < m) + { + q.resize (m, k - nj); + r.resize (k - nj, n - nj); + } + else + { + r.resize (k, n - nj); + } + + } +} + +void +ComplexQR::insert_row (const ComplexRowVector& u, octave_idx_type j) { octave_idx_type m = r.rows (); octave_idx_type n = r.columns (); + octave_idx_type k = std::min (m, n); if (! q.is_square () || u.length () != n) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > m) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - ComplexMatrix q1 (m+1, m+1); - ComplexMatrix r1 (m+1, n); + q.resize (m + 1, m + 1); + r.resize (m + 1, n); + ComplexRowVector utmp = u; + OCTAVE_LOCAL_BUFFER (double, rw, k); + F77_XFCN (zqrinr, ZQRINR, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + j + 1, utmp.fortran_vec (), rw)); - F77_XFCN (zqrinr, ZQRINR, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1, u.data ())); - - q = q1; - r = r1; } } @@ -259,19 +395,19 @@ octave_idx_type n = r.columns (); if (! q.is_square ()) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); else if (j < 0 || j > m-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - ComplexMatrix q1 (m-1, m-1); - ComplexMatrix r1 (m-1, n); + OCTAVE_LOCAL_BUFFER (Complex, w, m); + OCTAVE_LOCAL_BUFFER (double, rw, m); + F77_XFCN (zqrder, ZQRDER, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + w, rw)); - F77_XFCN (zqrder, ZQRDER, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1 )); - - q = q1; - r = r1; + q.resize (m - 1, m - 1); + r.resize (m - 1, n); } } @@ -283,22 +419,19 @@ 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"); + (*current_liboctave_error_handler) ("qrshift: index out of range"); else - F77_XFCN (zqrshc, ZQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (Complex, w, k); + OCTAVE_LOCAL_BUFFER (double, rw, k); + F77_XFCN (zqrshc, ZQRSHC, (m, n, k, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + i + 1, j + 1, w, rw)); + } } -void -ComplexQR::economize (void) -{ - octave_idx_type r_nc = r.columns (); - - if (r.rows () > r_nc) - { - q.resize (q.rows (), r_nc); - r.resize (r_nc, r_nc); - } -} +#endif /* ;;; Local Variables: *** diff -r 3d8a914c580e -r d66c9b6e506a liboctave/CmplxQR.h --- a/liboctave/CmplxQR.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/CmplxQR.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #if !defined (octave_ComplexQR_h) #define octave_ComplexQR_h 1 @@ -66,19 +65,27 @@ ComplexMatrix R (void) const { return r; } +#ifdef HAVE_QRUPDATE + + void update (const ComplexColumnVector& u, const ComplexColumnVector& v); + void update (const ComplexMatrix& u, const ComplexMatrix& v); - void insert_col (const ComplexMatrix& u, octave_idx_type j); + void insert_col (const ComplexColumnVector& u, octave_idx_type j); + + void insert_col (const ComplexMatrix& u, const Array& j); void delete_col (octave_idx_type j); - void insert_row (const ComplexMatrix& u, octave_idx_type j); + void delete_col (const Array& j); + + void insert_row (const ComplexRowVector& u, octave_idx_type j); void delete_row (octave_idx_type j); void shift_cols (octave_idx_type i, octave_idx_type j); - void economize(); +#endif friend std::ostream& operator << (std::ostream&, const ComplexQR&); diff -r 3d8a914c580e -r d66c9b6e506a liboctave/dbleCHOL.cc --- a/liboctave/dbleCHOL.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/dbleCHOL.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #ifdef HAVE_CONFIG_H #include #endif @@ -52,23 +51,30 @@ double*, const octave_idx_type&, const double&, double&, double*, octave_idx_type*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); +#ifdef HAVE_QRUPDATE + F77_RET_T - F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, double*, double*); + F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, const octave_idx_type&, + double*, double*); F77_RET_T - F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, + F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, const octave_idx_type&, + double*, double*, octave_idx_type&); + + F77_RET_T + F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, 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_FUNC (dchdex, DCHDEX) (const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, double*); 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&); + F77_FUNC (dchshx, DCHSHX) (const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + double*); +#endif } octave_idx_type @@ -186,26 +192,28 @@ (*current_liboctave_error_handler) ("CHOL requires square matrix"); } +#ifdef HAVE_QRUPDATE + void -CHOL::update (const Matrix& u) +CHOL::update (const ColumnVector& u) { octave_idx_type n = chol_mat.rows (); if (u.length () == n) { - Matrix tmp = u; + ColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w)); + F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); } octave_idx_type -CHOL::downdate (const Matrix& u) +CHOL::downdate (const ColumnVector& u) { octave_idx_type info = -1; @@ -213,38 +221,40 @@ if (u.length () == n) { - Matrix tmp = u; + ColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w, info)); + F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); return info; } octave_idx_type -CHOL::insert_sym (const Matrix& u, octave_idx_type j) +CHOL::insert_sym (const ColumnVector& 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"); + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("CHOL insert index out of range"); + (*current_liboctave_error_handler) ("cholinsert: index out of range"); else { - Matrix chol_mat1 (n+1, n+1); + ColumnVector utmp = u; + + OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dchinx, DCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), - j+1, u.data (), info)); + chol_mat.resize (n+1, n+1); - chol_mat = chol_mat1; + F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), w, info)); } return info; @@ -256,14 +266,15 @@ octave_idx_type n = chol_mat.rows (); if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("CHOL delete index out of range"); + (*current_liboctave_error_handler) ("choldelete: index out of range"); else { - Matrix chol_mat1 (n-1, n-1); + OCTAVE_LOCAL_BUFFER (double, w, n); - F77_XFCN (dchdex, DCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, w)); - chol_mat = chol_mat1; + chol_mat.resize (n-1, n-1); } } @@ -271,14 +282,20 @@ 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"); + (*current_liboctave_error_handler) ("cholshift: index out of range"); else - F77_XFCN (dqrshc, DQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (double, w, 2*n); + + F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w)); + } } +#endif + Matrix chol2inv (const Matrix& r) { diff -r 3d8a914c580e -r d66c9b6e506a liboctave/dbleCHOL.h --- a/liboctave/dbleCHOL.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/dbleCHOL.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,14 +22,13 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #if !defined (octave_CHOL_h) #define octave_CHOL_h 1 #include #include "dMatrix.h" +#include "dColVector.h" class OCTAVE_API @@ -64,16 +64,20 @@ void set (const Matrix& R); - void update (const Matrix& u); +#ifdef HAVE_QRUPDATE + + void update (const ColumnVector& u); - octave_idx_type downdate (const Matrix& u); + octave_idx_type downdate (const ColumnVector& u); - octave_idx_type insert_sym (const Matrix& u, octave_idx_type j); + octave_idx_type insert_sym (const ColumnVector& u, octave_idx_type j); void delete_sym (octave_idx_type j); void shift_sym (octave_idx_type i, octave_idx_type j); +#endif + friend OCTAVE_API std::ostream& operator << (std::ostream& os, const CHOL& a); private: diff -r 3d8a914c580e -r d66c9b6e506a liboctave/dbleQR.cc --- a/liboctave/dbleQR.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/dbleQR.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,7 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton -Copyright (C) 2008 Jaroslav Hajek +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -31,6 +31,7 @@ #include "lo-error.h" #include "Range.h" #include "idx-vector.h" +#include "oct-locbuf.h" extern "C" { @@ -42,33 +43,40 @@ 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 +#ifdef HAVE_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*); + double*, const octave_idx_type&, double*, const octave_idx_type&, + double*, double*, 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*); + double*, const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, const double*, 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&); + double*, const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, double*); 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*); + double*, const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, const double*, 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&); + double*, const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, double*); 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&); + double*, const octave_idx_type&, double*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + double*); + +#endif } QR::QR (const Matrix& a, QR::type qr_type) @@ -161,6 +169,26 @@ this->r = r_arg; } +#ifdef HAVE_QRUPDATE + +void +QR::update (const ColumnVector& u, const ColumnVector& v) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () == m && v.length () == n) + { + ColumnVector utmp = u, vtmp = v; + OCTAVE_LOCAL_BUFFER (double, w, 2*k); + F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w)); + } + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); +} + void QR::update (const Matrix& u, const Matrix& v) { @@ -168,32 +196,93 @@ octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); - if (u.length () == m && v.length () == n) - F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), - u.data (), v.data ())); + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + OCTAVE_LOCAL_BUFFER (double, w, 2*k); + for (octave_idx_type i = 0; i < u.cols (); i++) + { + ColumnVector utmp = u.column (i), vtmp = v.column (i); + F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w)); + } + } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void -QR::insert_col (const Matrix& u, octave_idx_type j) +QR::insert_col (const ColumnVector& u, octave_idx_type j) { octave_idx_type m = q.rows (); octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); if (u.length () != m) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - Matrix r1 (m, n+1); + if (k < m) + { + q.resize (m, k+1); + r.resize (k+1, n+1); + } + else + { + r.resize (k, n+1); + } + + ColumnVector utmp = u; + OCTAVE_LOCAL_BUFFER (double, w, k); + F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + utmp.data (), w)); + } +} + +void +QR::insert_col (const Matrix& u, const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); - F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1, u.data ())); + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); - r = r1; + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + octave_idx_type kmax = std::min (k + nj, m); + if (k < m) + { + q.resize (m, kmax); + r.resize (kmax, n + nj); + } + else + { + r.resize (k, n + nj); + } + + OCTAVE_LOCAL_BUFFER (double, w, kmax); + for (octave_idx_type i = 0; i < js.length (); i++) + { + ColumnVector utmp = u.column (jsi(i)); + F77_XFCN (dqrinc, DQRINC, (m, n + i, std::min (kmax, k + i), + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, + utmp.data (), w)); + } } } @@ -204,41 +293,87 @@ octave_idx_type k = r.rows (); octave_idx_type n = r.columns (); - if (k < m && k < n) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); - else if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - Matrix r1 (k, n-1); + OCTAVE_LOCAL_BUFFER (double, w, k); + F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, w)); - F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1)); - - r = r1; + if (k < m) + { + q.resize (m, k-1); + r.resize (k-1, n-1); + } + else + { + r.resize (k, n-1); + } } } void -QR::insert_row (const Matrix& u, octave_idx_type j) +QR::delete_col (const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + OCTAVE_LOCAL_BUFFER (double, w, k); + for (octave_idx_type i = 0; i < js.length (); i++) + { + F77_XFCN (dqrdec, DQRDEC, (m, n - i, k == m ? k : k - i, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, w)); + } + if (k < m) + { + q.resize (m, k - nj); + r.resize (k - nj, n - nj); + } + else + { + r.resize (k, n - nj); + } + + } +} + +void +QR::insert_row (const RowVector& u, octave_idx_type j) { octave_idx_type m = r.rows (); octave_idx_type n = r.columns (); + octave_idx_type k = std::min (m, n); if (! q.is_square () || u.length () != n) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > m) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - Matrix q1 (m+1, m+1); - Matrix r1 (m+1, n); + q.resize (m + 1, m + 1); + r.resize (m + 1, n); + RowVector utmp = u; + OCTAVE_LOCAL_BUFFER (double, w, k); + F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + j + 1, utmp.fortran_vec (), w)); - F77_XFCN (dqrinr, DQRINR, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1, u.data ())); - - q = q1; - r = r1; } } @@ -249,19 +384,18 @@ octave_idx_type n = r.columns (); if (! q.is_square ()) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); else if (j < 0 || j > m-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - Matrix q1 (m-1, m-1); - Matrix r1 (m-1, n); + OCTAVE_LOCAL_BUFFER (double, w, 2*m); + F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + w)); - F77_XFCN (dqrder, DQRDER, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1 )); - - q = q1; - r = r1; + q.resize (m - 1, m - 1); + r.resize (m - 1, n); } } @@ -273,23 +407,18 @@ 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"); + (*current_liboctave_error_handler) ("qrshift: 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 (); - - if (r.rows () > r_nc) { - q.resize (q.rows (), r_nc); - r.resize (r_nc, r_nc); + OCTAVE_LOCAL_BUFFER (double, w, 2*k); + F77_XFCN (dqrshc, DQRSHC, (m, n, k, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + i + 1, j + 1, w)); } } +#endif /* ;;; Local Variables: *** diff -r 3d8a914c580e -r d66c9b6e506a liboctave/dbleQR.h --- a/liboctave/dbleQR.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/dbleQR.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #if !defined (octave_QR_h) #define octave_QR_h 1 @@ -71,19 +70,27 @@ Matrix R (void) const { return r; } +#ifdef HAVE_QRUPDATE + + void update (const ColumnVector& u, const ColumnVector& v); + void update (const Matrix& u, const Matrix& v); - void insert_col (const Matrix& u, octave_idx_type j); + void insert_col (const ColumnVector& u, octave_idx_type j); + + void insert_col (const Matrix& u, const Array& j); void delete_col (octave_idx_type j); - void insert_row (const Matrix& u, octave_idx_type j); + void delete_col (const Array& j); + + void insert_row (const RowVector& u, octave_idx_type j); void delete_row (octave_idx_type j); void shift_cols (octave_idx_type i, octave_idx_type j); - void economize (void); +#endif friend std::ostream& operator << (std::ostream&, const QR&); diff -r 3d8a914c580e -r d66c9b6e506a liboctave/fCmplxCHOL.cc --- a/liboctave/fCmplxCHOL.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/fCmplxCHOL.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #ifdef HAVE_CONFIG_H #include #endif @@ -52,23 +51,30 @@ FloatComplex*, const octave_idx_type&, const float&, float&, FloatComplex*, float*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); +#ifdef HAVE_QRUPDATE + F77_RET_T - F77_FUNC (cch1up, CCH1UP) (const octave_idx_type&, FloatComplex*, FloatComplex*, float*); + F77_FUNC (cch1up, CCH1UP) (const octave_idx_type&, FloatComplex*, const octave_idx_type&, + FloatComplex*, float*); F77_RET_T - F77_FUNC (cch1dn, CCH1DN) (const octave_idx_type&, FloatComplex*, FloatComplex*, float*, + F77_FUNC (cch1dn, CCH1DN) (const octave_idx_type&, FloatComplex*, const octave_idx_type&, + FloatComplex*, float*, octave_idx_type&); + + F77_RET_T + F77_FUNC (cchinx, CCHINX) (const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, FloatComplex*, float*, octave_idx_type&); F77_RET_T - F77_FUNC (cqrshc, CQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, FloatComplex*, const octave_idx_type&, const octave_idx_type&); + F77_FUNC (cchdex, CCHDEX) (const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, float*); F77_RET_T - F77_FUNC (cchinx, CCHINX) (const octave_idx_type&, const FloatComplex*, FloatComplex*, const octave_idx_type&, - const FloatComplex*, octave_idx_type&); - - F77_RET_T - F77_FUNC (cchdex, CCHDEX) (const octave_idx_type&, const FloatComplex*, FloatComplex*, const octave_idx_type&); + F77_FUNC (cchshx, CCHSHX) (const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + FloatComplex*, float*); +#endif } octave_idx_type @@ -182,26 +188,28 @@ (*current_liboctave_error_handler) ("CHOL requires square matrix"); } +#ifdef HAVE_QRUPDATE + void -FloatComplexCHOL::update (const FloatComplexMatrix& u) +FloatComplexCHOL::update (const FloatComplexColumnVector& u) { octave_idx_type n = chol_mat.rows (); if (u.length () == n) { - FloatComplexMatrix tmp = u; + FloatComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, w, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w)); + F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), rw)); } else - (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); } octave_idx_type -FloatComplexCHOL::downdate (const FloatComplexMatrix& u) +FloatComplexCHOL::downdate (const FloatComplexColumnVector& u) { octave_idx_type info = -1; @@ -209,38 +217,40 @@ if (u.length () == n) { - FloatComplexMatrix tmp = u; + FloatComplexColumnVector utmp = u; - OCTAVE_LOCAL_BUFFER (float, w, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w, info)); + F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), rw, info)); } else - (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); return info; } octave_idx_type -FloatComplexCHOL::insert_sym (const FloatComplexMatrix& u, octave_idx_type j) +FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& 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"); + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("CHOL insert index out of range"); + (*current_liboctave_error_handler) ("cholinsert: index out of range"); else { - FloatComplexMatrix chol_mat1 (n+1, n+1); + FloatComplexColumnVector utmp = u; + + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cchinx, CCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), - j+1, u.data (), info)); + chol_mat.resize (n+1, n+1); - chol_mat = chol_mat1; + F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), rw, info)); } return info; @@ -252,14 +262,15 @@ octave_idx_type n = chol_mat.rows (); if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("CHOL delete index out of range"); + (*current_liboctave_error_handler) ("choldelete: index out of range"); else { - FloatComplexMatrix chol_mat1 (n-1, n-1); + OCTAVE_LOCAL_BUFFER (float, rw, n); - F77_XFCN (cchdex, CCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, rw)); - chol_mat = chol_mat1; + chol_mat.resize (n-1, n-1); } } @@ -267,14 +278,21 @@ FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j) { octave_idx_type n = chol_mat.rows (); - FloatComplex dummy; if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("CHOL shift index out of range"); + (*current_liboctave_error_handler) ("cholshift: index out of range"); else - F77_XFCN (cqrshc, CQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (FloatComplex, w, n); + OCTAVE_LOCAL_BUFFER (float, rw, n); + + F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w, rw)); + } } +#endif + FloatComplexMatrix chol2inv (const FloatComplexMatrix& r) { diff -r 3d8a914c580e -r d66c9b6e506a liboctave/fCmplxCHOL.h --- a/liboctave/fCmplxCHOL.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/fCmplxCHOL.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,14 +22,13 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #if !defined (octave_FloatComplexCHOL_h) #define octave_FloatComplexCHOL_h 1 #include #include "fCMatrix.h" +#include "fCColVector.h" class OCTAVE_API @@ -67,16 +67,20 @@ void set (const FloatComplexMatrix& R); - void update (const FloatComplexMatrix& u); +#ifdef HAVE_QRUPDATE + + void update (const FloatComplexColumnVector& u); - octave_idx_type downdate (const FloatComplexMatrix& u); + octave_idx_type downdate (const FloatComplexColumnVector& u); - octave_idx_type insert_sym (const FloatComplexMatrix& u, octave_idx_type j); + octave_idx_type insert_sym (const FloatComplexColumnVector& u, octave_idx_type j); void delete_sym (octave_idx_type j); void shift_sym (octave_idx_type i, octave_idx_type j); +#endif + friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatComplexCHOL& a); private: diff -r 3d8a914c580e -r d66c9b6e506a liboctave/fCmplxQR.cc --- a/liboctave/fCmplxQR.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/fCmplxQR.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #ifdef HAVE_CONFIG_H #include #endif @@ -32,6 +31,7 @@ #include "lo-error.h" #include "Range.h" #include "idx-vector.h" +#include "oct-locbuf.h" extern "C" { @@ -45,33 +45,40 @@ FloatComplex*, const octave_idx_type&, FloatComplex*, FloatComplex*, const octave_idx_type&, octave_idx_type&); - // these come from qrupdate +#ifdef HAVE_QRUPDATE F77_RET_T F77_FUNC (cqr1up, CQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, FloatComplex*, const FloatComplex*, const FloatComplex*); + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + FloatComplex*, FloatComplex*, FloatComplex*, float*); F77_RET_T F77_FUNC (cqrinc, CQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const FloatComplex*, FloatComplex*, const octave_idx_type&, const FloatComplex*); + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, const FloatComplex*, float*); F77_RET_T F77_FUNC (cqrdec, CQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const FloatComplex*, FloatComplex*, const octave_idx_type&); + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, float*); F77_RET_T F77_FUNC (cqrinr, CQRINR) (const octave_idx_type&, const octave_idx_type&, - const FloatComplex*, FloatComplex*, const FloatComplex*, FloatComplex*, - const octave_idx_type&, const FloatComplex*); + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, const FloatComplex*, float*); F77_RET_T F77_FUNC (cqrder, CQRDER) (const octave_idx_type&, const octave_idx_type&, - const FloatComplex*, FloatComplex*, const FloatComplex*, FloatComplex *, - const octave_idx_type&); + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, FloatComplex*, float*); F77_RET_T F77_FUNC (cqrshc, CQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, FloatComplex*, const octave_idx_type&, const octave_idx_type&); + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + FloatComplex*, float*); + +#endif } FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& a, QR::type qr_type) @@ -171,6 +178,27 @@ this->r = r_arg; } +#ifdef HAVE_QRUPDATE + +void +FloatComplexQR::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () == m && v.length () == n) + { + FloatComplexColumnVector utmp = u, vtmp = v; + OCTAVE_LOCAL_BUFFER (FloatComplex, w, k); + OCTAVE_LOCAL_BUFFER (float, rw, k); + F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w, rw)); + } + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); +} + void FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v) { @@ -178,32 +206,93 @@ octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); - if (u.length () == m && v.length () == n) - F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), - u.data (), v.data ())); + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + OCTAVE_LOCAL_BUFFER (FloatComplex, w, k); + OCTAVE_LOCAL_BUFFER (float, rw, k); + for (octave_idx_type i = 0; i < u.cols (); i++) + { + FloatComplexColumnVector utmp = u.column (i), vtmp = v.column (i); + F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w, rw)); + } + } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void -FloatComplexQR::insert_col (const FloatComplexMatrix& u, octave_idx_type j) +FloatComplexQR::insert_col (const FloatComplexColumnVector& u, octave_idx_type j) { octave_idx_type m = q.rows (); octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); if (u.length () != m) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - FloatComplexMatrix r1 (m,n+1); + if (k < m) + { + q.resize (m, k+1); + r.resize (k+1, n+1); + } + else + { + r.resize (k, n+1); + } + + FloatComplexColumnVector utmp = u; + OCTAVE_LOCAL_BUFFER (float, rw, k); + F77_XFCN (cqrinc, CQRINC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + utmp.data (), rw)); + } +} + +void +FloatComplexQR::insert_col (const FloatComplexMatrix& u, const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); - F77_XFCN (cqrinc, CQRINC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1, u.data ())); + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); - r = r1; + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + octave_idx_type kmax = std::min (k + nj, m); + if (k < m) + { + q.resize (m, kmax); + r.resize (kmax, n + nj); + } + else + { + r.resize (k, n + nj); + } + + OCTAVE_LOCAL_BUFFER (float, rw, kmax); + for (octave_idx_type i = 0; i < js.length (); i++) + { + F77_XFCN (cqrinc, CQRINC, (m, n + i, std::min (kmax, k + i), + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, + u.column (jsi(i)).data (), rw)); + } } } @@ -214,41 +303,87 @@ octave_idx_type k = r.rows (); octave_idx_type n = r.columns (); - if (k < m && k < n) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); - else if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - FloatComplexMatrix r1 (k, n-1); + OCTAVE_LOCAL_BUFFER (float, rw, k); + F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, rw)); - F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1)); - - r = r1; + if (k < m) + { + q.resize (m, k-1); + r.resize (k-1, n-1); + } + else + { + r.resize (k, n-1); + } } } void -FloatComplexQR::insert_row (const FloatComplexMatrix& u, octave_idx_type j) +FloatComplexQR::delete_col (const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + OCTAVE_LOCAL_BUFFER (float, rw, k); + for (octave_idx_type i = 0; i < js.length (); i++) + { + F77_XFCN (cqrdec, CQRDEC, (m, n - i, k == m ? k : k - i, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, rw)); + } + if (k < m) + { + q.resize (m, k - nj); + r.resize (k - nj, n - nj); + } + else + { + r.resize (k, n - nj); + } + + } +} + +void +FloatComplexQR::insert_row (const FloatComplexRowVector& u, octave_idx_type j) { octave_idx_type m = r.rows (); octave_idx_type n = r.columns (); + octave_idx_type k = std::min (m, n); if (! q.is_square () || u.length () != n) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > m) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - FloatComplexMatrix q1 (m+1, m+1); - FloatComplexMatrix r1 (m+1, n); + q.resize (m + 1, m + 1); + r.resize (m + 1, n); + FloatComplexRowVector utmp = u; + OCTAVE_LOCAL_BUFFER (float, rw, k); + F77_XFCN (cqrinr, CQRINR, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + j + 1, utmp.fortran_vec (), rw)); - F77_XFCN (cqrinr, CQRINR, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1, u.data ())); - - q = q1; - r = r1; } } @@ -259,19 +394,19 @@ octave_idx_type n = r.columns (); if (! q.is_square ()) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); else if (j < 0 || j > m-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - FloatComplexMatrix q1 (m-1, m-1); - FloatComplexMatrix r1 (m-1, n); + OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); + OCTAVE_LOCAL_BUFFER (float, rw, m); + F77_XFCN (cqrder, CQRDER, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + w, rw)); - F77_XFCN (cqrder, CQRDER, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1 )); - - q = q1; - r = r1; + q.resize (m - 1, m - 1); + r.resize (m - 1, n); } } @@ -283,22 +418,19 @@ 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"); + (*current_liboctave_error_handler) ("qrshift: index out of range"); else - F77_XFCN (cqrshc, CQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (FloatComplex, w, k); + OCTAVE_LOCAL_BUFFER (float, rw, k); + F77_XFCN (cqrshc, CQRSHC, (m, n, k, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + i + 1, j + 1, w, rw)); + } } -void -FloatComplexQR::economize (void) -{ - octave_idx_type r_nc = r.columns (); - - if (r.rows () > r_nc) - { - q.resize (q.rows (), r_nc); - r.resize (r_nc, r_nc); - } -} +#endif /* ;;; Local Variables: *** diff -r 3d8a914c580e -r d66c9b6e506a liboctave/fCmplxQR.h --- a/liboctave/fCmplxQR.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/fCmplxQR.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -65,19 +66,27 @@ FloatComplexMatrix R (void) const { return r; } +#ifdef HAVE_QRUPDATE + + void update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v); + void update (const FloatComplexMatrix& u, const FloatComplexMatrix& v); - void insert_col (const FloatComplexMatrix& u, octave_idx_type j); + void insert_col (const FloatComplexColumnVector& u, octave_idx_type j); + + void insert_col (const FloatComplexMatrix& u, const Array& j); void delete_col (octave_idx_type j); - void insert_row (const FloatComplexMatrix& u, octave_idx_type j); + void delete_col (const Array& j); + + void insert_row (const FloatComplexRowVector& u, octave_idx_type j); void delete_row (octave_idx_type j); void shift_cols (octave_idx_type i, octave_idx_type j); - void economize(); +#endif friend std::ostream& operator << (std::ostream&, const FloatComplexQR&); diff -r 3d8a914c580e -r d66c9b6e506a liboctave/floatCHOL.cc --- a/liboctave/floatCHOL.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/floatCHOL.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #ifdef HAVE_CONFIG_H #include #endif @@ -52,23 +51,30 @@ float*, const octave_idx_type&, const float&, float&, float*, octave_idx_type*, octave_idx_type& F77_CHAR_ARG_LEN_DECL); +#ifdef HAVE_QRUPDATE + F77_RET_T - F77_FUNC (sch1up, SCH1UP) (const octave_idx_type&, float*, float*, float*); + F77_FUNC (sch1up, SCH1UP) (const octave_idx_type&, float*, const octave_idx_type&, + float*, float*); F77_RET_T - F77_FUNC (sch1dn, SCH1DN) (const octave_idx_type&, float*, float*, float*, + F77_FUNC (sch1dn, SCH1DN) (const octave_idx_type&, float*, const octave_idx_type&, + float*, float*, octave_idx_type&); + + F77_RET_T + F77_FUNC (schinx, SCHINX) (const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, float*, float*, octave_idx_type&); F77_RET_T - F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - float*, float*, const octave_idx_type&, const octave_idx_type&); + F77_FUNC (schdex, SCHDEX) (const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, float*); F77_RET_T - F77_FUNC (schinx, SCHINX) (const octave_idx_type&, const float*, float*, const octave_idx_type&, - const float*, octave_idx_type&); - - F77_RET_T - F77_FUNC (schdex, SCHDEX) (const octave_idx_type&, const float*, float*, const octave_idx_type&); + F77_FUNC (schshx, SCHSHX) (const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + float*); +#endif } octave_idx_type @@ -186,26 +192,28 @@ (*current_liboctave_error_handler) ("FloatCHOL requires square matrix"); } +#ifdef HAVE_QRUPDATE + void -FloatCHOL::update (const FloatMatrix& u) +FloatCHOL::update (const FloatColumnVector& u) { octave_idx_type n = chol_mat.rows (); if (u.length () == n) { - FloatMatrix tmp = u; + FloatColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w)); + F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w)); } else - (*current_liboctave_error_handler) ("FloatCHOL update dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); } octave_idx_type -FloatCHOL::downdate (const FloatMatrix& u) +FloatCHOL::downdate (const FloatColumnVector& u) { octave_idx_type info = -1; @@ -213,38 +221,40 @@ if (u.length () == n) { - FloatMatrix tmp = u; + FloatColumnVector utmp = u; OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), - tmp.fortran_vec (), w, info)); + F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (), + utmp.fortran_vec (), w, info)); } else - (*current_liboctave_error_handler) ("FloatCHOL downdate dimension mismatch"); + (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); return info; } octave_idx_type -FloatCHOL::insert_sym (const FloatMatrix& u, octave_idx_type j) +FloatCHOL::insert_sym (const FloatColumnVector& 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) ("FloatCHOL insert dimension mismatch"); + if (u.length () != n + 1) + (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("FloatCHOL insert index out of range"); + (*current_liboctave_error_handler) ("cholinsert: index out of range"); else { - FloatMatrix chol_mat1 (n+1, n+1); + FloatColumnVector utmp = u; + + OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (schinx, SCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), - j+1, u.data (), info)); + chol_mat.resize (n+1, n+1); - chol_mat = chol_mat1; + F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, utmp.fortran_vec (), w, info)); } return info; @@ -256,14 +266,15 @@ octave_idx_type n = chol_mat.rows (); if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("FloatCHOL delete index out of range"); + (*current_liboctave_error_handler) ("choldelete: index out of range"); else { - FloatMatrix chol_mat1 (n-1, n-1); + OCTAVE_LOCAL_BUFFER (float, w, n); - F77_XFCN (schdex, SCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); + F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + j + 1, w)); - chol_mat = chol_mat1; + chol_mat.resize (n-1, n-1); } } @@ -271,14 +282,20 @@ FloatCHOL::shift_sym (octave_idx_type i, octave_idx_type j) { octave_idx_type n = chol_mat.rows (); - float dummy; if (i < 0 || i > n-1 || j < 0 || j > n-1) - (*current_liboctave_error_handler) ("FloatCHOL shift index out of range"); + (*current_liboctave_error_handler) ("cholshift: index out of range"); else - F77_XFCN (sqrshc, SQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (float, w, 2*n); + + F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), + i + 1, j + 1, w)); + } } +#endif + FloatMatrix chol2inv (const FloatMatrix& r) { diff -r 3d8a914c580e -r d66c9b6e506a liboctave/floatCHOL.h --- a/liboctave/floatCHOL.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/floatCHOL.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,14 +22,13 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #if !defined (octave_FloatCHOL_h) #define octave_FloatCHOL_h 1 #include #include "fMatrix.h" +#include "fColVector.h" class OCTAVE_API @@ -64,16 +64,20 @@ void set (const FloatMatrix& R); - void update (const FloatMatrix& u); +#ifdef HAVE_QRUPDATE + + void update (const FloatColumnVector& u); - octave_idx_type downdate (const FloatMatrix& u); + octave_idx_type downdate (const FloatColumnVector& u); - octave_idx_type insert_sym (const FloatMatrix& u, octave_idx_type j); + octave_idx_type insert_sym (const FloatColumnVector& u, octave_idx_type j); void delete_sym (octave_idx_type j); void shift_sym (octave_idx_type i, octave_idx_type j); +#endif + friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatCHOL& a); private: diff -r 3d8a914c580e -r d66c9b6e506a liboctave/floatQR.cc --- a/liboctave/floatQR.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/floatQR.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -30,6 +31,7 @@ #include "lo-error.h" #include "Range.h" #include "idx-vector.h" +#include "oct-locbuf.h" extern "C" { @@ -41,33 +43,40 @@ F77_FUNC (sorgqr, SORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, float*, const octave_idx_type&, float*, float*, const octave_idx_type&, octave_idx_type&); - // these come from qrupdate +#ifdef HAVE_QRUPDATE F77_RET_T F77_FUNC (sqr1up, SQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - float*, float*, const float*, const float*); + float*, const octave_idx_type&, float*, const octave_idx_type&, + float*, float*, float*); F77_RET_T F77_FUNC (sqrinc, SQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - float*, const float*, float*, const octave_idx_type&, const float*); + float*, const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, const float*, float*); F77_RET_T F77_FUNC (sqrdec, SQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - float*, const float*, float*, const octave_idx_type&); + float*, const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, float*); F77_RET_T F77_FUNC (sqrinr, SQRINR) (const octave_idx_type&, const octave_idx_type&, - const float*, float*, const float*, float*, - const octave_idx_type&, const float*); + float*, const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, const float*, float*); F77_RET_T F77_FUNC (sqrder, SQRDER) (const octave_idx_type&, const octave_idx_type&, - const float*, float*, const float*, float *, - const octave_idx_type&); + float*, const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, float*); F77_RET_T F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - float*, float*, const octave_idx_type&, const octave_idx_type&); + float*, const octave_idx_type&, float*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + float*); + +#endif } FloatQR::FloatQR (const FloatMatrix& a, QR::type qr_type) @@ -160,6 +169,26 @@ this->r = r_arg; } +#ifdef HAVE_QRUPDATE + +void +FloatQR::update (const FloatColumnVector& u, const FloatColumnVector& v) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () == m && v.length () == n) + { + FloatColumnVector utmp = u, vtmp = v; + OCTAVE_LOCAL_BUFFER (float, w, 2*k); + F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w)); + } + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); +} + void FloatQR::update (const FloatMatrix& u, const FloatMatrix& v) { @@ -167,32 +196,93 @@ octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); - if (u.length () == m && v.length () == n) - F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), - u.data (), v.data ())); + if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) + { + OCTAVE_LOCAL_BUFFER (float, w, 2*k); + for (octave_idx_type i = 0; i < u.cols (); i++) + { + FloatColumnVector utmp = u.column (i), vtmp = v.column (i); + F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, + utmp.fortran_vec (), vtmp.fortran_vec (), w)); + } + } else - (*current_liboctave_error_handler) ("QR update dimensions mismatch"); + (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); } void -FloatQR::insert_col (const FloatMatrix& u, octave_idx_type j) +FloatQR::insert_col (const FloatColumnVector& u, octave_idx_type j) { octave_idx_type m = q.rows (); octave_idx_type n = r.columns (); octave_idx_type k = q.columns (); if (u.length () != m) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > n) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - FloatMatrix r1 (m, n+1); + if (k < m) + { + q.resize (m, k+1); + r.resize (k+1, n+1); + } + else + { + r.resize (k, n+1); + } + + FloatColumnVector utmp = u; + OCTAVE_LOCAL_BUFFER (float, w, k); + F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + utmp.data (), w)); + } +} + +void +FloatQR::insert_col (const FloatMatrix& u, const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); - F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1, u.data ())); + Array jsi; + Array js = j.sort (jsi, ASCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); - r = r1; + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (u.length () != m || u.columns () != nj) + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); + else if (nj > 0 && (js(0) < 0 || js(nj-1) > n)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + octave_idx_type kmax = std::min (k + nj, m); + if (k < m) + { + q.resize (m, kmax); + r.resize (kmax, n + nj); + } + else + { + r.resize (k, n + nj); + } + + OCTAVE_LOCAL_BUFFER (float, w, kmax); + for (octave_idx_type i = 0; i < js.length (); i++) + { + FloatColumnVector utmp = u.column (jsi(i)); + F77_XFCN (sqrinc, SQRINC, (m, n + i, std::min (kmax, k + i), + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, + utmp.data (), w)); + } } } @@ -203,41 +293,87 @@ octave_idx_type k = r.rows (); octave_idx_type n = r.columns (); - if (k < m && k < n) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); - else if (j < 0 || j > n-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - FloatMatrix r1 (k, n-1); + OCTAVE_LOCAL_BUFFER (float, w, k); + F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, w)); - F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), r.data (), - r1.fortran_vec (), j+1)); - - r = r1; + if (k < m) + { + q.resize (m, k-1); + r.resize (k-1, n-1); + } + else + { + r.resize (k, n-1); + } } } void -FloatQR::insert_row (const FloatMatrix& u, octave_idx_type j) +FloatQR::delete_col (const Array& j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + Array jsi; + Array js = j.sort (jsi, DESCENDING); + octave_idx_type nj = js.length (); + bool dups = false; + for (octave_idx_type i = 0; i < nj - 1; i++) + dups = dups && js(i) == js(i+1); + + if (dups) + (*current_liboctave_error_handler) ("qrinsert: duplicate index detected"); + else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0)) + (*current_liboctave_error_handler) ("qrinsert: index out of range"); + else if (nj > 0) + { + OCTAVE_LOCAL_BUFFER (float, w, k); + for (octave_idx_type i = 0; i < js.length (); i++) + { + F77_XFCN (sqrdec, SQRDEC, (m, n - i, k == m ? k : k - i, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), js(i) + 1, w)); + } + if (k < m) + { + q.resize (m, k - nj); + r.resize (k - nj, n - nj); + } + else + { + r.resize (k, n - nj); + } + + } +} + +void +FloatQR::insert_row (const FloatRowVector& u, octave_idx_type j) { octave_idx_type m = r.rows (); octave_idx_type n = r.columns (); + octave_idx_type k = std::min (m, n); if (! q.is_square () || u.length () != n) - (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch"); else if (j < 0 || j > m) - (*current_liboctave_error_handler) ("QR insert index out of range"); + (*current_liboctave_error_handler) ("qrinsert: index out of range"); else { - FloatMatrix q1 (m+1, m+1); - FloatMatrix r1 (m+1, n); + q.resize (m + 1, m + 1); + r.resize (m + 1, n); + FloatRowVector utmp = u; + OCTAVE_LOCAL_BUFFER (float, w, k); + F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + j + 1, utmp.fortran_vec (), w)); - F77_XFCN (sqrinr, SQRINR, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1, u.data ())); - - q = q1; - r = r1; } } @@ -248,19 +384,18 @@ octave_idx_type n = r.columns (); if (! q.is_square ()) - (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch"); else if (j < 0 || j > m-1) - (*current_liboctave_error_handler) ("QR delete index out of range"); + (*current_liboctave_error_handler) ("qrdelete: index out of range"); else { - FloatMatrix q1 (m-1, m-1); - FloatMatrix r1 (m-1, n); + OCTAVE_LOCAL_BUFFER (float, w, 2*m); + F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), j + 1, + w)); - F77_XFCN (sqrder, SQRDER, (m, n, q.data (), q1.fortran_vec (), - r.data (), r1.fortran_vec (), j+1 )); - - q = q1; - r = r1; + q.resize (m - 1, m - 1); + r.resize (m - 1, n); } } @@ -272,22 +407,18 @@ 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"); + (*current_liboctave_error_handler) ("qrshift: index out of range"); else - F77_XFCN (sqrshc, SQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); + { + OCTAVE_LOCAL_BUFFER (float, w, 2*k); + F77_XFCN (sqrshc, SQRSHC, (m, n, k, + q.fortran_vec (), q.rows (), + r.fortran_vec (), r.rows (), + i + 1, j + 1, w)); + } } -void -FloatQR::economize (void) -{ - octave_idx_type r_nc = r.columns (); - - if (r.rows () > r_nc) - { - q.resize (q.rows (), r_nc); - r.resize (r_nc, r_nc); - } -} +#endif /* diff -r 3d8a914c580e -r d66c9b6e506a liboctave/floatQR.h --- a/liboctave/floatQR.h Tue Jan 20 21:15:17 2009 +0100 +++ b/liboctave/floatQR.h Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek This file is part of Octave. @@ -21,8 +22,6 @@ */ -// updating/downdating by Jaroslav Hajek 2008 - #if !defined (octave_FloatQR_h) #define octave_FloatQR_h 1 @@ -65,19 +64,27 @@ FloatMatrix R (void) const { return r; } +#ifdef HAVE_QRUPDATE + + void update (const FloatColumnVector& u, const FloatColumnVector& v); + void update (const FloatMatrix& u, const FloatMatrix& v); - void insert_col (const FloatMatrix& u, octave_idx_type j); + void insert_col (const FloatColumnVector& u, octave_idx_type j); + + void insert_col (const FloatMatrix& u, const Array& j); void delete_col (octave_idx_type j); - void insert_row (const FloatMatrix& u, octave_idx_type j); + void delete_col (const Array& j); + + void insert_row (const FloatRowVector& u, octave_idx_type j); void delete_row (octave_idx_type j); void shift_cols (octave_idx_type i, octave_idx_type j); - void economize (void); +#endif friend std::ostream& operator << (std::ostream&, const FloatQR&); diff -r 3d8a914c580e -r d66c9b6e506a src/ChangeLog --- a/src/ChangeLog Tue Jan 20 21:15:17 2009 +0100 +++ b/src/ChangeLog Tue Jan 20 21:16:42 2009 +0100 @@ -1,3 +1,10 @@ +2009-01-17 Jaroslav Hajek + + * DLD-FUNCTIONS/qr.cc (Fqrupdate, Fqrinsert, Fqrdelete, Fqrshift): + Reflect changes in liboctave. + * DLD-FUNCTIONS/chol.cc (Fcholupdate, Fcholinsert): + Reflect changes in liboctave. + 2009-01-19 Jaroslav Hajek * ov.h (octave_value::make_unique (int)): New method. diff -r 3d8a914c580e -r d66c9b6e506a src/DLD-FUNCTIONS/chol.cc --- a/src/DLD-FUNCTIONS/chol.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/src/DLD-FUNCTIONS/chol.cc Tue Jan 20 21:16:42 2009 +0100 @@ -2,6 +2,8 @@ Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek +Copyright (C) 2008, 2009 VZLU Prague This file is part of Octave. @@ -21,9 +23,6 @@ */ -// The cholupdate, cholinsert, choldelete and cholshift functions were -// written by Jaroslav Hajek , Copyright (C) 2008 -// VZLU Prague, a.s., Czech Republic. #ifdef HAVE_CONFIG_H #include @@ -577,6 +576,8 @@ return retval; } +#ifdef HAVE_QRUPDATE + DEFUN_DLD (cholupdate, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\ @@ -628,100 +629,84 @@ if (down || op == "+") if (argr.columns () == n && argu.rows () == n && argu.columns () == 1) { + int err = 0; if (argr.is_single_type () || argu.is_single_type ()) { - if (argr.is_real_matrix () && argu.is_real_matrix ()) + if (argr.is_real_type () && argu.is_real_type ()) { // real case FloatMatrix R = argr.float_matrix_value (); - FloatMatrix u = argu.float_matrix_value (); + FloatColumnVector u = argu.float_column_vector_value (); FloatCHOL fact; fact.set (R); - int err = 0; if (down) err = fact.downdate (u); else fact.update (u); - if (nargout > 1) - retval(1) = err; - else if (err) - error ("cholupdate: downdate violates positiveness"); - retval(0) = fact.chol_matrix (); } else { // complex case FloatComplexMatrix R = argr.float_complex_matrix_value (); - FloatComplexMatrix u = argu.float_complex_matrix_value (); + FloatComplexColumnVector u = argu.float_complex_column_vector_value (); FloatComplexCHOL fact; fact.set (R); - int err = 0; if (down) err = fact.downdate (u); else fact.update (u); - if (nargout > 1) - retval(1) = err; - else if (err) - error ("cholupdate: downdate violates positiveness"); - retval(0) = fact.chol_matrix (); } } else { - if (argr.is_real_matrix () && argu.is_real_matrix ()) + if (argr.is_real_type () && argu.is_real_type ()) { // real case Matrix R = argr.matrix_value (); - Matrix u = argu.matrix_value (); + ColumnVector u = argu.column_vector_value (); CHOL fact; fact.set (R); - int err = 0; if (down) err = fact.downdate (u); else fact.update (u); - if (nargout > 1) - retval(1) = err; - else if (err) - error ("cholupdate: downdate violates positiveness"); - retval(0) = fact.chol_matrix (); } else { // complex case ComplexMatrix R = argr.complex_matrix_value (); - ComplexMatrix u = argu.complex_matrix_value (); + ComplexColumnVector u = argu.complex_column_vector_value (); ComplexCHOL fact; fact.set (R); - int err = 0; if (down) err = fact.downdate (u); else fact.update (u); - if (nargout > 1) - retval(1) = err; - else if (err) - error ("cholupdate: downdate violates positiveness"); - retval(0) = fact.chol_matrix (); } } + + if (nargout > 1) + retval(1) = err; + else if (err == 1) + error ("cholupdate: downdate violates positiveness"); + else if (err == 2) + error ("cholupdate: singular matrix"); } else error ("cholupdate: dimension mismatch"); @@ -853,22 +838,18 @@ { if (j > 0 && j <= n+1) { + int err = 0; if (argr.is_single_type () || argu.is_single_type ()) { - if (argr.is_real_matrix () && argu.is_real_matrix ()) + if (argr.is_real_type () && argu.is_real_type ()) { // real case FloatMatrix R = argr.float_matrix_value (); - FloatMatrix u = argu.float_matrix_value (); + FloatColumnVector u = argu.float_column_vector_value (); FloatCHOL 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"); + err = fact.insert_sym (u, j-1); retval(0) = fact.chol_matrix (); } @@ -876,36 +857,26 @@ { // complex case FloatComplexMatrix R = argr.float_complex_matrix_value (); - FloatComplexMatrix u = argu.float_complex_matrix_value (); + FloatComplexColumnVector u = argu.float_complex_column_vector_value (); FloatComplexCHOL 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"); + err = fact.insert_sym (u, j-1); retval(0) = fact.chol_matrix (); } } else { - if (argr.is_real_matrix () && argu.is_real_matrix ()) + if (argr.is_real_type () && argu.is_real_type ()) { // real case Matrix R = argr.matrix_value (); - Matrix u = argu.matrix_value (); + ColumnVector u = argu.column_vector_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"); + err = fact.insert_sym (u, j-1); retval(0) = fact.chol_matrix (); } @@ -913,20 +884,24 @@ { // complex case ComplexMatrix R = argr.complex_matrix_value (); - ComplexMatrix u = argu.complex_matrix_value (); + ComplexColumnVector u = argu.complex_column_vector_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"); + err = fact.insert_sym (u, j-1); retval(0) = fact.chol_matrix (); } } + + if (nargout > 1) + retval(1) = err; + else if (err == 1) + error ("cholinsert: insertion violates positiveness"); + else if (err == 2) + error ("cholinsert: singular matrix"); + else if (err == 3) + error ("cholinsert: diagonal element must be real"); } else error ("cholinsert: index out of range"); @@ -1037,7 +1012,7 @@ { if (argr.is_single_type ()) { - if (argr.is_real_matrix ()) + if (argr.is_real_type ()) { // real case FloatMatrix R = argr.float_matrix_value (); @@ -1062,7 +1037,7 @@ } else { - if (argr.is_real_matrix ()) + if (argr.is_real_type ()) { // real case Matrix R = argr.matrix_value (); @@ -1178,7 +1153,7 @@ if (argr.is_single_type () && argi.is_single_type () && argj.is_single_type ()) { - if (argr.is_real_matrix ()) + if (argr.is_real_type ()) { // real case FloatMatrix R = argr.float_matrix_value (); @@ -1203,7 +1178,7 @@ } else { - if (argr.is_real_matrix ()) + if (argr.is_real_type ()) { // real case Matrix R = argr.matrix_value (); @@ -1301,6 +1276,8 @@ %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single')) */ +#endif + /* ;;; Local Variables: *** ;;; mode: C++ *** diff -r 3d8a914c580e -r d66c9b6e506a src/DLD-FUNCTIONS/qr.cc --- a/src/DLD-FUNCTIONS/qr.cc Tue Jan 20 21:15:17 2009 +0100 +++ b/src/DLD-FUNCTIONS/qr.cc Tue Jan 20 21:16:42 2009 +0100 @@ -1,6 +1,8 @@ /* Copyright (C) 1996, 1997, 1999, 2000, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008, 2009 Jaroslav Hajek +Copyright (C) 2008, 2009 VZLU Prague This file is part of Octave. @@ -20,10 +22,6 @@ */ -// The qrupdate, qrinsert, qrdelete and qrshift functions were written by -// Jaroslav Hajek , Copyright (C) 2008 VZLU -// Prague, a.s., Czech Republic. - #ifdef HAVE_CONFIG_H #include #endif @@ -741,6 +739,24 @@ */ +#ifdef HAVE_QRUPDATE + +static +bool check_qr_dims (const octave_value& q, const octave_value& r, + bool allow_ecf = false) +{ + octave_idx_type m = q.rows (), k = r.rows (), n = r.columns (); + return ((q.ndims () == 2 || r.ndims () == 2 && k == q.columns ()) + && (m == k || (allow_ecf && k == n && k < m))); +} + +static +bool check_index (const octave_value& i, bool vector_allowed = false) +{ + return ((i.is_real_type () || i.is_integer_type ()) + && (i.is_scalar_type () || vector_allowed)); +} + DEFUN_DLD (qrupdate, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{Q1}, @var{R1}] =} qrupdate (@var{Q}, @var{R}, @var{u}, @var{v})\n\ @@ -748,10 +764,14 @@ @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\ @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization\n\ of @w{@var{A} + @var{u}*@var{v}'}, where @var{u} and @var{v} are\n\ -column vectors (rank-1 update).\n\ +column vectors (rank-1 update) or matrices with equal number of columns\n\ +(rank-k update). Notice that the latter case is done as a sequence of rank-1 updates;\n\ +thus, for k large enough, it will be both faster and more accurate to recompute\n\ +the factorization from scratch.\n\ \n\ -If the matrix @var{Q} is not square, the matrix @var{A} is updated by\n\ -Q*Q'*u*v' instead of u*v'.\n\ +The QR factorization supplied may be either full\n\ +(Q is square) or economized (R is square).\n\ +\n\ @seealso{qr, qrinsert, qrdelete}\n\ @end deftypefn") { @@ -772,18 +792,12 @@ if (argq.is_numeric_type () && argr.is_numeric_type () && argu.is_numeric_type () && argv.is_numeric_type ()) { - octave_idx_type m = argq.rows (); - octave_idx_type n = argr.columns (); - octave_idx_type k = argq.columns (); - - if (argr.rows () == k - && argu.rows () == m && argu.columns () == 1 - && argv.rows () == n && argv.columns () == 1) + if (check_qr_dims (argq, argr, true)) { - if (argq.is_real_matrix () - && argr.is_real_matrix () - && argu.is_real_matrix () - && argv.is_real_matrix ()) + if (argq.is_real_type () + && argr.is_real_type () + && argu.is_real_type () + && argv.is_real_type ()) { // all real case if (argq.is_single_type () @@ -935,12 +949,19 @@ @code{\"row\"}).\n\ \n\ The default value of @var{orient} is @code{\"col\"}.\n\ +If @var{orient} is @code{\"col\"},\n\ +@var{u} may be a matrix and @var{j} an index vector\n\ +resulting in the QR@tie{}factorization of a matrix @var{B} such that\n\ +@w{B(:,@var{j})} gives @var{u} and @w{B(:,@var{j}) = []} gives @var{A}.\n\ +Notice that the latter case is done as a sequence of k insertions;\n\ +thus, for k large enough, it will be both faster and more accurate to recompute\n\ +the factorization from scratch.\n\ \n\ -If @var{orient} is @code{\"col\"} and the matrix @var{Q} is not square,\n\ -then what gets inserted is the projection of @var{u} onto the space\n\ -spanned by columns of @var{Q}, i.e. Q*Q'*u.\n\ +If @var{orient} is @code{\"col\"},\n\ +the QR factorization supplied may be either full\n\ +(Q is square) or economized (R is square).\n\ \n\ -If @var{orient} is @code{\"row\"}, @var{Q} must be square.\n\ +If @var{orient} is @code{\"row\"}, full factorization is needed.\n\ @seealso{qr, qrupdate, qrdelete}\n\ @end deftypefn") { @@ -959,30 +980,24 @@ octave_value argx = args(3); if (argq.is_numeric_type () && argr.is_numeric_type () - && argj.is_scalar_type () && argx.is_numeric_type () + && argx.is_numeric_type () && (nargin < 5 || args(4).is_string ())) { - octave_idx_type m = argq.rows (); - octave_idx_type n = argr.columns (); - octave_idx_type k = argq.columns (); - std::string orient = (nargin < 5) ? "col" : args(4).string_value (); - bool row = orient == "row"; + bool col = orient == "col"; - if (row || orient == "col") - if (argr.rows () == k - && (! row || m == k) - && argx.rows () == (row ? 1 : m) - && argx.columns () == (row ? n : 1)) + if (col || orient == "row") + if (check_qr_dims (argq, argr, col) + && (col || argx.rows () == 1)) { - octave_idx_type j = argj.idx_type_value (); - - if (j >= 1 && j <= (row ? n : m)+1) + if (check_index (argj, col)) { - if (argq.is_real_matrix () - && argr.is_real_matrix () - && argx.is_real_matrix ()) + MArray j = argj.int_vector_value (); + + if (argq.is_real_type () + && argr.is_real_type () + && argx.is_real_type ()) { // real case if (argq.is_single_type () @@ -995,10 +1010,10 @@ FloatQR fact (Q, R); - if (row) - fact.insert_row (x, j-1); + if (col) + fact.insert_col (x, j-1); else - fact.insert_col (x, j-1); + fact.insert_row (x.row (0), j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1012,10 +1027,10 @@ QR fact (Q, R); - if (row) - fact.insert_row (x, j-1); + if (col) + fact.insert_col (x, j-1); else - fact.insert_col (x, j-1); + fact.insert_row (x.row (0), j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1035,10 +1050,10 @@ FloatComplexQR fact (Q, R); - if (row) - fact.insert_row (x, j-1); + if (col) + fact.insert_col (x, j-1); else - fact.insert_col (x, j-1); + fact.insert_row (x.row (0), j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1051,10 +1066,10 @@ ComplexQR fact (Q, R); - if (row) - fact.insert_row (x, j-1); + if (col) + fact.insert_col (x, j-1); else - fact.insert_col (x, j-1); + fact.insert_row (x.row (0), j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1063,7 +1078,7 @@ } else - error ("qrinsert: index j out of range"); + error ("qrinsert: invalid index"); } else error ("qrinsert: dimension mismatch"); @@ -1150,18 +1165,19 @@ \n\ The default value of @var{orient} is \"col\".\n\ \n\ -If @var{orient} is \"col\", the matrix @var{Q} is not required to\n\ -be square.\n\ +If @var{orient} is @code{\"col\"},\n\ +@var{j} may be an index vector\n\ +resulting in the QR@tie{}factorization of a matrix @var{B} such that\n\ +@w{A(:,@var{j}) = []} gives @var{B}.\n\ +Notice that the latter case is done as a sequence of k deletions;\n\ +thus, for k large enough, it will be both faster and more accurate to recompute\n\ +the factorization from scratch.\n\ \n\ -For @sc{Matlab} compatibility, if @var{Q} is nonsquare on input, the\n\ -updated factorization is always stripped to the economical form, i.e.\n\ -@code{columns (Q) == rows (R) <= columns (R)}.\n\ +If @var{orient} is @code{\"col\"},\n\ +the QR factorization supplied may be either full\n\ +(Q is square) or economized (R is square).\n\ \n\ -To get the less intelligent but more natural behaviour when @var{Q}\n\ -retains it shape and @var{R} loses one column, set @var{orient} to\n\ -\"col+\" instead.\n\ -\n\ -If @var{orient} is \"row\", @var{Q} must be square.\n\ +If @var{orient} is @code{\"row\"}, full factorization is needed.\n\ @seealso{qr, qrinsert, qrupdate}\n\ @end deftypefn") { @@ -1179,27 +1195,21 @@ octave_value argj = args(2); if (argq.is_numeric_type () && argr.is_numeric_type () - && argj.is_scalar_type () && (nargin < 4 || args(3).is_string ())) { - octave_idx_type m = argq.rows (); - octave_idx_type k = argq.columns (); - octave_idx_type n = argr.columns (); - std::string orient = (nargin < 4) ? "col" : args(3).string_value (); - bool row = orient == "row"; - bool colp = orient == "col+"; + bool col = orient == "col"; - if (row || colp || orient == "col") - if (argr.rows () == k - && (! row || m == k)) + if (col || orient == "row") + if (check_qr_dims (argq, argr, col)) { - octave_idx_type j = argj.scalar_value (); - if (j >= 1 && j <= (row ? n : m)) + if (check_index (argj, col)) { - if (argq.is_real_matrix () - && argr.is_real_matrix ()) + MArray j = argj.int_vector_value (); + + if (argq.is_real_type () + && argr.is_real_type ()) { // real case if (argq.is_single_type () @@ -1210,15 +1220,10 @@ FloatQR fact (Q, R); - if (row) - fact.delete_row (j-1); + if (col) + fact.delete_col (j-1); else - { - fact.delete_col (j-1); - - if (! colp && k < m) - fact.economize (); - } + fact.delete_row (j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1230,15 +1235,10 @@ QR fact (Q, R); - if (row) - fact.delete_row (j-1); + if (col) + fact.delete_col (j-1); else - { - fact.delete_col (j-1); - - if (! colp && k < m) - fact.economize (); - } + fact.delete_row (j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1255,15 +1255,10 @@ FloatComplexQR fact (Q, R); - if (row) - fact.delete_row (j-1); + if (col) + fact.delete_col (j-1); else - { - fact.delete_col (j-1); - - if (! colp && k < m) - fact.economize (); - } + fact.delete_row (j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1275,15 +1270,10 @@ ComplexQR fact (Q, R); - if (row) - fact.delete_row (j-1); + if (col) + fact.delete_col (j-1); else - { - fact.delete_col (j-1); - - if (! colp && k < m) - fact.economize (); - } + fact.delete_row (j(0)-1); retval(1) = fact.R (); retval(0) = fact.Q (); @@ -1291,7 +1281,7 @@ } } else - error ("qrdelete: index j out of range"); + error ("qrdelete: invalid index"); } else error ("qrdelete: dimension mismatch"); @@ -1439,20 +1429,17 @@ octave_value argi = args(2); octave_value argj = args(3); - if (argq.is_numeric_type () && argr.is_numeric_type () - && argi.is_real_scalar () && argj.is_real_scalar ()) + if (argq.is_numeric_type () && argr.is_numeric_type ()) { - octave_idx_type n = argr.columns (); - octave_idx_type k = argq.columns (); + if (check_qr_dims (argq, argr, true)) + { + if (check_index (argi) && check_index (argj)) + { + octave_idx_type i = argi.int_value (); + octave_idx_type j = argj.int_value (); - 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 ()) + if (argq.is_real_type () + && argr.is_real_type ()) { // all real case if (argq.is_single_type () @@ -1508,7 +1495,7 @@ } } else - error ("qrshift: index out of range"); + error ("qrshift: invalid index"); } else error ("qrshift: dimensions mismatch"); @@ -1593,6 +1580,8 @@ %! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single')) */ +#endif + /* ;;; Local Variables: *** ;;; mode: C++ *** diff -r 3d8a914c580e -r d66c9b6e506a src/Makefile.in --- a/src/Makefile.in Tue Jan 20 21:15:17 2009 +0100 +++ b/src/Makefile.in Tue Jan 20 21:16:42 2009 +0100 @@ -301,7 +301,7 @@ -L../libcruft $(LIBCRUFT) -L../liboctave $(LIBOCTAVE) \ -L. $(LIBOCTINTERP) $(CHOLMOD_LIBS) $(UMFPACK_LIBS) $(AMD_LIBS) \ $(CAMD_LIBS) $(COLAMD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \ - $(FFTW_LIBS) $(ARPACK_LIBS) $(LIBS) $(FLIBS) + $(FFTW_LIBS) $(QRUPDATE_LIBS) $(ARPACK_LIBS) $(LIBS) $(FLIBS) BUILT_DISTFILES = DOCSTRINGS oct-gperf.h parse.cc lex.cc y.tab.h \ $(OPT_HANDLERS) $(BUILT_EXTRAS) @@ -372,7 +372,8 @@ $(OCTAVE_LIBS) \ $(LEXLIB) $(UMFPACK_LIBS) $(AMD_LIBS) $(CAMD_LIBS) $(COLAMD_LIBS) \ $(CHOLMOD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \ - $(FFTW_LIBS) $(ARPACK_LIBS) $(OPENGL_LIBS) $(LIBS) $(FLIBS) + $(FFTW_LIBS) $(QRUPDATE_LIBS) $(ARPACK_LIBS) $(OPENGL_LIBS) \ + $(LIBS) $(FLIBS) stmp-pic: pic @if [ -f stmp-pic ]; then \ @@ -652,6 +653,8 @@ __delaunayn__.oct: OCT_LINK_DEPS += $(QHULL_LIBS) __voronoi__.oct: OCT_LINK_DEPS += $(QHULL_LIBS) eigs.oct: OCT_LINK_DEPS += $(ARPACK_LIBS) +qr.oct: OCT_LINK_DEPS += $(QRUPDATE_LIBS) +chol.oct: OCT_LINK_DEPS += $(QRUPDATE_LIBS) regexp.oct: OCT_LINK_DEPS += $(REGEX_LIBS) urlwrite.oct: OCT_LINK_DEPS += $(CURL_LIBS) __glpk__.oct: OCT_LINK_DEPS += $(GLPK_LIBS)