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*> \date September 2012 127* 128*> \ingroup complex16OTHERauxiliary 129* 130*> \par Further Details: 131* ===================== 132*> 133*> \verbatim 134*> 135*> The matrix Q is represented as a product of nb elementary reflectors 136*> 137*> Q = H(1) H(2) . . . H(nb). 138*> 139*> Each H(i) has the form 140*> 141*> H(i) = I - tau * v * v**H 142*> 143*> where tau is a complex scalar, and v is a complex vector with 144*> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in 145*> A(i+k+1:n,i), and tau in TAU(i). 146*> 147*> The elements of the vectors v together form the (n-k+1)-by-nb matrix 148*> V which is needed, with T and Y, to apply the transformation to the 149*> unreduced part of the matrix, using an update of the form: 150*> A := (I - V*T*V**H) * (A - Y*V**H). 151*> 152*> The contents of A on exit are illustrated by the following example 153*> with n = 7, k = 3 and nb = 2: 154*> 155*> ( a a a a a ) 156*> ( a a a a a ) 157*> ( a a a a a ) 158*> ( h h a a a ) 159*> ( v1 h a a a ) 160*> ( v1 v2 a a a ) 161*> ( v1 v2 a a a ) 162*> 163*> where a denotes an element of the original matrix A, h denotes a 164*> modified element of the upper Hessenberg matrix H, and vi denotes an 165*> element of the vector defining H(i). 166*> 167*> This subroutine is a slight modification of LAPACK-3.0's DLAHRD 168*> incorporating improvements proposed by Quintana-Orti and Van de 169*> Gejin. Note that the entries of A(1:K,2:NB) differ from those 170*> returned by the original LAPACK-3.0's DLAHRD routine. (This 171*> subroutine is not backward compatible with LAPACK-3.0's DLAHRD.) 172*> \endverbatim 173* 174*> \par References: 175* ================ 176*> 177*> Gregorio Quintana-Orti and Robert van de Geijn, "Improving the 178*> performance of reduction to Hessenberg form," ACM Transactions on 179*> Mathematical Software, 32(2):180-194, June 2006. 180*> 181* ===================================================================== 182 SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) 183* 184* -- LAPACK auxiliary routine (version 3.4.2) -- 185* -- LAPACK is a software package provided by Univ. of Tennessee, -- 186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 187* September 2012 188* 189* .. Scalar Arguments .. 190 INTEGER K, LDA, LDT, LDY, N, NB 191* .. 192* .. Array Arguments .. 193 COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ), 194 $ Y( LDY, NB ) 195* .. 196* 197* ===================================================================== 198* 199* .. Parameters .. 200 COMPLEX*16 ZERO, ONE 201 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 202 $ ONE = ( 1.0D+0, 0.0D+0 ) ) 203* .. 204* .. Local Scalars .. 205 INTEGER I 206 COMPLEX*16 EI 207* .. 208* .. External Subroutines .. 209 EXTERNAL ZAXPY, ZCOPY, ZGEMM, ZGEMV, ZLACPY, 210 $ ZLARFG, ZSCAL, ZTRMM, ZTRMV, ZLACGV 211* .. 212* .. Intrinsic Functions .. 213 INTRINSIC MIN 214* .. 215* .. Executable Statements .. 216* 217* Quick return if possible 218* 219 IF( N.LE.1 ) 220 $ RETURN 221* 222 DO 10 I = 1, NB 223 IF( I.GT.1 ) THEN 224* 225* Update A(K+1:N,I) 226* 227* Update I-th column of A - Y * V**H 228* 229 CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) 230 CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY, 231 $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 ) 232 CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) 233* 234* Apply I - V * T**H * V**H to this column (call it b) from the 235* left, using the last column of T as workspace 236* 237* Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) 238* ( V2 ) ( b2 ) 239* 240* where V1 is unit lower triangular 241* 242* w := V1**H * b1 243* 244 CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) 245 CALL ZTRMV( 'Lower', 'Conjugate transpose', 'UNIT', 246 $ I-1, A( K+1, 1 ), 247 $ LDA, T( 1, NB ), 1 ) 248* 249* w := w + V2**H * b2 250* 251 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, 252 $ ONE, A( K+I, 1 ), 253 $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) 254* 255* w := T**H * w 256* 257 CALL ZTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT', 258 $ I-1, T, LDT, 259 $ T( 1, NB ), 1 ) 260* 261* b2 := b2 - V2*w 262* 263 CALL ZGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, 264 $ A( K+I, 1 ), 265 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) 266* 267* b1 := b1 - V1*w 268* 269 CALL ZTRMV( 'Lower', 'NO TRANSPOSE', 270 $ 'UNIT', I-1, 271 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) 272 CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) 273* 274 A( K+I-1, I-1 ) = EI 275 END IF 276* 277* Generate the elementary reflector H(I) to annihilate 278* A(K+I+1:N,I) 279* 280 CALL ZLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, 281 $ TAU( I ) ) 282 EI = A( K+I, I ) 283 A( K+I, I ) = ONE 284* 285* Compute Y(K+1:N,I) 286* 287 CALL ZGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, 288 $ ONE, A( K+1, I+1 ), 289 $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 ) 290 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, 291 $ ONE, A( K+I, 1 ), LDA, 292 $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) 293 CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, 294 $ Y( K+1, 1 ), LDY, 295 $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 ) 296 CALL ZSCAL( N-K, TAU( I ), Y( K+1, I ), 1 ) 297* 298* Compute T(1:I,I) 299* 300 CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) 301 CALL ZTRMV( 'Upper', 'No Transpose', 'NON-UNIT', 302 $ I-1, T, LDT, 303 $ T( 1, I ), 1 ) 304 T( I, I ) = TAU( I ) 305* 306 10 CONTINUE 307 A( K+NB, NB ) = EI 308* 309* Compute Y(1:K,1:NB) 310* 311 CALL ZLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY ) 312 CALL ZTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', 313 $ 'UNIT', K, NB, 314 $ ONE, A( K+1, 1 ), LDA, Y, LDY ) 315 IF( N.GT.K+NB ) 316 $ CALL ZGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, 317 $ NB, N-K-NB, ONE, 318 $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y, 319 $ LDY ) 320 CALL ZTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', 321 $ 'NON-UNIT', K, NB, 322 $ ONE, T, LDT, Y, LDY ) 323* 324 RETURN 325* 326* End of ZLAHR2 327* 328 END 329