![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLAHRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLAHRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLAHRD( 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 * COMPLEX A( LDA, * ), T( LDT, NB ), TAU( NB ), 00028 * $ Y( LDY, NB ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> CLAHRD reduces the first NB columns of a complex 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 a unitary similarity transformation 00040 *> Q**H * A * Q. The routine returns the matrices V and T which determine 00041 *> Q as a block reflector I - V*T*V**H, 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 CLAHR2 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 >= max(1,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 complexOTHERauxiliary 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**H 00143 *> 00144 *> where tau is a complex scalar, and v is a complex 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**H) * (A - Y*V**H). 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 CLAHRD( 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 COMPLEX A( LDA, * ), T( LDT, NB ), TAU( NB ), 00182 $ Y( LDY, NB ) 00183 * .. 00184 * 00185 * ===================================================================== 00186 * 00187 * .. Parameters .. 00188 COMPLEX ZERO, ONE 00189 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ), 00190 $ ONE = ( 1.0E+0, 0.0E+0 ) ) 00191 * .. 00192 * .. Local Scalars .. 00193 INTEGER I 00194 COMPLEX EI 00195 * .. 00196 * .. External Subroutines .. 00197 EXTERNAL CAXPY, CCOPY, CGEMV, CLACGV, CLARFG, CSCAL, 00198 $ CTRMV 00199 * .. 00200 * .. Intrinsic Functions .. 00201 INTRINSIC MIN 00202 * .. 00203 * .. Executable Statements .. 00204 * 00205 * Quick return if possible 00206 * 00207 IF( N.LE.1 ) 00208 $ RETURN 00209 * 00210 DO 10 I = 1, NB 00211 IF( I.GT.1 ) THEN 00212 * 00213 * Update A(1:n,i) 00214 * 00215 * Compute i-th column of A - Y * V**H 00216 * 00217 CALL CLACGV( I-1, A( K+I-1, 1 ), LDA ) 00218 CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, 00219 $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 ) 00220 CALL CLACGV( I-1, A( K+I-1, 1 ), LDA ) 00221 * 00222 * Apply I - V * T**H * V**H to this column (call it b) from the 00223 * left, using the last column of T as workspace 00224 * 00225 * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) 00226 * ( V2 ) ( b2 ) 00227 * 00228 * where V1 is unit lower triangular 00229 * 00230 * w := V1**H * b1 00231 * 00232 CALL CCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) 00233 CALL CTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1, 00234 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) 00235 * 00236 * w := w + V2**H *b2 00237 * 00238 CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE, 00239 $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE, 00240 $ T( 1, NB ), 1 ) 00241 * 00242 * w := T**H *w 00243 * 00244 CALL CTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1, 00245 $ T, LDT, T( 1, NB ), 1 ) 00246 * 00247 * b2 := b2 - V2*w 00248 * 00249 CALL CGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ), 00250 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) 00251 * 00252 * b1 := b1 - V1*w 00253 * 00254 CALL CTRMV( 'Lower', 'No transpose', 'Unit', I-1, 00255 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) 00256 CALL CAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) 00257 * 00258 A( K+I-1, I-1 ) = EI 00259 END IF 00260 * 00261 * Generate the elementary reflector H(i) to annihilate 00262 * A(k+i+1:n,i) 00263 * 00264 EI = A( K+I, I ) 00265 CALL CLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1, 00266 $ TAU( I ) ) 00267 A( K+I, I ) = ONE 00268 * 00269 * Compute Y(1:n,i) 00270 * 00271 CALL CGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA, 00272 $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 ) 00273 CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE, 00274 $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ), 00275 $ 1 ) 00276 CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1, 00277 $ ONE, Y( 1, I ), 1 ) 00278 CALL CSCAL( N, TAU( I ), Y( 1, I ), 1 ) 00279 * 00280 * Compute T(1:i,i) 00281 * 00282 CALL CSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) 00283 CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT, 00284 $ T( 1, I ), 1 ) 00285 T( I, I ) = TAU( I ) 00286 * 00287 10 CONTINUE 00288 A( K+NB, NB ) = EI 00289 * 00290 RETURN 00291 * 00292 * End of CLAHRD 00293 * 00294 END