![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGEHRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGEHRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgehrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgehrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgehrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER IHI, ILO, INFO, LDA, LWORK, N 00025 * .. 00026 * .. Array Arguments .. 00027 * COMPLEX A( LDA, * ), TAU( * ), WORK( * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> CGEHRD reduces a complex general matrix A to upper Hessenberg form H by 00037 *> an unitary similarity transformation: Q**H * A * Q = H . 00038 *> \endverbatim 00039 * 00040 * Arguments: 00041 * ========== 00042 * 00043 *> \param[in] N 00044 *> \verbatim 00045 *> N is INTEGER 00046 *> The order of the matrix A. N >= 0. 00047 *> \endverbatim 00048 *> 00049 *> \param[in] ILO 00050 *> \verbatim 00051 *> ILO is INTEGER 00052 *> \endverbatim 00053 *> 00054 *> \param[in] IHI 00055 *> \verbatim 00056 *> IHI is INTEGER 00057 *> 00058 *> It is assumed that A is already upper triangular in rows 00059 *> and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally 00060 *> set by a previous call to CGEBAL; otherwise they should be 00061 *> set to 1 and N respectively. See Further Details. 00062 *> 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 00063 *> \endverbatim 00064 *> 00065 *> \param[in,out] A 00066 *> \verbatim 00067 *> A is COMPLEX array, dimension (LDA,N) 00068 *> On entry, the N-by-N general matrix to be reduced. 00069 *> On exit, the upper triangle and the first subdiagonal of A 00070 *> are overwritten with the upper Hessenberg matrix H, and the 00071 *> elements below the first subdiagonal, with the array TAU, 00072 *> represent the unitary matrix Q as a product of elementary 00073 *> reflectors. See Further Details. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] LDA 00077 *> \verbatim 00078 *> LDA is INTEGER 00079 *> The leading dimension of the array A. LDA >= max(1,N). 00080 *> \endverbatim 00081 *> 00082 *> \param[out] TAU 00083 *> \verbatim 00084 *> TAU is COMPLEX array, dimension (N-1) 00085 *> The scalar factors of the elementary reflectors (see Further 00086 *> Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to 00087 *> zero. 00088 *> \endverbatim 00089 *> 00090 *> \param[out] WORK 00091 *> \verbatim 00092 *> WORK is COMPLEX array, dimension (LWORK) 00093 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] LWORK 00097 *> \verbatim 00098 *> LWORK is INTEGER 00099 *> The length of the array WORK. LWORK >= max(1,N). 00100 *> For optimum performance LWORK >= N*NB, where NB is the 00101 *> optimal blocksize. 00102 *> 00103 *> If LWORK = -1, then a workspace query is assumed; the routine 00104 *> only calculates the optimal size of the WORK array, returns 00105 *> this value as the first entry of the WORK array, and no error 00106 *> message related to LWORK is issued by XERBLA. 00107 *> \endverbatim 00108 *> 00109 *> \param[out] INFO 00110 *> \verbatim 00111 *> INFO is INTEGER 00112 *> = 0: successful exit 00113 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00114 *> \endverbatim 00115 * 00116 * Authors: 00117 * ======== 00118 * 00119 *> \author Univ. of Tennessee 00120 *> \author Univ. of California Berkeley 00121 *> \author Univ. of Colorado Denver 00122 *> \author NAG Ltd. 00123 * 00124 *> \date November 2011 00125 * 00126 *> \ingroup complexGEcomputational 00127 * 00128 *> \par Further Details: 00129 * ===================== 00130 *> 00131 *> \verbatim 00132 *> 00133 *> The matrix Q is represented as a product of (ihi-ilo) elementary 00134 *> reflectors 00135 *> 00136 *> Q = H(ilo) H(ilo+1) . . . H(ihi-1). 00137 *> 00138 *> Each H(i) has the form 00139 *> 00140 *> H(i) = I - tau * v * v**H 00141 *> 00142 *> where tau is a complex scalar, and v is a complex vector with 00143 *> v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on 00144 *> exit in A(i+2:ihi,i), and tau in TAU(i). 00145 *> 00146 *> The contents of A are illustrated by the following example, with 00147 *> n = 7, ilo = 2 and ihi = 6: 00148 *> 00149 *> on entry, on exit, 00150 *> 00151 *> ( a a a a a a a ) ( a a h h h h a ) 00152 *> ( a a a a a a ) ( a h h h h a ) 00153 *> ( a a a a a a ) ( h h h h h h ) 00154 *> ( a a a a a a ) ( v2 h h h h h ) 00155 *> ( a a a a a a ) ( v2 v3 h h h h ) 00156 *> ( a a a a a a ) ( v2 v3 v4 h h h ) 00157 *> ( a ) ( a ) 00158 *> 00159 *> where a denotes an element of the original matrix A, h denotes a 00160 *> modified element of the upper Hessenberg matrix H, and vi denotes an 00161 *> element of the vector defining H(i). 00162 *> 00163 *> This file is a slight modification of LAPACK-3.0's DGEHRD 00164 *> subroutine incorporating improvements proposed by Quintana-Orti and 00165 *> Van de Geijn (2006). (See DLAHR2.) 00166 *> \endverbatim 00167 *> 00168 * ===================================================================== 00169 SUBROUTINE CGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO ) 00170 * 00171 * -- LAPACK computational routine (version 3.4.0) -- 00172 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00174 * November 2011 00175 * 00176 * .. Scalar Arguments .. 00177 INTEGER IHI, ILO, INFO, LDA, LWORK, N 00178 * .. 00179 * .. Array Arguments .. 00180 COMPLEX A( LDA, * ), TAU( * ), WORK( * ) 00181 * .. 00182 * 00183 * ===================================================================== 00184 * 00185 * .. Parameters .. 00186 INTEGER NBMAX, LDT 00187 PARAMETER ( NBMAX = 64, LDT = NBMAX+1 ) 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 LOGICAL LQUERY 00194 INTEGER I, IB, IINFO, IWS, J, LDWORK, LWKOPT, NB, 00195 $ NBMIN, NH, NX 00196 COMPLEX EI 00197 * .. 00198 * .. Local Arrays .. 00199 COMPLEX T( LDT, NBMAX ) 00200 * .. 00201 * .. External Subroutines .. 00202 EXTERNAL CAXPY, CGEHD2, CGEMM, CLAHR2, CLARFB, CTRMM, 00203 $ XERBLA 00204 * .. 00205 * .. Intrinsic Functions .. 00206 INTRINSIC MAX, MIN 00207 * .. 00208 * .. External Functions .. 00209 INTEGER ILAENV 00210 EXTERNAL ILAENV 00211 * .. 00212 * .. Executable Statements .. 00213 * 00214 * Test the input parameters 00215 * 00216 INFO = 0 00217 NB = MIN( NBMAX, ILAENV( 1, 'CGEHRD', ' ', N, ILO, IHI, -1 ) ) 00218 LWKOPT = N*NB 00219 WORK( 1 ) = LWKOPT 00220 LQUERY = ( LWORK.EQ.-1 ) 00221 IF( N.LT.0 ) THEN 00222 INFO = -1 00223 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN 00224 INFO = -2 00225 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN 00226 INFO = -3 00227 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00228 INFO = -5 00229 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00230 INFO = -8 00231 END IF 00232 IF( INFO.NE.0 ) THEN 00233 CALL XERBLA( 'CGEHRD', -INFO ) 00234 RETURN 00235 ELSE IF( LQUERY ) THEN 00236 RETURN 00237 END IF 00238 * 00239 * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero 00240 * 00241 DO 10 I = 1, ILO - 1 00242 TAU( I ) = ZERO 00243 10 CONTINUE 00244 DO 20 I = MAX( 1, IHI ), N - 1 00245 TAU( I ) = ZERO 00246 20 CONTINUE 00247 * 00248 * Quick return if possible 00249 * 00250 NH = IHI - ILO + 1 00251 IF( NH.LE.1 ) THEN 00252 WORK( 1 ) = 1 00253 RETURN 00254 END IF 00255 * 00256 * Determine the block size 00257 * 00258 NB = MIN( NBMAX, ILAENV( 1, 'CGEHRD', ' ', N, ILO, IHI, -1 ) ) 00259 NBMIN = 2 00260 IWS = 1 00261 IF( NB.GT.1 .AND. NB.LT.NH ) THEN 00262 * 00263 * Determine when to cross over from blocked to unblocked code 00264 * (last block is always handled by unblocked code) 00265 * 00266 NX = MAX( NB, ILAENV( 3, 'CGEHRD', ' ', N, ILO, IHI, -1 ) ) 00267 IF( NX.LT.NH ) THEN 00268 * 00269 * Determine if workspace is large enough for blocked code 00270 * 00271 IWS = N*NB 00272 IF( LWORK.LT.IWS ) THEN 00273 * 00274 * Not enough workspace to use optimal NB: determine the 00275 * minimum value of NB, and reduce NB or force use of 00276 * unblocked code 00277 * 00278 NBMIN = MAX( 2, ILAENV( 2, 'CGEHRD', ' ', N, ILO, IHI, 00279 $ -1 ) ) 00280 IF( LWORK.GE.N*NBMIN ) THEN 00281 NB = LWORK / N 00282 ELSE 00283 NB = 1 00284 END IF 00285 END IF 00286 END IF 00287 END IF 00288 LDWORK = N 00289 * 00290 IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN 00291 * 00292 * Use unblocked code below 00293 * 00294 I = ILO 00295 * 00296 ELSE 00297 * 00298 * Use blocked code 00299 * 00300 DO 40 I = ILO, IHI - 1 - NX, NB 00301 IB = MIN( NB, IHI-I ) 00302 * 00303 * Reduce columns i:i+ib-1 to Hessenberg form, returning the 00304 * matrices V and T of the block reflector H = I - V*T*V**H 00305 * which performs the reduction, and also the matrix Y = A*V*T 00306 * 00307 CALL CLAHR2( IHI, I, IB, A( 1, I ), LDA, TAU( I ), T, LDT, 00308 $ WORK, LDWORK ) 00309 * 00310 * Apply the block reflector H to A(1:ihi,i+ib:ihi) from the 00311 * right, computing A := A - Y * V**H. V(i+ib,ib-1) must be set 00312 * to 1 00313 * 00314 EI = A( I+IB, I+IB-1 ) 00315 A( I+IB, I+IB-1 ) = ONE 00316 CALL CGEMM( 'No transpose', 'Conjugate transpose', 00317 $ IHI, IHI-I-IB+1, 00318 $ IB, -ONE, WORK, LDWORK, A( I+IB, I ), LDA, ONE, 00319 $ A( 1, I+IB ), LDA ) 00320 A( I+IB, I+IB-1 ) = EI 00321 * 00322 * Apply the block reflector H to A(1:i,i+1:i+ib-1) from the 00323 * right 00324 * 00325 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 00326 $ 'Unit', I, IB-1, 00327 $ ONE, A( I+1, I ), LDA, WORK, LDWORK ) 00328 DO 30 J = 0, IB-2 00329 CALL CAXPY( I, -ONE, WORK( LDWORK*J+1 ), 1, 00330 $ A( 1, I+J+1 ), 1 ) 00331 30 CONTINUE 00332 * 00333 * Apply the block reflector H to A(i+1:ihi,i+ib:n) from the 00334 * left 00335 * 00336 CALL CLARFB( 'Left', 'Conjugate transpose', 'Forward', 00337 $ 'Columnwise', 00338 $ IHI-I, N-I-IB+1, IB, A( I+1, I ), LDA, T, LDT, 00339 $ A( I+1, I+IB ), LDA, WORK, LDWORK ) 00340 40 CONTINUE 00341 END IF 00342 * 00343 * Use unblocked code to reduce the rest of the matrix 00344 * 00345 CALL CGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO ) 00346 WORK( 1 ) = IWS 00347 * 00348 RETURN 00349 * 00350 * End of CGEHRD 00351 * 00352 END