001:       SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
002:      $                   LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            NOINIT, RIGHTV
011:       INTEGER            INFO, LDB, LDH, N
012:       DOUBLE PRECISION   BIGNUM, EPS3, SMLNUM, WI, WR
013: *     ..
014: *     .. Array Arguments ..
015:       DOUBLE PRECISION   B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
016:      $                   WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  DLAEIN uses inverse iteration to find a right or left eigenvector
023: *  corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
024: *  matrix H.
025: *
026: *  Arguments
027: *  =========
028: *
029: *  RIGHTV   (input) LOGICAL
030: *          = .TRUE. : compute right eigenvector;
031: *          = .FALSE.: compute left eigenvector.
032: *
033: *  NOINIT   (input) LOGICAL
034: *          = .TRUE. : no initial vector supplied in (VR,VI).
035: *          = .FALSE.: initial vector supplied in (VR,VI).
036: *
037: *  N       (input) INTEGER
038: *          The order of the matrix H.  N >= 0.
039: *
040: *  H       (input) DOUBLE PRECISION array, dimension (LDH,N)
041: *          The upper Hessenberg matrix H.
042: *
043: *  LDH     (input) INTEGER
044: *          The leading dimension of the array H.  LDH >= max(1,N).
045: *
046: *  WR      (input) DOUBLE PRECISION
047: *  WI      (input) DOUBLE PRECISION
048: *          The real and imaginary parts of the eigenvalue of H whose
049: *          corresponding right or left eigenvector is to be computed.
050: *
051: *  VR      (input/output) DOUBLE PRECISION array, dimension (N)
052: *  VI      (input/output) DOUBLE PRECISION array, dimension (N)
053: *          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
054: *          a real starting vector for inverse iteration using the real
055: *          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
056: *          must contain the real and imaginary parts of a complex
057: *          starting vector for inverse iteration using the complex
058: *          eigenvalue (WR,WI); otherwise VR and VI need not be set.
059: *          On exit, if WI = 0.0 (real eigenvalue), VR contains the
060: *          computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
061: *          VR and VI contain the real and imaginary parts of the
062: *          computed complex eigenvector. The eigenvector is normalized
063: *          so that the component of largest magnitude has magnitude 1;
064: *          here the magnitude of a complex number (x,y) is taken to be
065: *          |x| + |y|.
066: *          VI is not referenced if WI = 0.0.
067: *
068: *  B       (workspace) DOUBLE PRECISION array, dimension (LDB,N)
069: *
070: *  LDB     (input) INTEGER
071: *          The leading dimension of the array B.  LDB >= N+1.
072: *
073: *  WORK   (workspace) DOUBLE PRECISION array, dimension (N)
074: *
075: *  EPS3    (input) DOUBLE PRECISION
076: *          A small machine-dependent value which is used to perturb
077: *          close eigenvalues, and to replace zero pivots.
078: *
079: *  SMLNUM  (input) DOUBLE PRECISION
080: *          A machine-dependent value close to the underflow threshold.
081: *
082: *  BIGNUM  (input) DOUBLE PRECISION
083: *          A machine-dependent value close to the overflow threshold.
084: *
085: *  INFO    (output) INTEGER
086: *          = 0:  successful exit
087: *          = 1:  inverse iteration did not converge; VR is set to the
088: *                last iterate, and so is VI if WI.ne.0.0.
089: *
090: *  =====================================================================
091: *
092: *     .. Parameters ..
093:       DOUBLE PRECISION   ZERO, ONE, TENTH
094:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
095: *     ..
096: *     .. Local Scalars ..
097:       CHARACTER          NORMIN, TRANS
098:       INTEGER            I, I1, I2, I3, IERR, ITS, J
099:       DOUBLE PRECISION   ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
100:      $                   REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
101:      $                   W1, X, XI, XR, Y
102: *     ..
103: *     .. External Functions ..
104:       INTEGER            IDAMAX
105:       DOUBLE PRECISION   DASUM, DLAPY2, DNRM2
106:       EXTERNAL           IDAMAX, DASUM, DLAPY2, DNRM2
107: *     ..
108: *     .. External Subroutines ..
109:       EXTERNAL           DLADIV, DLATRS, DSCAL
110: *     ..
111: *     .. Intrinsic Functions ..
112:       INTRINSIC          ABS, DBLE, MAX, SQRT
113: *     ..
114: *     .. Executable Statements ..
115: *
116:       INFO = 0
117: *
118: *     GROWTO is the threshold used in the acceptance test for an
119: *     eigenvector.
120: *
121:       ROOTN = SQRT( DBLE( N ) )
122:       GROWTO = TENTH / ROOTN
123:       NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
124: *
125: *     Form B = H - (WR,WI)*I (except that the subdiagonal elements and
126: *     the imaginary parts of the diagonal elements are not stored).
127: *
128:       DO 20 J = 1, N
129:          DO 10 I = 1, J - 1
130:             B( I, J ) = H( I, J )
131:    10    CONTINUE
132:          B( J, J ) = H( J, J ) - WR
133:    20 CONTINUE
134: *
135:       IF( WI.EQ.ZERO ) THEN
136: *
137: *        Real eigenvalue.
138: *
139:          IF( NOINIT ) THEN
140: *
141: *           Set initial vector.
142: *
143:             DO 30 I = 1, N
144:                VR( I ) = EPS3
145:    30       CONTINUE
146:          ELSE
147: *
148: *           Scale supplied initial vector.
149: *
150:             VNORM = DNRM2( N, VR, 1 )
151:             CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
152:      $                  1 )
153:          END IF
154: *
155:          IF( RIGHTV ) THEN
156: *
157: *           LU decomposition with partial pivoting of B, replacing zero
158: *           pivots by EPS3.
159: *
160:             DO 60 I = 1, N - 1
161:                EI = H( I+1, I )
162:                IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
163: *
164: *                 Interchange rows and eliminate.
165: *
166:                   X = B( I, I ) / EI
167:                   B( I, I ) = EI
168:                   DO 40 J = I + 1, N
169:                      TEMP = B( I+1, J )
170:                      B( I+1, J ) = B( I, J ) - X*TEMP
171:                      B( I, J ) = TEMP
172:    40             CONTINUE
173:                ELSE
174: *
175: *                 Eliminate without interchange.
176: *
177:                   IF( B( I, I ).EQ.ZERO )
178:      $               B( I, I ) = EPS3
179:                   X = EI / B( I, I )
180:                   IF( X.NE.ZERO ) THEN
181:                      DO 50 J = I + 1, N
182:                         B( I+1, J ) = B( I+1, J ) - X*B( I, J )
183:    50                CONTINUE
184:                   END IF
185:                END IF
186:    60       CONTINUE
187:             IF( B( N, N ).EQ.ZERO )
188:      $         B( N, N ) = EPS3
189: *
190:             TRANS = 'N'
191: *
192:          ELSE
193: *
194: *           UL decomposition with partial pivoting of B, replacing zero
195: *           pivots by EPS3.
196: *
197:             DO 90 J = N, 2, -1
198:                EJ = H( J, J-1 )
199:                IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
200: *
201: *                 Interchange columns and eliminate.
202: *
203:                   X = B( J, J ) / EJ
204:                   B( J, J ) = EJ
205:                   DO 70 I = 1, J - 1
206:                      TEMP = B( I, J-1 )
207:                      B( I, J-1 ) = B( I, J ) - X*TEMP
208:                      B( I, J ) = TEMP
209:    70             CONTINUE
210:                ELSE
211: *
212: *                 Eliminate without interchange.
213: *
214:                   IF( B( J, J ).EQ.ZERO )
215:      $               B( J, J ) = EPS3
216:                   X = EJ / B( J, J )
217:                   IF( X.NE.ZERO ) THEN
218:                      DO 80 I = 1, J - 1
219:                         B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
220:    80                CONTINUE
221:                   END IF
222:                END IF
223:    90       CONTINUE
224:             IF( B( 1, 1 ).EQ.ZERO )
225:      $         B( 1, 1 ) = EPS3
226: *
227:             TRANS = 'T'
228: *
229:          END IF
230: *
231:          NORMIN = 'N'
232:          DO 110 ITS = 1, N
233: *
234: *           Solve U*x = scale*v for a right eigenvector
235: *             or U'*x = scale*v for a left eigenvector,
236: *           overwriting x on v.
237: *
238:             CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
239:      $                   VR, SCALE, WORK, IERR )
240:             NORMIN = 'Y'
241: *
242: *           Test for sufficient growth in the norm of v.
243: *
244:             VNORM = DASUM( N, VR, 1 )
245:             IF( VNORM.GE.GROWTO*SCALE )
246:      $         GO TO 120
247: *
248: *           Choose new orthogonal starting vector and try again.
249: *
250:             TEMP = EPS3 / ( ROOTN+ONE )
251:             VR( 1 ) = EPS3
252:             DO 100 I = 2, N
253:                VR( I ) = TEMP
254:   100       CONTINUE
255:             VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
256:   110    CONTINUE
257: *
258: *        Failure to find eigenvector in N iterations.
259: *
260:          INFO = 1
261: *
262:   120    CONTINUE
263: *
264: *        Normalize eigenvector.
265: *
266:          I = IDAMAX( N, VR, 1 )
267:          CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
268:       ELSE
269: *
270: *        Complex eigenvalue.
271: *
272:          IF( NOINIT ) THEN
273: *
274: *           Set initial vector.
275: *
276:             DO 130 I = 1, N
277:                VR( I ) = EPS3
278:                VI( I ) = ZERO
279:   130       CONTINUE
280:          ELSE
281: *
282: *           Scale supplied initial vector.
283: *
284:             NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
285:             REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
286:             CALL DSCAL( N, REC, VR, 1 )
287:             CALL DSCAL( N, REC, VI, 1 )
288:          END IF
289: *
290:          IF( RIGHTV ) THEN
291: *
292: *           LU decomposition with partial pivoting of B, replacing zero
293: *           pivots by EPS3.
294: *
295: *           The imaginary part of the (i,j)-th element of U is stored in
296: *           B(j+1,i).
297: *
298:             B( 2, 1 ) = -WI
299:             DO 140 I = 2, N
300:                B( I+1, 1 ) = ZERO
301:   140       CONTINUE
302: *
303:             DO 170 I = 1, N - 1
304:                ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
305:                EI = H( I+1, I )
306:                IF( ABSBII.LT.ABS( EI ) ) THEN
307: *
308: *                 Interchange rows and eliminate.
309: *
310:                   XR = B( I, I ) / EI
311:                   XI = B( I+1, I ) / EI
312:                   B( I, I ) = EI
313:                   B( I+1, I ) = ZERO
314:                   DO 150 J = I + 1, N
315:                      TEMP = B( I+1, J )
316:                      B( I+1, J ) = B( I, J ) - XR*TEMP
317:                      B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
318:                      B( I, J ) = TEMP
319:                      B( J+1, I ) = ZERO
320:   150             CONTINUE
321:                   B( I+2, I ) = -WI
322:                   B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
323:                   B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
324:                ELSE
325: *
326: *                 Eliminate without interchanging rows.
327: *
328:                   IF( ABSBII.EQ.ZERO ) THEN
329:                      B( I, I ) = EPS3
330:                      B( I+1, I ) = ZERO
331:                      ABSBII = EPS3
332:                   END IF
333:                   EI = ( EI / ABSBII ) / ABSBII
334:                   XR = B( I, I )*EI
335:                   XI = -B( I+1, I )*EI
336:                   DO 160 J = I + 1, N
337:                      B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
338:      $                             XI*B( J+1, I )
339:                      B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
340:   160             CONTINUE
341:                   B( I+2, I+1 ) = B( I+2, I+1 ) - WI
342:                END IF
343: *
344: *              Compute 1-norm of offdiagonal elements of i-th row.
345: *
346:                WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
347:      $                     DASUM( N-I, B( I+2, I ), 1 )
348:   170       CONTINUE
349:             IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
350:      $         B( N, N ) = EPS3
351:             WORK( N ) = ZERO
352: *
353:             I1 = N
354:             I2 = 1
355:             I3 = -1
356:          ELSE
357: *
358: *           UL decomposition with partial pivoting of conjg(B),
359: *           replacing zero pivots by EPS3.
360: *
361: *           The imaginary part of the (i,j)-th element of U is stored in
362: *           B(j+1,i).
363: *
364:             B( N+1, N ) = WI
365:             DO 180 J = 1, N - 1
366:                B( N+1, J ) = ZERO
367:   180       CONTINUE
368: *
369:             DO 210 J = N, 2, -1
370:                EJ = H( J, J-1 )
371:                ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
372:                IF( ABSBJJ.LT.ABS( EJ ) ) THEN
373: *
374: *                 Interchange columns and eliminate
375: *
376:                   XR = B( J, J ) / EJ
377:                   XI = B( J+1, J ) / EJ
378:                   B( J, J ) = EJ
379:                   B( J+1, J ) = ZERO
380:                   DO 190 I = 1, J - 1
381:                      TEMP = B( I, J-1 )
382:                      B( I, J-1 ) = B( I, J ) - XR*TEMP
383:                      B( J, I ) = B( J+1, I ) - XI*TEMP
384:                      B( I, J ) = TEMP
385:                      B( J+1, I ) = ZERO
386:   190             CONTINUE
387:                   B( J+1, J-1 ) = WI
388:                   B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
389:                   B( J, J-1 ) = B( J, J-1 ) - XR*WI
390:                ELSE
391: *
392: *                 Eliminate without interchange.
393: *
394:                   IF( ABSBJJ.EQ.ZERO ) THEN
395:                      B( J, J ) = EPS3
396:                      B( J+1, J ) = ZERO
397:                      ABSBJJ = EPS3
398:                   END IF
399:                   EJ = ( EJ / ABSBJJ ) / ABSBJJ
400:                   XR = B( J, J )*EJ
401:                   XI = -B( J+1, J )*EJ
402:                   DO 200 I = 1, J - 1
403:                      B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
404:      $                             XI*B( J+1, I )
405:                      B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
406:   200             CONTINUE
407:                   B( J, J-1 ) = B( J, J-1 ) + WI
408:                END IF
409: *
410: *              Compute 1-norm of offdiagonal elements of j-th column.
411: *
412:                WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
413:      $                     DASUM( J-1, B( J+1, 1 ), LDB )
414:   210       CONTINUE
415:             IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
416:      $         B( 1, 1 ) = EPS3
417:             WORK( 1 ) = ZERO
418: *
419:             I1 = 1
420:             I2 = N
421:             I3 = 1
422:          END IF
423: *
424:          DO 270 ITS = 1, N
425:             SCALE = ONE
426:             VMAX = ONE
427:             VCRIT = BIGNUM
428: *
429: *           Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
430: *             or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector,
431: *           overwriting (xr,xi) on (vr,vi).
432: *
433:             DO 250 I = I1, I2, I3
434: *
435:                IF( WORK( I ).GT.VCRIT ) THEN
436:                   REC = ONE / VMAX
437:                   CALL DSCAL( N, REC, VR, 1 )
438:                   CALL DSCAL( N, REC, VI, 1 )
439:                   SCALE = SCALE*REC
440:                   VMAX = ONE
441:                   VCRIT = BIGNUM
442:                END IF
443: *
444:                XR = VR( I )
445:                XI = VI( I )
446:                IF( RIGHTV ) THEN
447:                   DO 220 J = I + 1, N
448:                      XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
449:                      XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
450:   220             CONTINUE
451:                ELSE
452:                   DO 230 J = 1, I - 1
453:                      XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
454:                      XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
455:   230             CONTINUE
456:                END IF
457: *
458:                W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
459:                IF( W.GT.SMLNUM ) THEN
460:                   IF( W.LT.ONE ) THEN
461:                      W1 = ABS( XR ) + ABS( XI )
462:                      IF( W1.GT.W*BIGNUM ) THEN
463:                         REC = ONE / W1
464:                         CALL DSCAL( N, REC, VR, 1 )
465:                         CALL DSCAL( N, REC, VI, 1 )
466:                         XR = VR( I )
467:                         XI = VI( I )
468:                         SCALE = SCALE*REC
469:                         VMAX = VMAX*REC
470:                      END IF
471:                   END IF
472: *
473: *                 Divide by diagonal element of B.
474: *
475:                   CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
476:      $                         VI( I ) )
477:                   VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
478:                   VCRIT = BIGNUM / VMAX
479:                ELSE
480:                   DO 240 J = 1, N
481:                      VR( J ) = ZERO
482:                      VI( J ) = ZERO
483:   240             CONTINUE
484:                   VR( I ) = ONE
485:                   VI( I ) = ONE
486:                   SCALE = ZERO
487:                   VMAX = ONE
488:                   VCRIT = BIGNUM
489:                END IF
490:   250       CONTINUE
491: *
492: *           Test for sufficient growth in the norm of (VR,VI).
493: *
494:             VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
495:             IF( VNORM.GE.GROWTO*SCALE )
496:      $         GO TO 280
497: *
498: *           Choose a new orthogonal starting vector and try again.
499: *
500:             Y = EPS3 / ( ROOTN+ONE )
501:             VR( 1 ) = EPS3
502:             VI( 1 ) = ZERO
503: *
504:             DO 260 I = 2, N
505:                VR( I ) = Y
506:                VI( I ) = ZERO
507:   260       CONTINUE
508:             VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
509:   270    CONTINUE
510: *
511: *        Failure to find eigenvector in N iterations
512: *
513:          INFO = 1
514: *
515:   280    CONTINUE
516: *
517: *        Normalize eigenvector.
518: *
519:          VNORM = ZERO
520:          DO 290 I = 1, N
521:             VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
522:   290    CONTINUE
523:          CALL DSCAL( N, ONE / VNORM, VR, 1 )
524:          CALL DSCAL( N, ONE / VNORM, VI, 1 )
525: *
526:       END IF
527: *
528:       RETURN
529: *
530: *     End of DLAEIN
531: *
532:       END
533: