001:       SUBROUTINE SHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
002:      $                   LDZ, WORK, LWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
010:       CHARACTER          COMPZ, JOB
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               H( LDH, * ), WI( * ), WORK( * ), WR( * ),
014:      $                   Z( LDZ, * )
015: *     ..
016: *     Purpose
017: *     =======
018: *
019: *     SHSEQR computes the eigenvalues of a Hessenberg matrix H
020: *     and, optionally, the matrices T and Z from the Schur decomposition
021: *     H = Z T Z**T, where T is an upper quasi-triangular matrix (the
022: *     Schur form), and Z is the orthogonal matrix of Schur vectors.
023: *
024: *     Optionally Z may be postmultiplied into an input orthogonal
025: *     matrix Q so that this routine can give the Schur factorization
026: *     of a matrix A which has been reduced to the Hessenberg form H
027: *     by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
028: *
029: *     Arguments
030: *     =========
031: *
032: *     JOB   (input) CHARACTER*1
033: *           = 'E':  compute eigenvalues only;
034: *           = 'S':  compute eigenvalues and the Schur form T.
035: *
036: *     COMPZ (input) CHARACTER*1
037: *           = 'N':  no Schur vectors are computed;
038: *           = 'I':  Z is initialized to the unit matrix and the matrix Z
039: *                   of Schur vectors of H is returned;
040: *           = 'V':  Z must contain an orthogonal matrix Q on entry, and
041: *                   the product Q*Z is returned.
042: *
043: *     N     (input) INTEGER
044: *           The order of the matrix H.  N .GE. 0.
045: *
046: *     ILO   (input) INTEGER
047: *     IHI   (input) INTEGER
048: *           It is assumed that H is already upper triangular in rows
049: *           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
050: *           set by a previous call to SGEBAL, and then passed to SGEHRD
051: *           when the matrix output by SGEBAL is reduced to Hessenberg
052: *           form. Otherwise ILO and IHI should be set to 1 and N
053: *           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
054: *           If N = 0, then ILO = 1 and IHI = 0.
055: *
056: *     H     (input/output) REAL array, dimension (LDH,N)
057: *           On entry, the upper Hessenberg matrix H.
058: *           On exit, if INFO = 0 and JOB = 'S', then H contains the
059: *           upper quasi-triangular matrix T from the Schur decomposition
060: *           (the Schur form); 2-by-2 diagonal blocks (corresponding to
061: *           complex conjugate pairs of eigenvalues) are returned in
062: *           standard form, with H(i,i) = H(i+1,i+1) and
063: *           H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
064: *           contents of H are unspecified on exit.  (The output value of
065: *           H when INFO.GT.0 is given under the description of INFO
066: *           below.)
067: *
068: *           Unlike earlier versions of SHSEQR, this subroutine may
069: *           explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
070: *           or j = IHI+1, IHI+2, ... N.
071: *
072: *     LDH   (input) INTEGER
073: *           The leading dimension of the array H. LDH .GE. max(1,N).
074: *
075: *     WR    (output) REAL array, dimension (N)
076: *     WI    (output) REAL array, dimension (N)
077: *           The real and imaginary parts, respectively, of the computed
078: *           eigenvalues. If two eigenvalues are computed as a complex
079: *           conjugate pair, they are stored in consecutive elements of
080: *           WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
081: *           WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
082: *           the same order as on the diagonal of the Schur form returned
083: *           in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
084: *           diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
085: *           WI(i+1) = -WI(i).
086: *
087: *     Z     (input/output) REAL array, dimension (LDZ,N)
088: *           If COMPZ = 'N', Z is not referenced.
089: *           If COMPZ = 'I', on entry Z need not be set and on exit,
090: *           if INFO = 0, Z contains the orthogonal matrix Z of the Schur
091: *           vectors of H.  If COMPZ = 'V', on entry Z must contain an
092: *           N-by-N matrix Q, which is assumed to be equal to the unit
093: *           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
094: *           if INFO = 0, Z contains Q*Z.
095: *           Normally Q is the orthogonal matrix generated by SORGHR
096: *           after the call to SGEHRD which formed the Hessenberg matrix
097: *           H. (The output value of Z when INFO.GT.0 is given under
098: *           the description of INFO below.)
099: *
100: *     LDZ   (input) INTEGER
101: *           The leading dimension of the array Z.  if COMPZ = 'I' or
102: *           COMPZ = 'V', then LDZ.GE.MAX(1,N).  Otherwize, LDZ.GE.1.
103: *
104: *     WORK  (workspace/output) REAL array, dimension (LWORK)
105: *           On exit, if INFO = 0, WORK(1) returns an estimate of
106: *           the optimal value for LWORK.
107: *
108: *     LWORK (input) INTEGER
109: *           The dimension of the array WORK.  LWORK .GE. max(1,N)
110: *           is sufficient and delivers very good and sometimes
111: *           optimal performance.  However, LWORK as large as 11*N
112: *           may be required for optimal performance.  A workspace
113: *           query is recommended to determine the optimal workspace
114: *           size.
115: *
116: *           If LWORK = -1, then SHSEQR does a workspace query.
117: *           In this case, SHSEQR checks the input parameters and
118: *           estimates the optimal workspace size for the given
119: *           values of N, ILO and IHI.  The estimate is returned
120: *           in WORK(1).  No error message related to LWORK is
121: *           issued by XERBLA.  Neither H nor Z are accessed.
122: *
123: *
124: *     INFO  (output) INTEGER
125: *             =  0:  successful exit
126: *           .LT. 0:  if INFO = -i, the i-th argument had an illegal
127: *                    value
128: *           .GT. 0:  if INFO = i, SHSEQR failed to compute all of
129: *                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
130: *                and WI contain those eigenvalues which have been
131: *                successfully computed.  (Failures are rare.)
132: *
133: *                If INFO .GT. 0 and JOB = 'E', then on exit, the
134: *                remaining unconverged eigenvalues are the eigen-
135: *                values of the upper Hessenberg matrix rows and
136: *                columns ILO through INFO of the final, output
137: *                value of H.
138: *
139: *                If INFO .GT. 0 and JOB   = 'S', then on exit
140: *
141: *           (*)  (initial value of H)*U  = U*(final value of H)
142: *
143: *                where U is an orthogonal matrix.  The final
144: *                value of H is upper Hessenberg and quasi-triangular
145: *                in rows and columns INFO+1 through IHI.
146: *
147: *                If INFO .GT. 0 and COMPZ = 'V', then on exit
148: *
149: *                  (final value of Z)  =  (initial value of Z)*U
150: *
151: *                where U is the orthogonal matrix in (*) (regard-
152: *                less of the value of JOB.)
153: *
154: *                If INFO .GT. 0 and COMPZ = 'I', then on exit
155: *                      (final value of Z)  = U
156: *                where U is the orthogonal matrix in (*) (regard-
157: *                less of the value of JOB.)
158: *
159: *                If INFO .GT. 0 and COMPZ = 'N', then Z is not
160: *                accessed.
161: *
162: *     ================================================================
163: *             Default values supplied by
164: *             ILAENV(ISPEC,'SHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
165: *             It is suggested that these defaults be adjusted in order
166: *             to attain best performance in each particular
167: *             computational environment.
168: *
169: *            ISPEC=12: The SLAHQR vs SLAQR0 crossover point.
170: *                      Default: 75. (Must be at least 11.)
171: *
172: *            ISPEC=13: Recommended deflation window size.
173: *                      This depends on ILO, IHI and NS.  NS is the
174: *                      number of simultaneous shifts returned
175: *                      by ILAENV(ISPEC=15).  (See ISPEC=15 below.)
176: *                      The default for (IHI-ILO+1).LE.500 is NS.
177: *                      The default for (IHI-ILO+1).GT.500 is 3*NS/2.
178: *
179: *            ISPEC=14: Nibble crossover point. (See IPARMQ for
180: *                      details.)  Default: 14% of deflation window
181: *                      size.
182: *
183: *            ISPEC=15: Number of simultaneous shifts in a multishift
184: *                      QR iteration.
185: *
186: *                      If IHI-ILO+1 is ...
187: *
188: *                      greater than      ...but less    ... the
189: *                      or equal to ...      than        default is
190: *
191: *                           1               30          NS =   2(+)
192: *                          30               60          NS =   4(+)
193: *                          60              150          NS =  10(+)
194: *                         150              590          NS =  **
195: *                         590             3000          NS =  64
196: *                        3000             6000          NS = 128
197: *                        6000             infinity      NS = 256
198: *
199: *                  (+)  By default some or all matrices of this order
200: *                       are passed to the implicit double shift routine
201: *                       SLAHQR and this parameter is ignored.  See
202: *                       ISPEC=12 above and comments in IPARMQ for
203: *                       details.
204: *
205: *                 (**)  The asterisks (**) indicate an ad-hoc
206: *                       function of N increasing from 10 to 64.
207: *
208: *            ISPEC=16: Select structured matrix multiply.
209: *                      If the number of simultaneous shifts (specified
210: *                      by ISPEC=15) is less than 14, then the default
211: *                      for ISPEC=16 is 0.  Otherwise the default for
212: *                      ISPEC=16 is 2.
213: *
214: *     ================================================================
215: *     Based on contributions by
216: *        Karen Braman and Ralph Byers, Department of Mathematics,
217: *        University of Kansas, USA
218: *
219: *     ================================================================
220: *     References:
221: *       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
222: *       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
223: *       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
224: *       929--947, 2002.
225: *
226: *       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
227: *       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
228: *       of Matrix Analysis, volume 23, pages 948--973, 2002.
229: *
230: *     ================================================================
231: *     .. Parameters ..
232: *
233: *     ==== Matrices of order NTINY or smaller must be processed by
234: *     .    SLAHQR because of insufficient subdiagonal scratch space.
235: *     .    (This is a hard limit.) ====
236:       INTEGER            NTINY
237:       PARAMETER          ( NTINY = 11 )
238: *
239: *     ==== NL allocates some local workspace to help small matrices
240: *     .    through a rare SLAHQR failure.  NL .GT. NTINY = 11 is
241: *     .    required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
242: *     .    mended.  (The default value of NMIN is 75.)  Using NL = 49
243: *     .    allows up to six simultaneous shifts and a 16-by-16
244: *     .    deflation window.  ====
245:       INTEGER            NL
246:       PARAMETER          ( NL = 49 )
247:       REAL               ZERO, ONE
248:       PARAMETER          ( ZERO = 0.0e0, ONE = 1.0e0 )
249: *     ..
250: *     .. Local Arrays ..
251:       REAL               HL( NL, NL ), WORKL( NL )
252: *     ..
253: *     .. Local Scalars ..
254:       INTEGER            I, KBOT, NMIN
255:       LOGICAL            INITZ, LQUERY, WANTT, WANTZ
256: *     ..
257: *     .. External Functions ..
258:       INTEGER            ILAENV
259:       LOGICAL            LSAME
260:       EXTERNAL           ILAENV, LSAME
261: *     ..
262: *     .. External Subroutines ..
263:       EXTERNAL           SLACPY, SLAHQR, SLAQR0, SLASET, XERBLA
264: *     ..
265: *     .. Intrinsic Functions ..
266:       INTRINSIC          MAX, MIN, REAL
267: *     ..
268: *     .. Executable Statements ..
269: *
270: *     ==== Decode and check the input parameters. ====
271: *
272:       WANTT = LSAME( JOB, 'S' )
273:       INITZ = LSAME( COMPZ, 'I' )
274:       WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
275:       WORK( 1 ) = REAL( MAX( 1, N ) )
276:       LQUERY = LWORK.EQ.-1
277: *
278:       INFO = 0
279:       IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
280:          INFO = -1
281:       ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
282:          INFO = -2
283:       ELSE IF( N.LT.0 ) THEN
284:          INFO = -3
285:       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
286:          INFO = -4
287:       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
288:          INFO = -5
289:       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
290:          INFO = -7
291:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
292:          INFO = -11
293:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
294:          INFO = -13
295:       END IF
296: *
297:       IF( INFO.NE.0 ) THEN
298: *
299: *        ==== Quick return in case of invalid argument. ====
300: *
301:          CALL XERBLA( 'SHSEQR', -INFO )
302:          RETURN
303: *
304:       ELSE IF( N.EQ.0 ) THEN
305: *
306: *        ==== Quick return in case N = 0; nothing to do. ====
307: *
308:          RETURN
309: *
310:       ELSE IF( LQUERY ) THEN
311: *
312: *        ==== Quick return in case of a workspace query ====
313: *
314:          CALL SLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
315:      $                IHI, Z, LDZ, WORK, LWORK, INFO )
316: *        ==== Ensure reported workspace size is backward-compatible with
317: *        .    previous LAPACK versions. ====
318:          WORK( 1 ) = MAX( REAL( MAX( 1, N ) ), WORK( 1 ) )
319:          RETURN
320: *
321:       ELSE
322: *
323: *        ==== copy eigenvalues isolated by SGEBAL ====
324: *
325:          DO 10 I = 1, ILO - 1
326:             WR( I ) = H( I, I )
327:             WI( I ) = ZERO
328:    10    CONTINUE
329:          DO 20 I = IHI + 1, N
330:             WR( I ) = H( I, I )
331:             WI( I ) = ZERO
332:    20    CONTINUE
333: *
334: *        ==== Initialize Z, if requested ====
335: *
336:          IF( INITZ )
337:      $      CALL SLASET( 'A', N, N, ZERO, ONE, Z, LDZ )
338: *
339: *        ==== Quick return if possible ====
340: *
341:          IF( ILO.EQ.IHI ) THEN
342:             WR( ILO ) = H( ILO, ILO )
343:             WI( ILO ) = ZERO
344:             RETURN
345:          END IF
346: *
347: *        ==== SLAHQR/SLAQR0 crossover point ====
348: *
349:          NMIN = ILAENV( 12, 'SHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N,
350:      $          ILO, IHI, LWORK )
351:          NMIN = MAX( NTINY, NMIN )
352: *
353: *        ==== SLAQR0 for big matrices; SLAHQR for small ones ====
354: *
355:          IF( N.GT.NMIN ) THEN
356:             CALL SLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
357:      $                   IHI, Z, LDZ, WORK, LWORK, INFO )
358:          ELSE
359: *
360: *           ==== Small matrix ====
361: *
362:             CALL SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
363:      $                   IHI, Z, LDZ, INFO )
364: *
365:             IF( INFO.GT.0 ) THEN
366: *
367: *              ==== A rare SLAHQR failure!  SLAQR0 sometimes succeeds
368: *              .    when SLAHQR fails. ====
369: *
370:                KBOT = INFO
371: *
372:                IF( N.GE.NL ) THEN
373: *
374: *                 ==== Larger matrices have enough subdiagonal scratch
375: *                 .    space to call SLAQR0 directly. ====
376: *
377:                   CALL SLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, WR,
378:      $                         WI, ILO, IHI, Z, LDZ, WORK, LWORK, INFO )
379: *
380:                ELSE
381: *
382: *                 ==== Tiny matrices don't have enough subdiagonal
383: *                 .    scratch space to benefit from SLAQR0.  Hence,
384: *                 .    tiny matrices must be copied into a larger
385: *                 .    array before calling SLAQR0. ====
386: *
387:                   CALL SLACPY( 'A', N, N, H, LDH, HL, NL )
388:                   HL( N+1, N ) = ZERO
389:                   CALL SLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ),
390:      $                         NL )
391:                   CALL SLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, WR,
392:      $                         WI, ILO, IHI, Z, LDZ, WORKL, NL, INFO )
393:                   IF( WANTT .OR. INFO.NE.0 )
394:      $               CALL SLACPY( 'A', N, N, HL, NL, H, LDH )
395:                END IF
396:             END IF
397:          END IF
398: *
399: *        ==== Clear out the trash, if necessary. ====
400: *
401:          IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 )
402:      $      CALL SLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH )
403: *
404: *        ==== Ensure reported workspace size is backward-compatible with
405: *        .    previous LAPACK versions. ====
406: *
407:          WORK( 1 ) = MAX( REAL( MAX( 1, N ) ), WORK( 1 ) )
408:       END IF
409: *
410: *     ==== End of SHSEQR ====
411: *
412:       END
413: