1      SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
2     $                   LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
3*
4*  -- LAPACK auxiliary routine (instrumented to count operations) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     September 30, 1994
8*
9*     .. Scalar Arguments ..
10      LOGICAL            NOINIT, RIGHTV
11      INTEGER            INFO, LDB, LDH, N
12      DOUBLE PRECISION   BIGNUM, EPS3, SMLNUM, WI, WR
13*     ..
14*     .. Array Arguments ..
15      DOUBLE PRECISION   B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
16     $                   WORK( * )
17*     ..
18*     Common block to return operation count.
19*     .. Common blocks ..
20      COMMON             / LATIME / OPS, ITCNT
21*     ..
22*     .. Scalars in Common ..
23      DOUBLE PRECISION   ITCNT, OPS
24*     ..
25*
26*  Purpose
27*  =======
28*
29*  DLAEIN uses inverse iteration to find a right or left eigenvector
30*  corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
31*  matrix H.
32*
33*  Arguments
34*  =========
35*
36*  RIGHTV   (input) LOGICAL
37*          = .TRUE. : compute right eigenvector;
38*          = .FALSE.: compute left eigenvector.
39*
40*  NOINIT   (input) LOGICAL
41*          = .TRUE. : no initial vector supplied in (VR,VI).
42*          = .FALSE.: initial vector supplied in (VR,VI).
43*
44*  N       (input) INTEGER
45*          The order of the matrix H.  N >= 0.
46*
47*  H       (input) DOUBLE PRECISION array, dimension (LDH,N)
48*          The upper Hessenberg matrix H.
49*
50*  LDH     (input) INTEGER
51*          The leading dimension of the array H.  LDH >= max(1,N).
52*
53*  WR      (input) DOUBLE PRECISION
54*  WI      (input) DOUBLE PRECISION
55*          The real and imaginary parts of the eigenvalue of H whose
56*          corresponding right or left eigenvector is to be computed.
57*
58*  VR      (input/output) DOUBLE PRECISION array, dimension (N)
59*  VI      (input/output) DOUBLE PRECISION array, dimension (N)
60*          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
61*          a real starting vector for inverse iteration using the real
62*          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
63*          must contain the real and imaginary parts of a complex
64*          starting vector for inverse iteration using the complex
65*          eigenvalue (WR,WI); otherwise VR and VI need not be set.
66*          On exit, if WI = 0.0 (real eigenvalue), VR contains the
67*          computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
68*          VR and VI contain the real and imaginary parts of the
69*          computed complex eigenvector. The eigenvector is normalized
70*          so that the component of largest magnitude has magnitude 1;
71*          here the magnitude of a complex number (x,y) is taken to be
72*          |x| + |y|.
73*          VI is not referenced if WI = 0.0.
74*
75*  B       (workspace) DOUBLE PRECISION array, dimension (LDB,N)
76*
77*  LDB     (input) INTEGER
78*          The leading dimension of the array B.  LDB >= N+1.
79*
80*  WORK   (workspace) DOUBLE PRECISION array, dimension (N)
81*
82*  EPS3    (input) DOUBLE PRECISION
83*          A small machine-dependent value which is used to perturb
84*          close eigenvalues, and to replace zero pivots.
85*
86*  SMLNUM  (input) DOUBLE PRECISION
87*          A machine-dependent value close to the underflow threshold.
88*
89*  BIGNUM  (input) DOUBLE PRECISION
90*          A machine-dependent value close to the overflow threshold.
91*
92*  INFO    (output) INTEGER
93*          = 0:  successful exit
94*          = 1:  inverse iteration did not converge; VR is set to the
95*                last iterate, and so is VI if WI.ne.0.0.
96*
97*  =====================================================================
98*
99*     .. Parameters ..
100      DOUBLE PRECISION   ZERO, ONE, TENTH
101      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
102*     ..
103*     .. Local Scalars ..
104      CHARACTER          NORMIN, TRANS
105      INTEGER            I, I1, I2, I3, IERR, ITS, J
106      DOUBLE PRECISION   ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
107     $                   OPST, REC, ROOTN, SCALE, TEMP, VCRIT, VMAX,
108     $                   VNORM, W, W1, X, XI, XR, Y
109*     ..
110*     .. External Functions ..
111      INTEGER            IDAMAX
112      DOUBLE PRECISION   DASUM, DLAPY2, DNRM2
113      EXTERNAL           IDAMAX, DASUM, DLAPY2, DNRM2
114*     ..
115*     .. External Subroutines ..
116      EXTERNAL           DLADIV, DLATRS, DSCAL
117*     ..
118*     .. Intrinsic Functions ..
119      INTRINSIC          ABS, DBLE, MAX, SQRT
120*     ..
121*     .. Executable Statements ..
122*
123      INFO = 0
124***
125*     Initialize
126      OPST = 0
127***
128*
129*     GROWTO is the threshold used in the acceptance test for an
130*     eigenvector.
131*
132      ROOTN = SQRT( DBLE( N ) )
133      GROWTO = TENTH / ROOTN
134      NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
135***
136*        Increment op count for computing ROOTN, GROWTO and NRMSML
137      OPST = OPST + 4
138***
139*
140*     Form B = H - (WR,WI)*I (except that the subdiagonal elements and
141*     the imaginary parts of the diagonal elements are not stored).
142*
143      DO 20 J = 1, N
144         DO 10 I = 1, J - 1
145            B( I, J ) = H( I, J )
146   10    CONTINUE
147         B( J, J ) = H( J, J ) - WR
148   20 CONTINUE
149***
150      OPST = OPST + N
151***
152*
153      IF( WI.EQ.ZERO ) THEN
154*
155*        Real eigenvalue.
156*
157         IF( NOINIT ) THEN
158*
159*           Set initial vector.
160*
161            DO 30 I = 1, N
162               VR( I ) = EPS3
163   30       CONTINUE
164         ELSE
165*
166*           Scale supplied initial vector.
167*
168            VNORM = DNRM2( N, VR, 1 )
169            CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
170     $                  1 )
171***
172            OPST = OPST + ( 3*N+2 )
173***
174         END IF
175*
176         IF( RIGHTV ) THEN
177*
178*           LU decomposition with partial pivoting of B, replacing zero
179*           pivots by EPS3.
180*
181            DO 60 I = 1, N - 1
182               EI = H( I+1, I )
183               IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
184*
185*                 Interchange rows and eliminate.
186*
187                  X = B( I, I ) / EI
188                  B( I, I ) = EI
189                  DO 40 J = I + 1, N
190                     TEMP = B( I+1, J )
191                     B( I+1, J ) = B( I, J ) - X*TEMP
192                     B( I, J ) = TEMP
193   40             CONTINUE
194               ELSE
195*
196*                 Eliminate without interchange.
197*
198                  IF( B( I, I ).EQ.ZERO )
199     $               B( I, I ) = EPS3
200                  X = EI / B( I, I )
201                  IF( X.NE.ZERO ) THEN
202                     DO 50 J = I + 1, N
203                        B( I+1, J ) = B( I+1, J ) - X*B( I, J )
204   50                CONTINUE
205                  END IF
206               END IF
207   60       CONTINUE
208            IF( B( N, N ).EQ.ZERO )
209     $         B( N, N ) = EPS3
210***
211*           Increment op count for LU decomposition
212            OPS = OPS + ( N-1 )*( N+1 )
213***
214*
215            TRANS = 'N'
216*
217         ELSE
218*
219*           UL decomposition with partial pivoting of B, replacing zero
220*           pivots by EPS3.
221*
222            DO 90 J = N, 2, -1
223               EJ = H( J, J-1 )
224               IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
225*
226*                 Interchange columns and eliminate.
227*
228                  X = B( J, J ) / EJ
229                  B( J, J ) = EJ
230                  DO 70 I = 1, J - 1
231                     TEMP = B( I, J-1 )
232                     B( I, J-1 ) = B( I, J ) - X*TEMP
233                     B( I, J ) = TEMP
234   70             CONTINUE
235               ELSE
236*
237*                 Eliminate without interchange.
238*
239                  IF( B( J, J ).EQ.ZERO )
240     $               B( J, J ) = EPS3
241                  X = EJ / B( J, J )
242                  IF( X.NE.ZERO ) THEN
243                     DO 80 I = 1, J - 1
244                        B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
245   80                CONTINUE
246                  END IF
247               END IF
248   90       CONTINUE
249            IF( B( 1, 1 ).EQ.ZERO )
250     $         B( 1, 1 ) = EPS3
251***
252*           Increment op count for UL decomposition
253            OPS = OPS + ( N-1 )*( N+1 )
254***
255*
256            TRANS = 'T'
257*
258         END IF
259*
260         NORMIN = 'N'
261         DO 110 ITS = 1, N
262*
263*           Solve U*x = scale*v for a right eigenvector
264*             or U'*x = scale*v for a left eigenvector,
265*           overwriting x on v.
266*
267            CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
268     $                   VR, SCALE, WORK, IERR )
269***
270*           Increment opcount for triangular solver, assuming that
271*           ops DLATRS = ops DTRSV, with no scaling in DLATRS.
272            OPS = OPS + N*N
273***
274            NORMIN = 'Y'
275*
276*           Test for sufficient growth in the norm of v.
277*
278            VNORM = DASUM( N, VR, 1 )
279***
280            OPST = OPST + N
281***
282            IF( VNORM.GE.GROWTO*SCALE )
283     $         GO TO 120
284*
285*           Choose new orthogonal starting vector and try again.
286*
287            TEMP = EPS3 / ( ROOTN+ONE )
288            VR( 1 ) = EPS3
289            DO 100 I = 2, N
290               VR( I ) = TEMP
291  100       CONTINUE
292            VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
293***
294            OPST = OPST + 4
295***
296  110    CONTINUE
297*
298*        Failure to find eigenvector in N iterations.
299*
300         INFO = 1
301*
302  120    CONTINUE
303*
304*        Normalize eigenvector.
305*
306         I = IDAMAX( N, VR, 1 )
307         CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
308***
309         OPST = OPST + ( 2*N+1 )
310***
311      ELSE
312*
313*        Complex eigenvalue.
314*
315         IF( NOINIT ) THEN
316*
317*           Set initial vector.
318*
319            DO 130 I = 1, N
320               VR( I ) = EPS3
321               VI( I ) = ZERO
322  130       CONTINUE
323         ELSE
324*
325*           Scale supplied initial vector.
326*
327            NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
328            REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
329            CALL DSCAL( N, REC, VR, 1 )
330            CALL DSCAL( N, REC, VI, 1 )
331***
332            OPST = OPST + ( 6*N+5 )
333***
334         END IF
335*
336         IF( RIGHTV ) THEN
337*
338*           LU decomposition with partial pivoting of B, replacing zero
339*           pivots by EPS3.
340*
341*           The imaginary part of the (i,j)-th element of U is stored in
342*           B(j+1,i).
343*
344            B( 2, 1 ) = -WI
345            DO 140 I = 2, N
346               B( I+1, 1 ) = ZERO
347  140       CONTINUE
348*
349            DO 170 I = 1, N - 1
350               ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
351               EI = H( I+1, I )
352               IF( ABSBII.LT.ABS( EI ) ) THEN
353*
354*                 Interchange rows and eliminate.
355*
356                  XR = B( I, I ) / EI
357                  XI = B( I+1, I ) / EI
358                  B( I, I ) = EI
359                  B( I+1, I ) = ZERO
360                  DO 150 J = I + 1, N
361                     TEMP = B( I+1, J )
362                     B( I+1, J ) = B( I, J ) - XR*TEMP
363                     B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
364                     B( I, J ) = TEMP
365                     B( J+1, I ) = ZERO
366  150             CONTINUE
367                  B( I+2, I ) = -WI
368                  B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
369                  B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
370***
371                  OPST = OPST + ( 4*( N-I )+6 )
372***
373               ELSE
374*
375*                 Eliminate without interchanging rows.
376*
377                  IF( ABSBII.EQ.ZERO ) THEN
378                     B( I, I ) = EPS3
379                     B( I+1, I ) = ZERO
380                     ABSBII = EPS3
381                  END IF
382                  EI = ( EI / ABSBII ) / ABSBII
383                  XR = B( I, I )*EI
384                  XI = -B( I+1, I )*EI
385                  DO 160 J = I + 1, N
386                     B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
387     $                             XI*B( J+1, I )
388                     B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
389  160             CONTINUE
390                  B( I+2, I+1 ) = B( I+2, I+1 ) - WI
391***
392                  OPST = OPST + ( 7*( N-I )+4 )
393***
394               END IF
395*
396*              Compute 1-norm of offdiagonal elements of i-th row.
397*
398               WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
399     $                     DASUM( N-I, B( I+2, I ), 1 )
400***
401               OPST = OPST + ( 2*( N-I )+4 )
402***
403  170       CONTINUE
404            IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
405     $         B( N, N ) = EPS3
406            WORK( N ) = ZERO
407*
408            I1 = N
409            I2 = 1
410            I3 = -1
411         ELSE
412*
413*           UL decomposition with partial pivoting of conjg(B),
414*           replacing zero pivots by EPS3.
415*
416*           The imaginary part of the (i,j)-th element of U is stored in
417*           B(j+1,i).
418*
419            B( N+1, N ) = WI
420            DO 180 J = 1, N - 1
421               B( N+1, J ) = ZERO
422  180       CONTINUE
423*
424            DO 210 J = N, 2, -1
425               EJ = H( J, J-1 )
426               ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
427               IF( ABSBJJ.LT.ABS( EJ ) ) THEN
428*
429*                 Interchange columns and eliminate
430*
431                  XR = B( J, J ) / EJ
432                  XI = B( J+1, J ) / EJ
433                  B( J, J ) = EJ
434                  B( J+1, J ) = ZERO
435                  DO 190 I = 1, J - 1
436                     TEMP = B( I, J-1 )
437                     B( I, J-1 ) = B( I, J ) - XR*TEMP
438                     B( J, I ) = B( J+1, I ) - XI*TEMP
439                     B( I, J ) = TEMP
440                     B( J+1, I ) = ZERO
441  190             CONTINUE
442                  B( J+1, J-1 ) = WI
443                  B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
444                  B( J, J-1 ) = B( J, J-1 ) - XR*WI
445***
446                  OPST = OPST + ( 4*( J-1 )+6 )
447***
448               ELSE
449*
450*                 Eliminate without interchange.
451*
452                  IF( ABSBJJ.EQ.ZERO ) THEN
453                     B( J, J ) = EPS3
454                     B( J+1, J ) = ZERO
455                     ABSBJJ = EPS3
456                  END IF
457                  EJ = ( EJ / ABSBJJ ) / ABSBJJ
458                  XR = B( J, J )*EJ
459                  XI = -B( J+1, J )*EJ
460                  DO 200 I = 1, J - 1
461                     B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
462     $                             XI*B( J+1, I )
463                     B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
464  200             CONTINUE
465                  B( J, J-1 ) = B( J, J-1 ) + WI
466***
467                  OPST = OPST + ( 7*( J-1 )+4 )
468***
469               END IF
470*
471*              Compute 1-norm of offdiagonal elements of j-th column.
472*
473               WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
474     $                     DASUM( J-1, B( J+1, 1 ), LDB )
475***
476               OPST = OPST + ( 2*( J-1 )+4 )
477***
478  210       CONTINUE
479            IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
480     $         B( 1, 1 ) = EPS3
481            WORK( 1 ) = ZERO
482*
483            I1 = 1
484            I2 = N
485            I3 = 1
486         END IF
487*
488         DO 270 ITS = 1, N
489            SCALE = ONE
490            VMAX = ONE
491            VCRIT = BIGNUM
492*
493*           Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
494*             or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector,
495*           overwriting (xr,xi) on (vr,vi).
496*
497            DO 250 I = I1, I2, I3
498*
499               IF( WORK( I ).GT.VCRIT ) THEN
500                  REC = ONE / VMAX
501                  CALL DSCAL( N, REC, VR, 1 )
502                  CALL DSCAL( N, REC, VI, 1 )
503                  SCALE = SCALE*REC
504                  VMAX = ONE
505                  VCRIT = BIGNUM
506               END IF
507*
508               XR = VR( I )
509               XI = VI( I )
510               IF( RIGHTV ) THEN
511                  DO 220 J = I + 1, N
512                     XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
513                     XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
514  220             CONTINUE
515               ELSE
516                  DO 230 J = 1, I - 1
517                     XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
518                     XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
519  230             CONTINUE
520               END IF
521*
522               W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
523               IF( W.GT.SMLNUM ) THEN
524                  IF( W.LT.ONE ) THEN
525                     W1 = ABS( XR ) + ABS( XI )
526                     IF( W1.GT.W*BIGNUM ) THEN
527                        REC = ONE / W1
528                        CALL DSCAL( N, REC, VR, 1 )
529                        CALL DSCAL( N, REC, VI, 1 )
530                        XR = VR( I )
531                        XI = VI( I )
532                        SCALE = SCALE*REC
533                        VMAX = VMAX*REC
534                     END IF
535                  END IF
536*
537*                 Divide by diagonal element of B.
538*
539                  CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
540     $                         VI( I ) )
541                  VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
542                  VCRIT = BIGNUM / VMAX
543***
544                  OPST = OPST + 9
545***
546               ELSE
547                  DO 240 J = 1, N
548                     VR( J ) = ZERO
549                     VI( J ) = ZERO
550  240             CONTINUE
551                  VR( I ) = ONE
552                  VI( I ) = ONE
553                  SCALE = ZERO
554                  VMAX = ONE
555                  VCRIT = BIGNUM
556               END IF
557  250       CONTINUE
558***
559*           Increment op count for loop 260, assuming no scaling
560            OPS = OPS + 4*N*( N-1 )
561***
562*
563*           Test for sufficient growth in the norm of (VR,VI).
564*
565            VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
566***
567            OPST = OPST + 2*N
568***
569            IF( VNORM.GE.GROWTO*SCALE )
570     $         GO TO 280
571*
572*           Choose a new orthogonal starting vector and try again.
573*
574            Y = EPS3 / ( ROOTN+ONE )
575            VR( 1 ) = EPS3
576            VI( 1 ) = ZERO
577*
578            DO 260 I = 2, N
579               VR( I ) = Y
580               VI( I ) = ZERO
581  260       CONTINUE
582            VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
583***
584            OPST = OPST + 4
585***
586  270    CONTINUE
587*
588*        Failure to find eigenvector in N iterations
589*
590         INFO = 1
591*
592  280    CONTINUE
593*
594*        Normalize eigenvector.
595*
596         VNORM = ZERO
597         DO 290 I = 1, N
598            VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
599  290    CONTINUE
600         CALL DSCAL( N, ONE / VNORM, VR, 1 )
601         CALL DSCAL( N, ONE / VNORM, VI, 1 )
602***
603         OPST = OPST + ( 4*N+1 )
604***
605*
606      END IF
607*
608***
609*     Compute final op count
610      OPS = OPS + OPST
611***
612      RETURN
613*
614*     End of DLAEIN
615*
616      END
617