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