LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slaein.f
Go to the documentation of this file.
1*> \brief \b SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLAEIN + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaein.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaein.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaein.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
22* LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
23*
24* .. Scalar Arguments ..
25* LOGICAL NOINIT, RIGHTV
26* INTEGER INFO, LDB, LDH, N
27* REAL BIGNUM, EPS3, SMLNUM, WI, WR
28* ..
29* .. Array Arguments ..
30* REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
31* $ WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> SLAEIN uses inverse iteration to find a right or left eigenvector
41*> corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
42*> matrix H.
43*> \endverbatim
44*
45* Arguments:
46* ==========
47*
48*> \param[in] RIGHTV
49*> \verbatim
50*> RIGHTV is LOGICAL
51*> = .TRUE. : compute right eigenvector;
52*> = .FALSE.: compute left eigenvector.
53*> \endverbatim
54*>
55*> \param[in] NOINIT
56*> \verbatim
57*> NOINIT is LOGICAL
58*> = .TRUE. : no initial vector supplied in (VR,VI).
59*> = .FALSE.: initial vector supplied in (VR,VI).
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*> N is INTEGER
65*> The order of the matrix H. N >= 0.
66*> \endverbatim
67*>
68*> \param[in] H
69*> \verbatim
70*> H is REAL array, dimension (LDH,N)
71*> The upper Hessenberg matrix H.
72*> \endverbatim
73*>
74*> \param[in] LDH
75*> \verbatim
76*> LDH is INTEGER
77*> The leading dimension of the array H. LDH >= max(1,N).
78*> \endverbatim
79*>
80*> \param[in] WR
81*> \verbatim
82*> WR is REAL
83*> \endverbatim
84*>
85*> \param[in] WI
86*> \verbatim
87*> WI is REAL
88*> The real and imaginary parts of the eigenvalue of H whose
89*> corresponding right or left eigenvector is to be computed.
90*> \endverbatim
91*>
92*> \param[in,out] VR
93*> \verbatim
94*> VR is REAL array, dimension (N)
95*> \endverbatim
96*>
97*> \param[in,out] VI
98*> \verbatim
99*> VI is REAL array, dimension (N)
100*> On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
101*> a real starting vector for inverse iteration using the real
102*> eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
103*> must contain the real and imaginary parts of a complex
104*> starting vector for inverse iteration using the complex
105*> eigenvalue (WR,WI); otherwise VR and VI need not be set.
106*> On exit, if WI = 0.0 (real eigenvalue), VR contains the
107*> computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
108*> VR and VI contain the real and imaginary parts of the
109*> computed complex eigenvector. The eigenvector is normalized
110*> so that the component of largest magnitude has magnitude 1;
111*> here the magnitude of a complex number (x,y) is taken to be
112*> |x| + |y|.
113*> VI is not referenced if WI = 0.0.
114*> \endverbatim
115*>
116*> \param[out] B
117*> \verbatim
118*> B is REAL array, dimension (LDB,N)
119*> \endverbatim
120*>
121*> \param[in] LDB
122*> \verbatim
123*> LDB is INTEGER
124*> The leading dimension of the array B. LDB >= N+1.
125*> \endverbatim
126*>
127*> \param[out] WORK
128*> \verbatim
129*> WORK is REAL array, dimension (N)
130*> \endverbatim
131*>
132*> \param[in] EPS3
133*> \verbatim
134*> EPS3 is REAL
135*> A small machine-dependent value which is used to perturb
136*> close eigenvalues, and to replace zero pivots.
137*> \endverbatim
138*>
139*> \param[in] SMLNUM
140*> \verbatim
141*> SMLNUM is REAL
142*> A machine-dependent value close to the underflow threshold.
143*> \endverbatim
144*>
145*> \param[in] BIGNUM
146*> \verbatim
147*> BIGNUM is REAL
148*> A machine-dependent value close to the overflow threshold.
149*> \endverbatim
150*>
151*> \param[out] INFO
152*> \verbatim
153*> INFO is INTEGER
154*> = 0: successful exit
155*> = 1: inverse iteration did not converge; VR is set to the
156*> last iterate, and so is VI if WI.ne.0.0.
157*> \endverbatim
158*
159* Authors:
160* ========
161*
162*> \author Univ. of Tennessee
163*> \author Univ. of California Berkeley
164*> \author Univ. of Colorado Denver
165*> \author NAG Ltd.
166*
167*> \ingroup laein
168*
169* =====================================================================
170 SUBROUTINE slaein( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
171 $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
172*
173* -- LAPACK auxiliary routine --
174* -- LAPACK is a software package provided by Univ. of Tennessee, --
175* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176*
177* .. Scalar Arguments ..
178 LOGICAL NOINIT, RIGHTV
179 INTEGER INFO, LDB, LDH, N
180 REAL BIGNUM, EPS3, SMLNUM, WI, WR
181* ..
182* .. Array Arguments ..
183 REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
184 $ work( * )
185* ..
186*
187* =====================================================================
188*
189* .. Parameters ..
190 REAL ZERO, ONE, TENTH
191 parameter( zero = 0.0e+0, one = 1.0e+0, tenth = 1.0e-1 )
192* ..
193* .. Local Scalars ..
194 CHARACTER NORMIN, TRANS
195 INTEGER I, I1, I2, I3, IERR, ITS, J
196 REAL ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
197 $ rec, rootn, scale, temp, vcrit, vmax, vnorm, w,
198 $ w1, x, xi, xr, y
199* ..
200* .. External Functions ..
201 INTEGER ISAMAX
202 REAL SASUM, SLAPY2, SNRM2
203 EXTERNAL isamax, sasum, slapy2, snrm2
204* ..
205* .. External Subroutines ..
206 EXTERNAL sladiv, slatrs, sscal
207* ..
208* .. Intrinsic Functions ..
209 INTRINSIC abs, max, real, sqrt
210* ..
211* .. Executable Statements ..
212*
213 info = 0
214*
215* GROWTO is the threshold used in the acceptance test for an
216* eigenvector.
217*
218 rootn = sqrt( real( n ) )
219 growto = tenth / rootn
220 nrmsml = max( one, eps3*rootn )*smlnum
221*
222* Form B = H - (WR,WI)*I (except that the subdiagonal elements and
223* the imaginary parts of the diagonal elements are not stored).
224*
225 DO 20 j = 1, n
226 DO 10 i = 1, j - 1
227 b( i, j ) = h( i, j )
228 10 CONTINUE
229 b( j, j ) = h( j, j ) - wr
230 20 CONTINUE
231*
232 IF( wi.EQ.zero ) THEN
233*
234* Real eigenvalue.
235*
236 IF( noinit ) THEN
237*
238* Set initial vector.
239*
240 DO 30 i = 1, n
241 vr( i ) = eps3
242 30 CONTINUE
243 ELSE
244*
245* Scale supplied initial vector.
246*
247 vnorm = snrm2( n, vr, 1 )
248 CALL sscal( n, ( eps3*rootn ) / max( vnorm, nrmsml ), vr,
249 $ 1 )
250 END IF
251*
252 IF( rightv ) THEN
253*
254* LU decomposition with partial pivoting of B, replacing zero
255* pivots by EPS3.
256*
257 DO 60 i = 1, n - 1
258 ei = h( i+1, i )
259 IF( abs( b( i, i ) ).LT.abs( ei ) ) THEN
260*
261* Interchange rows and eliminate.
262*
263 x = b( i, i ) / ei
264 b( i, i ) = ei
265 DO 40 j = i + 1, n
266 temp = b( i+1, j )
267 b( i+1, j ) = b( i, j ) - x*temp
268 b( i, j ) = temp
269 40 CONTINUE
270 ELSE
271*
272* Eliminate without interchange.
273*
274 IF( b( i, i ).EQ.zero )
275 $ b( i, i ) = eps3
276 x = ei / b( i, i )
277 IF( x.NE.zero ) THEN
278 DO 50 j = i + 1, n
279 b( i+1, j ) = b( i+1, j ) - x*b( i, j )
280 50 CONTINUE
281 END IF
282 END IF
283 60 CONTINUE
284 IF( b( n, n ).EQ.zero )
285 $ b( n, n ) = eps3
286*
287 trans = 'N'
288*
289 ELSE
290*
291* UL decomposition with partial pivoting of B, replacing zero
292* pivots by EPS3.
293*
294 DO 90 j = n, 2, -1
295 ej = h( j, j-1 )
296 IF( abs( b( j, j ) ).LT.abs( ej ) ) THEN
297*
298* Interchange columns and eliminate.
299*
300 x = b( j, j ) / ej
301 b( j, j ) = ej
302 DO 70 i = 1, j - 1
303 temp = b( i, j-1 )
304 b( i, j-1 ) = b( i, j ) - x*temp
305 b( i, j ) = temp
306 70 CONTINUE
307 ELSE
308*
309* Eliminate without interchange.
310*
311 IF( b( j, j ).EQ.zero )
312 $ b( j, j ) = eps3
313 x = ej / b( j, j )
314 IF( x.NE.zero ) THEN
315 DO 80 i = 1, j - 1
316 b( i, j-1 ) = b( i, j-1 ) - x*b( i, j )
317 80 CONTINUE
318 END IF
319 END IF
320 90 CONTINUE
321 IF( b( 1, 1 ).EQ.zero )
322 $ b( 1, 1 ) = eps3
323*
324 trans = 'T'
325*
326 END IF
327*
328 normin = 'N'
329 DO 110 its = 1, n
330*
331* Solve U*x = scale*v for a right eigenvector
332* or U**T*x = scale*v for a left eigenvector,
333* overwriting x on v.
334*
335 CALL slatrs( 'Upper', trans, 'Nonunit', normin, n, b, ldb,
336 $ vr, scale, work, ierr )
337 normin = 'Y'
338*
339* Test for sufficient growth in the norm of v.
340*
341 vnorm = sasum( n, vr, 1 )
342 IF( vnorm.GE.growto*scale )
343 $ GO TO 120
344*
345* Choose new orthogonal starting vector and try again.
346*
347 temp = eps3 / ( rootn+one )
348 vr( 1 ) = eps3
349 DO 100 i = 2, n
350 vr( i ) = temp
351 100 CONTINUE
352 vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
353 110 CONTINUE
354*
355* Failure to find eigenvector in N iterations.
356*
357 info = 1
358*
359 120 CONTINUE
360*
361* Normalize eigenvector.
362*
363 i = isamax( n, vr, 1 )
364 CALL sscal( n, one / abs( vr( i ) ), vr, 1 )
365 ELSE
366*
367* Complex eigenvalue.
368*
369 IF( noinit ) THEN
370*
371* Set initial vector.
372*
373 DO 130 i = 1, n
374 vr( i ) = eps3
375 vi( i ) = zero
376 130 CONTINUE
377 ELSE
378*
379* Scale supplied initial vector.
380*
381 norm = slapy2( snrm2( n, vr, 1 ), snrm2( n, vi, 1 ) )
382 rec = ( eps3*rootn ) / max( norm, nrmsml )
383 CALL sscal( n, rec, vr, 1 )
384 CALL sscal( n, rec, vi, 1 )
385 END IF
386*
387 IF( rightv ) THEN
388*
389* LU decomposition with partial pivoting of B, replacing zero
390* pivots by EPS3.
391*
392* The imaginary part of the (i,j)-th element of U is stored in
393* B(j+1,i).
394*
395 b( 2, 1 ) = -wi
396 DO 140 i = 2, n
397 b( i+1, 1 ) = zero
398 140 CONTINUE
399*
400 DO 170 i = 1, n - 1
401 absbii = slapy2( b( i, i ), b( i+1, i ) )
402 ei = h( i+1, i )
403 IF( absbii.LT.abs( ei ) ) THEN
404*
405* Interchange rows and eliminate.
406*
407 xr = b( i, i ) / ei
408 xi = b( i+1, i ) / ei
409 b( i, i ) = ei
410 b( i+1, i ) = zero
411 DO 150 j = i + 1, n
412 temp = b( i+1, j )
413 b( i+1, j ) = b( i, j ) - xr*temp
414 b( j+1, i+1 ) = b( j+1, i ) - xi*temp
415 b( i, j ) = temp
416 b( j+1, i ) = zero
417 150 CONTINUE
418 b( i+2, i ) = -wi
419 b( i+1, i+1 ) = b( i+1, i+1 ) - xi*wi
420 b( i+2, i+1 ) = b( i+2, i+1 ) + xr*wi
421 ELSE
422*
423* Eliminate without interchanging rows.
424*
425 IF( absbii.EQ.zero ) THEN
426 b( i, i ) = eps3
427 b( i+1, i ) = zero
428 absbii = eps3
429 END IF
430 ei = ( ei / absbii ) / absbii
431 xr = b( i, i )*ei
432 xi = -b( i+1, i )*ei
433 DO 160 j = i + 1, n
434 b( i+1, j ) = b( i+1, j ) - xr*b( i, j ) +
435 $ xi*b( j+1, i )
436 b( j+1, i+1 ) = -xr*b( j+1, i ) - xi*b( i, j )
437 160 CONTINUE
438 b( i+2, i+1 ) = b( i+2, i+1 ) - wi
439 END IF
440*
441* Compute 1-norm of offdiagonal elements of i-th row.
442*
443 work( i ) = sasum( n-i, b( i, i+1 ), ldb ) +
444 $ sasum( n-i, b( i+2, i ), 1 )
445 170 CONTINUE
446 IF( b( n, n ).EQ.zero .AND. b( n+1, n ).EQ.zero )
447 $ b( n, n ) = eps3
448 work( n ) = zero
449*
450 i1 = n
451 i2 = 1
452 i3 = -1
453 ELSE
454*
455* UL decomposition with partial pivoting of conjg(B),
456* replacing zero pivots by EPS3.
457*
458* The imaginary part of the (i,j)-th element of U is stored in
459* B(j+1,i).
460*
461 b( n+1, n ) = wi
462 DO 180 j = 1, n - 1
463 b( n+1, j ) = zero
464 180 CONTINUE
465*
466 DO 210 j = n, 2, -1
467 ej = h( j, j-1 )
468 absbjj = slapy2( b( j, j ), b( j+1, j ) )
469 IF( absbjj.LT.abs( ej ) ) THEN
470*
471* Interchange columns and eliminate
472*
473 xr = b( j, j ) / ej
474 xi = b( j+1, j ) / ej
475 b( j, j ) = ej
476 b( j+1, j ) = zero
477 DO 190 i = 1, j - 1
478 temp = b( i, j-1 )
479 b( i, j-1 ) = b( i, j ) - xr*temp
480 b( j, i ) = b( j+1, i ) - xi*temp
481 b( i, j ) = temp
482 b( j+1, i ) = zero
483 190 CONTINUE
484 b( j+1, j-1 ) = wi
485 b( j-1, j-1 ) = b( j-1, j-1 ) + xi*wi
486 b( j, j-1 ) = b( j, j-1 ) - xr*wi
487 ELSE
488*
489* Eliminate without interchange.
490*
491 IF( absbjj.EQ.zero ) THEN
492 b( j, j ) = eps3
493 b( j+1, j ) = zero
494 absbjj = eps3
495 END IF
496 ej = ( ej / absbjj ) / absbjj
497 xr = b( j, j )*ej
498 xi = -b( j+1, j )*ej
499 DO 200 i = 1, j - 1
500 b( i, j-1 ) = b( i, j-1 ) - xr*b( i, j ) +
501 $ xi*b( j+1, i )
502 b( j, i ) = -xr*b( j+1, i ) - xi*b( i, j )
503 200 CONTINUE
504 b( j, j-1 ) = b( j, j-1 ) + wi
505 END IF
506*
507* Compute 1-norm of offdiagonal elements of j-th column.
508*
509 work( j ) = sasum( j-1, b( 1, j ), 1 ) +
510 $ sasum( j-1, b( j+1, 1 ), ldb )
511 210 CONTINUE
512 IF( b( 1, 1 ).EQ.zero .AND. b( 2, 1 ).EQ.zero )
513 $ b( 1, 1 ) = eps3
514 work( 1 ) = zero
515*
516 i1 = 1
517 i2 = n
518 i3 = 1
519 END IF
520*
521 DO 270 its = 1, n
522 scale = one
523 vmax = one
524 vcrit = bignum
525*
526* Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
527* or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
528* overwriting (xr,xi) on (vr,vi).
529*
530 DO 250 i = i1, i2, i3
531*
532 IF( work( i ).GT.vcrit ) THEN
533 rec = one / vmax
534 CALL sscal( n, rec, vr, 1 )
535 CALL sscal( n, rec, vi, 1 )
536 scale = scale*rec
537 vmax = one
538 vcrit = bignum
539 END IF
540*
541 xr = vr( i )
542 xi = vi( i )
543 IF( rightv ) THEN
544 DO 220 j = i + 1, n
545 xr = xr - b( i, j )*vr( j ) + b( j+1, i )*vi( j )
546 xi = xi - b( i, j )*vi( j ) - b( j+1, i )*vr( j )
547 220 CONTINUE
548 ELSE
549 DO 230 j = 1, i - 1
550 xr = xr - b( j, i )*vr( j ) + b( i+1, j )*vi( j )
551 xi = xi - b( j, i )*vi( j ) - b( i+1, j )*vr( j )
552 230 CONTINUE
553 END IF
554*
555 w = abs( b( i, i ) ) + abs( b( i+1, i ) )
556 IF( w.GT.smlnum ) THEN
557 IF( w.LT.one ) THEN
558 w1 = abs( xr ) + abs( xi )
559 IF( w1.GT.w*bignum ) THEN
560 rec = one / w1
561 CALL sscal( n, rec, vr, 1 )
562 CALL sscal( n, rec, vi, 1 )
563 xr = vr( i )
564 xi = vi( i )
565 scale = scale*rec
566 vmax = vmax*rec
567 END IF
568 END IF
569*
570* Divide by diagonal element of B.
571*
572 CALL sladiv( xr, xi, b( i, i ), b( i+1, i ), vr( i ),
573 $ vi( i ) )
574 vmax = max( abs( vr( i ) )+abs( vi( i ) ), vmax )
575 vcrit = bignum / vmax
576 ELSE
577 DO 240 j = 1, n
578 vr( j ) = zero
579 vi( j ) = zero
580 240 CONTINUE
581 vr( i ) = one
582 vi( i ) = one
583 scale = zero
584 vmax = one
585 vcrit = bignum
586 END IF
587 250 CONTINUE
588*
589* Test for sufficient growth in the norm of (VR,VI).
590*
591 vnorm = sasum( n, vr, 1 ) + sasum( n, vi, 1 )
592 IF( vnorm.GE.growto*scale )
593 $ GO TO 280
594*
595* Choose a new orthogonal starting vector and try again.
596*
597 y = eps3 / ( rootn+one )
598 vr( 1 ) = eps3
599 vi( 1 ) = zero
600*
601 DO 260 i = 2, n
602 vr( i ) = y
603 vi( i ) = zero
604 260 CONTINUE
605 vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
606 270 CONTINUE
607*
608* Failure to find eigenvector in N iterations
609*
610 info = 1
611*
612 280 CONTINUE
613*
614* Normalize eigenvector.
615*
616 vnorm = zero
617 DO 290 i = 1, n
618 vnorm = max( vnorm, abs( vr( i ) )+abs( vi( i ) ) )
619 290 CONTINUE
620 CALL sscal( n, one / vnorm, vr, 1 )
621 CALL sscal( n, one / vnorm, vi, 1 )
622*
623 END IF
624*
625 RETURN
626*
627* End of SLAEIN
628*
629 END
subroutine sladiv(a, b, c, d, p, q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition sladiv.f:91
subroutine slaein(rightv, noinit, n, h, ldh, wr, wi, vr, vi, b, ldb, work, eps3, smlnum, bignum, info)
SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition slaein.f:172
subroutine slatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
SLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition slatrs.f:238
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79