1*> \brief \b DLAHR2 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 DLAHR2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahr2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahr2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahr2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAHR2( 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* DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), 28* $ Y( LDY, NB ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> DLAHR2 reduces the first NB columns of A real 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 orthogonal similarity transformation 40*> Q**T * A * Q. The routine returns the matrices V and T which determine 41*> Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T. 42*> 43*> This is an auxiliary routine called by DGEHRD. 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERauxiliary 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**T 142*> 143*> where tau is a real scalar, and v is a real 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**T) * (A - Y*V**T). 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 DLAHR2( 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 DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), 194 $ Y( LDY, NB ) 195* .. 196* 197* ===================================================================== 198* 199* .. Parameters .. 200 DOUBLE PRECISION ZERO, ONE 201 PARAMETER ( ZERO = 0.0D+0, 202 $ ONE = 1.0D+0 ) 203* .. 204* .. Local Scalars .. 205 INTEGER I 206 DOUBLE PRECISION EI 207* .. 208* .. External Subroutines .. 209 EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DLACPY, 210 $ DLARFG, DSCAL, DTRMM, DTRMV 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**T 228* 229 CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY, 230 $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 ) 231* 232* Apply I - V * T**T * V**T to this column (call it b) from the 233* left, using the last column of T as workspace 234* 235* Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) 236* ( V2 ) ( b2 ) 237* 238* where V1 is unit lower triangular 239* 240* w := V1**T * b1 241* 242 CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) 243 CALL DTRMV( 'Lower', 'Transpose', 'UNIT', 244 $ I-1, A( K+1, 1 ), 245 $ LDA, T( 1, NB ), 1 ) 246* 247* w := w + V2**T * b2 248* 249 CALL DGEMV( 'Transpose', N-K-I+1, I-1, 250 $ ONE, A( K+I, 1 ), 251 $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) 252* 253* w := T**T * w 254* 255 CALL DTRMV( 'Upper', 'Transpose', 'NON-UNIT', 256 $ I-1, T, LDT, 257 $ T( 1, NB ), 1 ) 258* 259* b2 := b2 - V2*w 260* 261 CALL DGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, 262 $ A( K+I, 1 ), 263 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) 264* 265* b1 := b1 - V1*w 266* 267 CALL DTRMV( 'Lower', 'NO TRANSPOSE', 268 $ 'UNIT', I-1, 269 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) 270 CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) 271* 272 A( K+I-1, I-1 ) = EI 273 END IF 274* 275* Generate the elementary reflector H(I) to annihilate 276* A(K+I+1:N,I) 277* 278 CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, 279 $ TAU( I ) ) 280 EI = A( K+I, I ) 281 A( K+I, I ) = ONE 282* 283* Compute Y(K+1:N,I) 284* 285 CALL DGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, 286 $ ONE, A( K+1, I+1 ), 287 $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 ) 288 CALL DGEMV( 'Transpose', N-K-I+1, I-1, 289 $ ONE, A( K+I, 1 ), LDA, 290 $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) 291 CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, 292 $ Y( K+1, 1 ), LDY, 293 $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 ) 294 CALL DSCAL( N-K, TAU( I ), Y( K+1, I ), 1 ) 295* 296* Compute T(1:I,I) 297* 298 CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) 299 CALL DTRMV( 'Upper', 'No Transpose', 'NON-UNIT', 300 $ I-1, T, LDT, 301 $ T( 1, I ), 1 ) 302 T( I, I ) = TAU( I ) 303* 304 10 CONTINUE 305 A( K+NB, NB ) = EI 306* 307* Compute Y(1:K,1:NB) 308* 309 CALL DLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY ) 310 CALL DTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', 311 $ 'UNIT', K, NB, 312 $ ONE, A( K+1, 1 ), LDA, Y, LDY ) 313 IF( N.GT.K+NB ) 314 $ CALL DGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, 315 $ NB, N-K-NB, ONE, 316 $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y, 317 $ LDY ) 318 CALL DTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', 319 $ 'NON-UNIT', K, NB, 320 $ ONE, T, LDT, Y, LDY ) 321* 322 RETURN 323* 324* End of DLAHR2 325* 326 END 327