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