![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAHRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLAHRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER K, LDA, LDT, LDY, N, NB 00025 * .. 00026 * .. Array Arguments .. 00027 * DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), 00028 * $ Y( LDY, NB ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> DLAHRD reduces the first NB columns of a real general n-by-(n-k+1) 00038 *> matrix A so that elements below the k-th subdiagonal are zero. The 00039 *> reduction is performed by an orthogonal similarity transformation 00040 *> Q**T * A * Q. The routine returns the matrices V and T which determine 00041 *> Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T. 00042 *> 00043 *> This is an OBSOLETE auxiliary routine. 00044 *> This routine will be 'deprecated' in a future release. 00045 *> Please use the new routine DLAHR2 instead. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] N 00052 *> \verbatim 00053 *> N is INTEGER 00054 *> The order of the matrix A. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] K 00058 *> \verbatim 00059 *> K is INTEGER 00060 *> The offset for the reduction. Elements below the k-th 00061 *> subdiagonal in the first NB columns are reduced to zero. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] NB 00065 *> \verbatim 00066 *> NB is INTEGER 00067 *> The number of columns to be reduced. 00068 *> \endverbatim 00069 *> 00070 *> \param[in,out] A 00071 *> \verbatim 00072 *> A is DOUBLE PRECISION array, dimension (LDA,N-K+1) 00073 *> On entry, the n-by-(n-k+1) general matrix A. 00074 *> On exit, the elements on and above the k-th subdiagonal in 00075 *> the first NB columns are overwritten with the corresponding 00076 *> elements of the reduced matrix; the elements below the k-th 00077 *> subdiagonal, with the array TAU, represent the matrix Q as a 00078 *> product of elementary reflectors. The other columns of A are 00079 *> unchanged. See Further Details. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] LDA 00083 *> \verbatim 00084 *> LDA is INTEGER 00085 *> The leading dimension of the array A. LDA >= max(1,N). 00086 *> \endverbatim 00087 *> 00088 *> \param[out] TAU 00089 *> \verbatim 00090 *> TAU is DOUBLE PRECISION array, dimension (NB) 00091 *> The scalar factors of the elementary reflectors. See Further 00092 *> Details. 00093 *> \endverbatim 00094 *> 00095 *> \param[out] T 00096 *> \verbatim 00097 *> T is DOUBLE PRECISION array, dimension (LDT,NB) 00098 *> The upper triangular matrix T. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDT 00102 *> \verbatim 00103 *> LDT is INTEGER 00104 *> The leading dimension of the array T. LDT >= NB. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] Y 00108 *> \verbatim 00109 *> Y is DOUBLE PRECISION array, dimension (LDY,NB) 00110 *> The n-by-nb matrix Y. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] LDY 00114 *> \verbatim 00115 *> LDY is INTEGER 00116 *> The leading dimension of the array Y. LDY >= N. 00117 *> \endverbatim 00118 * 00119 * Authors: 00120 * ======== 00121 * 00122 *> \author Univ. of Tennessee 00123 *> \author Univ. of California Berkeley 00124 *> \author Univ. of Colorado Denver 00125 *> \author NAG Ltd. 00126 * 00127 *> \date November 2011 00128 * 00129 *> \ingroup doubleOTHERauxiliary 00130 * 00131 *> \par Further Details: 00132 * ===================== 00133 *> 00134 *> \verbatim 00135 *> 00136 *> The matrix Q is represented as a product of nb elementary reflectors 00137 *> 00138 *> Q = H(1) H(2) . . . H(nb). 00139 *> 00140 *> Each H(i) has the form 00141 *> 00142 *> H(i) = I - tau * v * v**T 00143 *> 00144 *> where tau is a real scalar, and v is a real vector with 00145 *> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in 00146 *> A(i+k+1:n,i), and tau in TAU(i). 00147 *> 00148 *> The elements of the vectors v together form the (n-k+1)-by-nb matrix 00149 *> V which is needed, with T and Y, to apply the transformation to the 00150 *> unreduced part of the matrix, using an update of the form: 00151 *> A := (I - V*T*V**T) * (A - Y*V**T). 00152 *> 00153 *> The contents of A on exit are illustrated by the following example 00154 *> with n = 7, k = 3 and nb = 2: 00155 *> 00156 *> ( a h a a a ) 00157 *> ( a h a a a ) 00158 *> ( a h a a a ) 00159 *> ( h h a a a ) 00160 *> ( v1 h a a a ) 00161 *> ( v1 v2 a a a ) 00162 *> ( v1 v2 a a a ) 00163 *> 00164 *> where a denotes an element of the original matrix A, h denotes a 00165 *> modified element of the upper Hessenberg matrix H, and vi denotes an 00166 *> element of the vector defining H(i). 00167 *> \endverbatim 00168 *> 00169 * ===================================================================== 00170 SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) 00171 * 00172 * -- LAPACK auxiliary routine (version 3.4.0) -- 00173 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00174 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00175 * November 2011 00176 * 00177 * .. Scalar Arguments .. 00178 INTEGER K, LDA, LDT, LDY, N, NB 00179 * .. 00180 * .. Array Arguments .. 00181 DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), 00182 $ Y( LDY, NB ) 00183 * .. 00184 * 00185 * ===================================================================== 00186 * 00187 * .. Parameters .. 00188 DOUBLE PRECISION ZERO, ONE 00189 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00190 * .. 00191 * .. Local Scalars .. 00192 INTEGER I 00193 DOUBLE PRECISION EI 00194 * .. 00195 * .. External Subroutines .. 00196 EXTERNAL DAXPY, DCOPY, DGEMV, DLARFG, DSCAL, DTRMV 00197 * .. 00198 * .. Intrinsic Functions .. 00199 INTRINSIC MIN 00200 * .. 00201 * .. Executable Statements .. 00202 * 00203 * Quick return if possible 00204 * 00205 IF( N.LE.1 ) 00206 $ RETURN 00207 * 00208 DO 10 I = 1, NB 00209 IF( I.GT.1 ) THEN 00210 * 00211 * Update A(1:n,i) 00212 * 00213 * Compute i-th column of A - Y * V**T 00214 * 00215 CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, 00216 $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 ) 00217 * 00218 * Apply I - V * T**T * V**T to this column (call it b) from the 00219 * left, using the last column of T as workspace 00220 * 00221 * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) 00222 * ( V2 ) ( b2 ) 00223 * 00224 * where V1 is unit lower triangular 00225 * 00226 * w := V1**T * b1 00227 * 00228 CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) 00229 CALL DTRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ), 00230 $ LDA, T( 1, NB ), 1 ) 00231 * 00232 * w := w + V2**T *b2 00233 * 00234 CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), 00235 $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) 00236 * 00237 * w := T**T *w 00238 * 00239 CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT, 00240 $ T( 1, NB ), 1 ) 00241 * 00242 * b2 := b2 - V2*w 00243 * 00244 CALL DGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ), 00245 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) 00246 * 00247 * b1 := b1 - V1*w 00248 * 00249 CALL DTRMV( 'Lower', 'No transpose', 'Unit', I-1, 00250 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) 00251 CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) 00252 * 00253 A( K+I-1, I-1 ) = EI 00254 END IF 00255 * 00256 * Generate the elementary reflector H(i) to annihilate 00257 * A(k+i+1:n,i) 00258 * 00259 CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, 00260 $ TAU( I ) ) 00261 EI = A( K+I, I ) 00262 A( K+I, I ) = ONE 00263 * 00264 * Compute Y(1:n,i) 00265 * 00266 CALL DGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA, 00267 $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 ) 00268 CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA, 00269 $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) 00270 CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1, 00271 $ ONE, Y( 1, I ), 1 ) 00272 CALL DSCAL( N, TAU( I ), Y( 1, I ), 1 ) 00273 * 00274 * Compute T(1:i,i) 00275 * 00276 CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) 00277 CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT, 00278 $ T( 1, I ), 1 ) 00279 T( I, I ) = TAU( I ) 00280 * 00281 10 CONTINUE 00282 A( K+NB, NB ) = EI 00283 * 00284 RETURN 00285 * 00286 * End of DLAHRD 00287 * 00288 END