![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHSEQR 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHSEQR + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhseqr.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhseqr.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhseqr.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, 00022 * WORK, LWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N 00026 * CHARACTER COMPZ, JOB 00027 * .. 00028 * .. Array Arguments .. 00029 * COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> ZHSEQR computes the eigenvalues of a Hessenberg matrix H 00039 *> and, optionally, the matrices T and Z from the Schur decomposition 00040 *> H = Z T Z**H, where T is an upper triangular matrix (the 00041 *> Schur form), and Z is the unitary matrix of Schur vectors. 00042 *> 00043 *> Optionally Z may be postmultiplied into an input unitary 00044 *> matrix Q so that this routine can give the Schur factorization 00045 *> of a matrix A which has been reduced to the Hessenberg form H 00046 *> by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] JOB 00053 *> \verbatim 00054 *> JOB is CHARACTER*1 00055 *> = 'E': compute eigenvalues only; 00056 *> = 'S': compute eigenvalues and the Schur form T. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] COMPZ 00060 *> \verbatim 00061 *> COMPZ is CHARACTER*1 00062 *> = 'N': no Schur vectors are computed; 00063 *> = 'I': Z is initialized to the unit matrix and the matrix Z 00064 *> of Schur vectors of H is returned; 00065 *> = 'V': Z must contain an unitary matrix Q on entry, and 00066 *> the product Q*Z is returned. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] N 00070 *> \verbatim 00071 *> N is INTEGER 00072 *> The order of the matrix H. N .GE. 0. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] ILO 00076 *> \verbatim 00077 *> ILO is INTEGER 00078 *> \endverbatim 00079 *> 00080 *> \param[in] IHI 00081 *> \verbatim 00082 *> IHI is INTEGER 00083 *> 00084 *> It is assumed that H is already upper triangular in rows 00085 *> and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally 00086 *> set by a previous call to ZGEBAL, and then passed to ZGEHRD 00087 *> when the matrix output by ZGEBAL is reduced to Hessenberg 00088 *> form. Otherwise ILO and IHI should be set to 1 and N 00089 *> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. 00090 *> If N = 0, then ILO = 1 and IHI = 0. 00091 *> \endverbatim 00092 *> 00093 *> \param[in,out] H 00094 *> \verbatim 00095 *> H is COMPLEX*16 array, dimension (LDH,N) 00096 *> On entry, the upper Hessenberg matrix H. 00097 *> On exit, if INFO = 0 and JOB = 'S', H contains the upper 00098 *> triangular matrix T from the Schur decomposition (the 00099 *> Schur form). If INFO = 0 and JOB = 'E', the contents of 00100 *> H are unspecified on exit. (The output value of H when 00101 *> INFO.GT.0 is given under the description of INFO below.) 00102 *> 00103 *> Unlike earlier versions of ZHSEQR, this subroutine may 00104 *> explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1 00105 *> or j = IHI+1, IHI+2, ... N. 00106 *> \endverbatim 00107 *> 00108 *> \param[in] LDH 00109 *> \verbatim 00110 *> LDH is INTEGER 00111 *> The leading dimension of the array H. LDH .GE. max(1,N). 00112 *> \endverbatim 00113 *> 00114 *> \param[out] W 00115 *> \verbatim 00116 *> W is COMPLEX*16 array, dimension (N) 00117 *> The computed eigenvalues. If JOB = 'S', the eigenvalues are 00118 *> stored in the same order as on the diagonal of the Schur 00119 *> form returned in H, with W(i) = H(i,i). 00120 *> \endverbatim 00121 *> 00122 *> \param[in,out] Z 00123 *> \verbatim 00124 *> Z is COMPLEX*16 array, dimension (LDZ,N) 00125 *> If COMPZ = 'N', Z is not referenced. 00126 *> If COMPZ = 'I', on entry Z need not be set and on exit, 00127 *> if INFO = 0, Z contains the unitary matrix Z of the Schur 00128 *> vectors of H. If COMPZ = 'V', on entry Z must contain an 00129 *> N-by-N matrix Q, which is assumed to be equal to the unit 00130 *> matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit, 00131 *> if INFO = 0, Z contains Q*Z. 00132 *> Normally Q is the unitary matrix generated by ZUNGHR 00133 *> after the call to ZGEHRD which formed the Hessenberg matrix 00134 *> H. (The output value of Z when INFO.GT.0 is given under 00135 *> the description of INFO below.) 00136 *> \endverbatim 00137 *> 00138 *> \param[in] LDZ 00139 *> \verbatim 00140 *> LDZ is INTEGER 00141 *> The leading dimension of the array Z. if COMPZ = 'I' or 00142 *> COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1. 00143 *> \endverbatim 00144 *> 00145 *> \param[out] WORK 00146 *> \verbatim 00147 *> WORK is COMPLEX*16 array, dimension (LWORK) 00148 *> On exit, if INFO = 0, WORK(1) returns an estimate of 00149 *> the optimal value for LWORK. 00150 *> \endverbatim 00151 *> 00152 *> \param[in] LWORK 00153 *> \verbatim 00154 *> LWORK is INTEGER 00155 *> The dimension of the array WORK. LWORK .GE. max(1,N) 00156 *> is sufficient and delivers very good and sometimes 00157 *> optimal performance. However, LWORK as large as 11*N 00158 *> may be required for optimal performance. A workspace 00159 *> query is recommended to determine the optimal workspace 00160 *> size. 00161 *> 00162 *> If LWORK = -1, then ZHSEQR does a workspace query. 00163 *> In this case, ZHSEQR checks the input parameters and 00164 *> estimates the optimal workspace size for the given 00165 *> values of N, ILO and IHI. The estimate is returned 00166 *> in WORK(1). No error message related to LWORK is 00167 *> issued by XERBLA. Neither H nor Z are accessed. 00168 *> \endverbatim 00169 *> 00170 *> \param[out] INFO 00171 *> \verbatim 00172 *> INFO is INTEGER 00173 *> = 0: successful exit 00174 *> .LT. 0: if INFO = -i, the i-th argument had an illegal 00175 *> value 00176 *> .GT. 0: if INFO = i, ZHSEQR failed to compute all of 00177 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR 00178 *> and WI contain those eigenvalues which have been 00179 *> successfully computed. (Failures are rare.) 00180 *> 00181 *> If INFO .GT. 0 and JOB = 'E', then on exit, the 00182 *> remaining unconverged eigenvalues are the eigen- 00183 *> values of the upper Hessenberg matrix rows and 00184 *> columns ILO through INFO of the final, output 00185 *> value of H. 00186 *> 00187 *> If INFO .GT. 0 and JOB = 'S', then on exit 00188 *> 00189 *> (*) (initial value of H)*U = U*(final value of H) 00190 *> 00191 *> where U is a unitary matrix. The final 00192 *> value of H is upper Hessenberg and triangular in 00193 *> rows and columns INFO+1 through IHI. 00194 *> 00195 *> If INFO .GT. 0 and COMPZ = 'V', then on exit 00196 *> 00197 *> (final value of Z) = (initial value of Z)*U 00198 *> 00199 *> where U is the unitary matrix in (*) (regard- 00200 *> less of the value of JOB.) 00201 *> 00202 *> If INFO .GT. 0 and COMPZ = 'I', then on exit 00203 *> (final value of Z) = U 00204 *> where U is the unitary matrix in (*) (regard- 00205 *> less of the value of JOB.) 00206 *> 00207 *> If INFO .GT. 0 and COMPZ = 'N', then Z is not 00208 *> accessed. 00209 *> \endverbatim 00210 * 00211 * Authors: 00212 * ======== 00213 * 00214 *> \author Univ. of Tennessee 00215 *> \author Univ. of California Berkeley 00216 *> \author Univ. of Colorado Denver 00217 *> \author NAG Ltd. 00218 * 00219 *> \date November 2011 00220 * 00221 *> \ingroup complex16OTHERcomputational 00222 * 00223 *> \par Contributors: 00224 * ================== 00225 *> 00226 *> Karen Braman and Ralph Byers, Department of Mathematics, 00227 *> University of Kansas, USA 00228 * 00229 *> \par Further Details: 00230 * ===================== 00231 *> 00232 *> \verbatim 00233 *> 00234 *> Default values supplied by 00235 *> ILAENV(ISPEC,'ZHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK). 00236 *> It is suggested that these defaults be adjusted in order 00237 *> to attain best performance in each particular 00238 *> computational environment. 00239 *> 00240 *> ISPEC=12: The ZLAHQR vs ZLAQR0 crossover point. 00241 *> Default: 75. (Must be at least 11.) 00242 *> 00243 *> ISPEC=13: Recommended deflation window size. 00244 *> This depends on ILO, IHI and NS. NS is the 00245 *> number of simultaneous shifts returned 00246 *> by ILAENV(ISPEC=15). (See ISPEC=15 below.) 00247 *> The default for (IHI-ILO+1).LE.500 is NS. 00248 *> The default for (IHI-ILO+1).GT.500 is 3*NS/2. 00249 *> 00250 *> ISPEC=14: Nibble crossover point. (See IPARMQ for 00251 *> details.) Default: 14% of deflation window 00252 *> size. 00253 *> 00254 *> ISPEC=15: Number of simultaneous shifts in a multishift 00255 *> QR iteration. 00256 *> 00257 *> If IHI-ILO+1 is ... 00258 *> 00259 *> greater than ...but less ... the 00260 *> or equal to ... than default is 00261 *> 00262 *> 1 30 NS = 2(+) 00263 *> 30 60 NS = 4(+) 00264 *> 60 150 NS = 10(+) 00265 *> 150 590 NS = ** 00266 *> 590 3000 NS = 64 00267 *> 3000 6000 NS = 128 00268 *> 6000 infinity NS = 256 00269 *> 00270 *> (+) By default some or all matrices of this order 00271 *> are passed to the implicit double shift routine 00272 *> ZLAHQR and this parameter is ignored. See 00273 *> ISPEC=12 above and comments in IPARMQ for 00274 *> details. 00275 *> 00276 *> (**) The asterisks (**) indicate an ad-hoc 00277 *> function of N increasing from 10 to 64. 00278 *> 00279 *> ISPEC=16: Select structured matrix multiply. 00280 *> If the number of simultaneous shifts (specified 00281 *> by ISPEC=15) is less than 14, then the default 00282 *> for ISPEC=16 is 0. Otherwise the default for 00283 *> ISPEC=16 is 2. 00284 *> \endverbatim 00285 * 00286 *> \par References: 00287 * ================ 00288 *> 00289 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 00290 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 00291 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages 00292 *> 929--947, 2002. 00293 *> \n 00294 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 00295 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal 00296 *> of Matrix Analysis, volume 23, pages 948--973, 2002. 00297 * 00298 * ===================================================================== 00299 SUBROUTINE ZHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, 00300 $ WORK, LWORK, INFO ) 00301 * 00302 * -- LAPACK computational routine (version 3.4.0) -- 00303 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00304 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00305 * November 2011 00306 * 00307 * .. Scalar Arguments .. 00308 INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N 00309 CHARACTER COMPZ, JOB 00310 * .. 00311 * .. Array Arguments .. 00312 COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) 00313 * .. 00314 * 00315 * ===================================================================== 00316 * 00317 * .. Parameters .. 00318 * 00319 * ==== Matrices of order NTINY or smaller must be processed by 00320 * . ZLAHQR because of insufficient subdiagonal scratch space. 00321 * . (This is a hard limit.) ==== 00322 INTEGER NTINY 00323 PARAMETER ( NTINY = 11 ) 00324 * 00325 * ==== NL allocates some local workspace to help small matrices 00326 * . through a rare ZLAHQR failure. NL .GT. NTINY = 11 is 00327 * . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom- 00328 * . mended. (The default value of NMIN is 75.) Using NL = 49 00329 * . allows up to six simultaneous shifts and a 16-by-16 00330 * . deflation window. ==== 00331 INTEGER NL 00332 PARAMETER ( NL = 49 ) 00333 COMPLEX*16 ZERO, ONE 00334 PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ), 00335 $ ONE = ( 1.0d0, 0.0d0 ) ) 00336 DOUBLE PRECISION RZERO 00337 PARAMETER ( RZERO = 0.0d0 ) 00338 * .. 00339 * .. Local Arrays .. 00340 COMPLEX*16 HL( NL, NL ), WORKL( NL ) 00341 * .. 00342 * .. Local Scalars .. 00343 INTEGER KBOT, NMIN 00344 LOGICAL INITZ, LQUERY, WANTT, WANTZ 00345 * .. 00346 * .. External Functions .. 00347 INTEGER ILAENV 00348 LOGICAL LSAME 00349 EXTERNAL ILAENV, LSAME 00350 * .. 00351 * .. External Subroutines .. 00352 EXTERNAL XERBLA, ZCOPY, ZLACPY, ZLAHQR, ZLAQR0, ZLASET 00353 * .. 00354 * .. Intrinsic Functions .. 00355 INTRINSIC DBLE, DCMPLX, MAX, MIN 00356 * .. 00357 * .. Executable Statements .. 00358 * 00359 * ==== Decode and check the input parameters. ==== 00360 * 00361 WANTT = LSAME( JOB, 'S' ) 00362 INITZ = LSAME( COMPZ, 'I' ) 00363 WANTZ = INITZ .OR. LSAME( COMPZ, 'V' ) 00364 WORK( 1 ) = DCMPLX( DBLE( MAX( 1, N ) ), RZERO ) 00365 LQUERY = LWORK.EQ.-1 00366 * 00367 INFO = 0 00368 IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN 00369 INFO = -1 00370 ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN 00371 INFO = -2 00372 ELSE IF( N.LT.0 ) THEN 00373 INFO = -3 00374 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN 00375 INFO = -4 00376 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN 00377 INFO = -5 00378 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN 00379 INFO = -7 00380 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN 00381 INFO = -10 00382 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00383 INFO = -12 00384 END IF 00385 * 00386 IF( INFO.NE.0 ) THEN 00387 * 00388 * ==== Quick return in case of invalid argument. ==== 00389 * 00390 CALL XERBLA( 'ZHSEQR', -INFO ) 00391 RETURN 00392 * 00393 ELSE IF( N.EQ.0 ) THEN 00394 * 00395 * ==== Quick return in case N = 0; nothing to do. ==== 00396 * 00397 RETURN 00398 * 00399 ELSE IF( LQUERY ) THEN 00400 * 00401 * ==== Quick return in case of a workspace query ==== 00402 * 00403 CALL ZLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI, Z, 00404 $ LDZ, WORK, LWORK, INFO ) 00405 * ==== Ensure reported workspace size is backward-compatible with 00406 * . previous LAPACK versions. ==== 00407 WORK( 1 ) = DCMPLX( MAX( DBLE( WORK( 1 ) ), DBLE( MAX( 1, 00408 $ N ) ) ), RZERO ) 00409 RETURN 00410 * 00411 ELSE 00412 * 00413 * ==== copy eigenvalues isolated by ZGEBAL ==== 00414 * 00415 IF( ILO.GT.1 ) 00416 $ CALL ZCOPY( ILO-1, H, LDH+1, W, 1 ) 00417 IF( IHI.LT.N ) 00418 $ CALL ZCOPY( N-IHI, H( IHI+1, IHI+1 ), LDH+1, W( IHI+1 ), 1 ) 00419 * 00420 * ==== Initialize Z, if requested ==== 00421 * 00422 IF( INITZ ) 00423 $ CALL ZLASET( 'A', N, N, ZERO, ONE, Z, LDZ ) 00424 * 00425 * ==== Quick return if possible ==== 00426 * 00427 IF( ILO.EQ.IHI ) THEN 00428 W( ILO ) = H( ILO, ILO ) 00429 RETURN 00430 END IF 00431 * 00432 * ==== ZLAHQR/ZLAQR0 crossover point ==== 00433 * 00434 NMIN = ILAENV( 12, 'ZHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N, 00435 $ ILO, IHI, LWORK ) 00436 NMIN = MAX( NTINY, NMIN ) 00437 * 00438 * ==== ZLAQR0 for big matrices; ZLAHQR for small ones ==== 00439 * 00440 IF( N.GT.NMIN ) THEN 00441 CALL ZLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI, 00442 $ Z, LDZ, WORK, LWORK, INFO ) 00443 ELSE 00444 * 00445 * ==== Small matrix ==== 00446 * 00447 CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI, 00448 $ Z, LDZ, INFO ) 00449 * 00450 IF( INFO.GT.0 ) THEN 00451 * 00452 * ==== A rare ZLAHQR failure! ZLAQR0 sometimes succeeds 00453 * . when ZLAHQR fails. ==== 00454 * 00455 KBOT = INFO 00456 * 00457 IF( N.GE.NL ) THEN 00458 * 00459 * ==== Larger matrices have enough subdiagonal scratch 00460 * . space to call ZLAQR0 directly. ==== 00461 * 00462 CALL ZLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, W, 00463 $ ILO, IHI, Z, LDZ, WORK, LWORK, INFO ) 00464 * 00465 ELSE 00466 * 00467 * ==== Tiny matrices don't have enough subdiagonal 00468 * . scratch space to benefit from ZLAQR0. Hence, 00469 * . tiny matrices must be copied into a larger 00470 * . array before calling ZLAQR0. ==== 00471 * 00472 CALL ZLACPY( 'A', N, N, H, LDH, HL, NL ) 00473 HL( N+1, N ) = ZERO 00474 CALL ZLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ), 00475 $ NL ) 00476 CALL ZLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, W, 00477 $ ILO, IHI, Z, LDZ, WORKL, NL, INFO ) 00478 IF( WANTT .OR. INFO.NE.0 ) 00479 $ CALL ZLACPY( 'A', N, N, HL, NL, H, LDH ) 00480 END IF 00481 END IF 00482 END IF 00483 * 00484 * ==== Clear out the trash, if necessary. ==== 00485 * 00486 IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 ) 00487 $ CALL ZLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH ) 00488 * 00489 * ==== Ensure reported workspace size is backward-compatible with 00490 * . previous LAPACK versions. ==== 00491 * 00492 WORK( 1 ) = DCMPLX( MAX( DBLE( MAX( 1, N ) ), 00493 $ DBLE( WORK( 1 ) ) ), RZERO ) 00494 END IF 00495 * 00496 * ==== End of ZHSEQR ==== 00497 * 00498 END