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 cqrqh(m,n,R,ldr,c,s) 22c purpose: brings an upper trapezoidal matrix R into upper 23c Hessenberg form using min(m-1,n) Givens rotations. 24c (complex version) 25c arguments: 26c m (in) number of rows of the matrix R 27c n (in) number of columns of the matrix R 28c R (io) on entry, the upper Hessenberg matrix R 29c on exit, the updated upper trapezoidal matrix 30c ldr (in) leading dimension of R, >= m 31c c(in) rotation cosines, size at least min(m-1,n) 32c s(in) rotation sines, size at least min(m-1,n) 33c 34 integer m,n,ldr 35 complex R(ldr,*),s(*) 36 real c(*) 37 external xerbla 38 complex t 39 integer info,i,ii,j 40c quick return if possible. 41 if (m == 0 .or. m == 1 .or. n == 0) return 42c check arguments. 43 info = 0 44 if (m < 0) then 45 info = 1 46 else if (n < 0) then 47 info = 2 48 else if (ldr < m) then 49 info = 4 50 end if 51 if (info /= 0) then 52 call xerbla('CQRQH',info) 53 return 54 end if 55 do i = 1,n 56c apply stored rotations, column-wise 57 ii = min(m-1,i) 58 t = R(ii+1,i) 59 do j = ii,1,-1 60 R(j+1,i) = c(j)*t - conjg(s(j))*R(j,i) 61 t = c(j)*R(j,i) + s(j)*t 62 end do 63 R(1,i) = t 64 end do 65 end subroutine 66