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