1      SUBROUTINE SLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
2     $                   LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
3*
4*  -- LAPACK auxiliary routine (version 3.0) --
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      REAL               BIGNUM, EPS3, SMLNUM, WI, WR
13*     ..
14*     .. Array Arguments ..
15      REAL               B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
16     $                   WORK( * )
17*     ..
18*
19*  Purpose
20*  =======
21*
22*  SLAEIN uses inverse iteration to find a right or left eigenvector
23*  corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
24*  matrix H.
25*
26*  Arguments
27*  =========
28*
29*  RIGHTV   (input) LOGICAL
30*          = .TRUE. : compute right eigenvector;
31*          = .FALSE.: compute left eigenvector.
32*
33*  NOINIT   (input) LOGICAL
34*          = .TRUE. : no initial vector supplied in (VR,VI).
35*          = .FALSE.: initial vector supplied in (VR,VI).
36*
37*  N       (input) INTEGER
38*          The order of the matrix H.  N >= 0.
39*
40*  H       (input) REAL array, dimension (LDH,N)
41*          The upper Hessenberg matrix H.
42*
43*  LDH     (input) INTEGER
44*          The leading dimension of the array H.  LDH >= max(1,N).
45*
46*  WR      (input) REAL
47*  WI      (input) REAL
48*          The real and imaginary parts of the eigenvalue of H whose
49*          corresponding right or left eigenvector is to be computed.
50*
51*  VR      (input/output) REAL array, dimension (N)
52*  VI      (input/output) REAL array, dimension (N)
53*          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
54*          a real starting vector for inverse iteration using the real
55*          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
56*          must contain the real and imaginary parts of a complex
57*          starting vector for inverse iteration using the complex
58*          eigenvalue (WR,WI); otherwise VR and VI need not be set.
59*          On exit, if WI = 0.0 (real eigenvalue), VR contains the
60*          computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
61*          VR and VI contain the real and imaginary parts of the
62*          computed complex eigenvector. The eigenvector is normalized
63*          so that the component of largest magnitude has magnitude 1;
64*          here the magnitude of a complex number (x,y) is taken to be
65*          |x| + |y|.
66*          VI is not referenced if WI = 0.0.
67*
68*  B       (workspace) REAL array, dimension (LDB,N)
69*
70*  LDB     (input) INTEGER
71*          The leading dimension of the array B.  LDB >= N+1.
72*
73*  WORK   (workspace) REAL array, dimension (N)
74*
75*  EPS3    (input) REAL
76*          A small machine-dependent value which is used to perturb
77*          close eigenvalues, and to replace zero pivots.
78*
79*  SMLNUM  (input) REAL
80*          A machine-dependent value close to the underflow threshold.
81*
82*  BIGNUM  (input) REAL
83*          A machine-dependent value close to the overflow threshold.
84*
85*  INFO    (output) INTEGER
86*          = 0:  successful exit
87*          = 1:  inverse iteration did not converge; VR is set to the
88*                last iterate, and so is VI if WI.ne.0.0.
89*
90*  =====================================================================
91*
92*     .. Parameters ..
93      REAL               ZERO, ONE, TENTH
94      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TENTH = 1.0E-1 )
95*     ..
96*     .. Local Scalars ..
97      CHARACTER          NORMIN, TRANS
98      INTEGER            I, I1, I2, I3, IERR, ITS, J
99      REAL               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            ISAMAX
105      REAL               SASUM, SLAPY2, SNRM2
106      EXTERNAL           ISAMAX, SASUM, SLAPY2, SNRM2
107*     ..
108*     .. External Subroutines ..
109      EXTERNAL           SLADIV, SLATRS, SSCAL
110*     ..
111*     .. Intrinsic Functions ..
112      INTRINSIC          ABS, MAX, REAL, 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( REAL( 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 = SNRM2( N, VR, 1 )
151            CALL SSCAL( 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 SLATRS( '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 = SASUM( 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 = ISAMAX( N, VR, 1 )
267         CALL SSCAL( 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 = SLAPY2( SNRM2( N, VR, 1 ), SNRM2( N, VI, 1 ) )
285            REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
286            CALL SSCAL( N, REC, VR, 1 )
287            CALL SSCAL( 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 = SLAPY2( 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 ) = SASUM( N-I, B( I, I+1 ), LDB ) +
347     $                     SASUM( 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 = SLAPY2( 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 ) = SASUM( J-1, B( 1, J ), 1 ) +
413     $                     SASUM( 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 SSCAL( N, REC, VR, 1 )
438                  CALL SSCAL( 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 SSCAL( N, REC, VR, 1 )
465                        CALL SSCAL( 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 SLADIV( 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 = SASUM( N, VR, 1 ) + SASUM( 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 SSCAL( N, ONE / VNORM, VR, 1 )
524         CALL SSCAL( N, ONE / VNORM, VI, 1 )
525*
526      END IF
527*
528      RETURN
529*
530*     End of SLAEIN
531*
532      END
533