1*> \brief \b ZLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLAHR2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahr2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahr2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahr2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) 22* 23* .. Scalar Arguments .. 24* INTEGER K, LDA, LDT, LDY, N, NB 25* .. 26* .. Array Arguments .. 27* COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ), 28* $ Y( LDY, NB ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) 38*> matrix A so that elements below the k-th subdiagonal are zero. The 39*> reduction is performed by an unitary similarity transformation 40*> Q**H * A * Q. The routine returns the matrices V and T which determine 41*> Q as a block reflector I - V*T*V**H, and also the matrix Y = A * V * T. 42*> 43*> This is an auxiliary routine called by ZGEHRD. 44*> \endverbatim 45* 46* Arguments: 47* ========== 48* 49*> \param[in] N 50*> \verbatim 51*> N is INTEGER 52*> The order of the matrix A. 53*> \endverbatim 54*> 55*> \param[in] K 56*> \verbatim 57*> K is INTEGER 58*> The offset for the reduction. Elements below the k-th 59*> subdiagonal in the first NB columns are reduced to zero. 60*> K < N. 61*> \endverbatim 62*> 63*> \param[in] NB 64*> \verbatim 65*> NB is INTEGER 66*> The number of columns to be reduced. 67*> \endverbatim 68*> 69*> \param[in,out] A 70*> \verbatim 71*> A is COMPLEX*16 array, dimension (LDA,N-K+1) 72*> On entry, the n-by-(n-k+1) general matrix A. 73*> On exit, the elements on and above the k-th subdiagonal in 74*> the first NB columns are overwritten with the corresponding 75*> elements of the reduced matrix; the elements below the k-th 76*> subdiagonal, with the array TAU, represent the matrix Q as a 77*> product of elementary reflectors. The other columns of A are 78*> unchanged. See Further Details. 79*> \endverbatim 80*> 81*> \param[in] LDA 82*> \verbatim 83*> LDA is INTEGER 84*> The leading dimension of the array A. LDA >= max(1,N). 85*> \endverbatim 86*> 87*> \param[out] TAU 88*> \verbatim 89*> TAU is COMPLEX*16 array, dimension (NB) 90*> The scalar factors of the elementary reflectors. See Further 91*> Details. 92*> \endverbatim 93*> 94*> \param[out] T 95*> \verbatim 96*> T is COMPLEX*16 array, dimension (LDT,NB) 97*> The upper triangular matrix T. 98*> \endverbatim 99*> 100*> \param[in] LDT 101*> \verbatim 102*> LDT is INTEGER 103*> The leading dimension of the array T. LDT >= NB. 104*> \endverbatim 105*> 106*> \param[out] Y 107*> \verbatim 108*> Y is COMPLEX*16 array, dimension (LDY,NB) 109*> The n-by-nb matrix Y. 110*> \endverbatim 111*> 112*> \param[in] LDY 113*> \verbatim 114*> LDY is INTEGER 115*> The leading dimension of the array Y. LDY >= N. 116*> \endverbatim 117* 118* Authors: 119* ======== 120* 121*> \author Univ. of Tennessee 122*> \author Univ. of California Berkeley 123*> \author Univ. of Colorado Denver 124*> \author NAG Ltd. 125* 126*> \ingroup complex16OTHERauxiliary 127* 128*> \par Further Details: 129* ===================== 130*> 131*> \verbatim 132*> 133*> The matrix Q is represented as a product of nb elementary reflectors 134*> 135*> Q = H(1) H(2) . . . H(nb). 136*> 137*> Each H(i) has the form 138*> 139*> H(i) = I - tau * v * v**H 140*> 141*> where tau is a complex scalar, and v is a complex vector with 142*> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in 143*> A(i+k+1:n,i), and tau in TAU(i). 144*> 145*> The elements of the vectors v together form the (n-k+1)-by-nb matrix 146*> V which is needed, with T and Y, to apply the transformation to the 147*> unreduced part of the matrix, using an update of the form: 148*> A := (I - V*T*V**H) * (A - Y*V**H). 149*> 150*> The contents of A on exit are illustrated by the following example 151*> with n = 7, k = 3 and nb = 2: 152*> 153*> ( a a a a a ) 154*> ( a a a a a ) 155*> ( a a a a a ) 156*> ( h h a a a ) 157*> ( v1 h a a a ) 158*> ( v1 v2 a a a ) 159*> ( v1 v2 a a a ) 160*> 161*> where a denotes an element of the original matrix A, h denotes a 162*> modified element of the upper Hessenberg matrix H, and vi denotes an 163*> element of the vector defining H(i). 164*> 165*> This subroutine is a slight modification of LAPACK-3.0's ZLAHRD 166*> incorporating improvements proposed by Quintana-Orti and Van de 167*> Gejin. Note that the entries of A(1:K,2:NB) differ from those 168*> returned by the original LAPACK-3.0's ZLAHRD routine. (This 169*> subroutine is not backward compatible with LAPACK-3.0's ZLAHRD.) 170*> \endverbatim 171* 172*> \par References: 173* ================ 174*> 175*> Gregorio Quintana-Orti and Robert van de Geijn, "Improving the 176*> performance of reduction to Hessenberg form," ACM Transactions on 177*> Mathematical Software, 32(2):180-194, June 2006. 178*> 179* ===================================================================== 180 SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) 181* 182* -- LAPACK auxiliary routine -- 183* -- LAPACK is a software package provided by Univ. of Tennessee, -- 184* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 185* 186* .. Scalar Arguments .. 187 INTEGER K, LDA, LDT, LDY, N, NB 188* .. 189* .. Array Arguments .. 190 COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ), 191 $ Y( LDY, NB ) 192* .. 193* 194* ===================================================================== 195* 196* .. Parameters .. 197 COMPLEX*16 ZERO, ONE 198 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 199 $ ONE = ( 1.0D+0, 0.0D+0 ) ) 200* .. 201* .. Local Scalars .. 202 INTEGER I 203 COMPLEX*16 EI 204* .. 205* .. External Subroutines .. 206 EXTERNAL ZAXPY, ZCOPY, ZGEMM, ZGEMV, ZLACPY, 207 $ ZLARFG, ZSCAL, ZTRMM, ZTRMV, ZLACGV 208* .. 209* .. Intrinsic Functions .. 210 INTRINSIC MIN 211* .. 212* .. Executable Statements .. 213* 214* Quick return if possible 215* 216 IF( N.LE.1 ) 217 $ RETURN 218* 219 DO 10 I = 1, NB 220 IF( I.GT.1 ) THEN 221* 222* Update A(K+1:N,I) 223* 224* Update I-th column of A - Y * V**H 225* 226 CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) 227 CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY, 228 $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 ) 229 CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) 230* 231* Apply I - V * T**H * V**H to this column (call it b) from the 232* left, using the last column of T as workspace 233* 234* Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) 235* ( V2 ) ( b2 ) 236* 237* where V1 is unit lower triangular 238* 239* w := V1**H * b1 240* 241 CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) 242 CALL ZTRMV( 'Lower', 'Conjugate transpose', 'UNIT', 243 $ I-1, A( K+1, 1 ), 244 $ LDA, T( 1, NB ), 1 ) 245* 246* w := w + V2**H * b2 247* 248 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, 249 $ ONE, A( K+I, 1 ), 250 $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) 251* 252* w := T**H * w 253* 254 CALL ZTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT', 255 $ I-1, T, LDT, 256 $ T( 1, NB ), 1 ) 257* 258* b2 := b2 - V2*w 259* 260 CALL ZGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, 261 $ A( K+I, 1 ), 262 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) 263* 264* b1 := b1 - V1*w 265* 266 CALL ZTRMV( 'Lower', 'NO TRANSPOSE', 267 $ 'UNIT', I-1, 268 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) 269 CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) 270* 271 A( K+I-1, I-1 ) = EI 272 END IF 273* 274* Generate the elementary reflector H(I) to annihilate 275* A(K+I+1:N,I) 276* 277 CALL ZLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, 278 $ TAU( I ) ) 279 EI = A( K+I, I ) 280 A( K+I, I ) = ONE 281* 282* Compute Y(K+1:N,I) 283* 284 CALL ZGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, 285 $ ONE, A( K+1, I+1 ), 286 $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 ) 287 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, 288 $ ONE, A( K+I, 1 ), LDA, 289 $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) 290 CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, 291 $ Y( K+1, 1 ), LDY, 292 $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 ) 293 CALL ZSCAL( N-K, TAU( I ), Y( K+1, I ), 1 ) 294* 295* Compute T(1:I,I) 296* 297 CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) 298 CALL ZTRMV( 'Upper', 'No Transpose', 'NON-UNIT', 299 $ I-1, T, LDT, 300 $ T( 1, I ), 1 ) 301 T( I, I ) = TAU( I ) 302* 303 10 CONTINUE 304 A( K+NB, NB ) = EI 305* 306* Compute Y(1:K,1:NB) 307* 308 CALL ZLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY ) 309 CALL ZTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', 310 $ 'UNIT', K, NB, 311 $ ONE, A( K+1, 1 ), LDA, Y, LDY ) 312 IF( N.GT.K+NB ) 313 $ CALL ZGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, 314 $ NB, N-K-NB, ONE, 315 $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y, 316 $ LDY ) 317 CALL ZTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', 318 $ 'NON-UNIT', K, NB, 319 $ ONE, T, LDT, Y, LDY ) 320* 321 RETURN 322* 323* End of ZLAHR2 324* 325 END 326