1c Copyright (C) 2008, 2009 VZLU Prague, a.s., Czech Republic 2c 3c Author: Jaroslav Hajek <highegg@gmail.com> 4c 5c This file is part of qrupdate. 6c 7c qrupdate is free software; you can redistribute it and/or modify 8c it under the terms of the GNU General Public License as published by 9c the Free Software Foundation; either version 3 of the License, or 10c (at your option) any later version. 11c 12c This program is distributed in the hope that it will be useful, 13c but WITHOUT ANY WARRANTY; without even the implied warranty of 14c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15c GNU General Public License for more details. 16c 17c You should have received a copy of the GNU General Public License 18c along with this software; see the file COPYING. If not, see 19c <http://www.gnu.org/licenses/>. 20c 21 subroutine cqrder(m,n,Q,ldq,R,ldr,j,w,rw) 22c purpose: updates a QR factorization after deleting a row. 23c i.e., given an m-by-m unitary matrix Q, an m-by-n 24c upper trapezoidal matrix R and index j in the range 25c 1:m, this subroutine updates Q ->Q1 and an R -> R1 26c so that Q1 is again unitary, R1 upper trapezoidal, 27c and Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R. 28c (complex version) 29c 30c arguments: 31c m (in) number of rows of the matrix Q. 32c n (in) number of columns of the matrix R. 33c Q (io) on entry, the unitary matrix Q. 34c on exit, the updated matrix Q1. 35c ldq (in) leading dimension of Q. ldq >= m. 36c R (io) on entry, the original matrix R. 37c on exit, the updated matrix R1. 38c ldr (in) leading dimension of R. ldr >= m. 39c j (in) the position of the deleted row. 40c w (out) a workspace vector of size m. 41c rw (out) a real workspace vector of size m. 42c 43 integer m,n,j,ldq,ldr 44 complex Q(ldq,*),R(ldr,*),w(*) 45 real rw(*) 46 external xerbla,ccopy,cqrtv1,cqrot,cqrqh 47 integer info,i,k 48c quick return if possible 49 if (m == 1) return 50c check arguments 51 info = 0 52 if (m < 1) then 53 info = 1 54 else if (j < 1 .or. j > m) then 55 info = 7 56 end if 57 if (info /= 0) then 58 call xerbla('CQRDER',info) 59 return 60 end if 61c eliminate Q(j,2:m). 62 do k = 1,m 63 w(k) = conjg(Q(j,k)) 64 end do 65 call cqrtv1(m,w,rw) 66c apply rotations to Q. 67 call cqrot('B',m,m,Q,ldq,rw,w(2)) 68c form Q1. 69 do k = 1,m-1 70 if (j > 1) call ccopy(j-1,Q(1,k+1),1,Q(1,k),1) 71 if (j < m) call ccopy(m-j,Q(j+1,k+1),1,Q(j,k),1) 72 end do 73c apply rotations to R. 74 call cqrqh(m,n,R,ldr,rw,w(2)) 75c form R1. 76 do k = 1,n 77 do i = 1,m-1 78 R(i,k) = R(i+1,k) 79 end do 80 end do 81 end subroutine 82