LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
claein.f
Go to the documentation of this file.
1*> \brief \b CLAEIN 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 CLAEIN + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claein.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claein.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claein.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
22* EPS3, SMLNUM, INFO )
23*
24* .. Scalar Arguments ..
25* LOGICAL NOINIT, RIGHTV
26* INTEGER INFO, LDB, LDH, N
27* REAL EPS3, SMLNUM
28* COMPLEX W
29* ..
30* .. Array Arguments ..
31* REAL RWORK( * )
32* COMPLEX B( LDB, * ), H( LDH, * ), V( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> CLAEIN uses inverse iteration to find a right or left eigenvector
42*> corresponding to the eigenvalue W of a complex upper Hessenberg
43*> matrix H.
44*> \endverbatim
45*
46* Arguments:
47* ==========
48*
49*> \param[in] RIGHTV
50*> \verbatim
51*> RIGHTV is LOGICAL
52*> = .TRUE. : compute right eigenvector;
53*> = .FALSE.: compute left eigenvector.
54*> \endverbatim
55*>
56*> \param[in] NOINIT
57*> \verbatim
58*> NOINIT is LOGICAL
59*> = .TRUE. : no initial vector supplied in V
60*> = .FALSE.: initial vector supplied in V.
61*> \endverbatim
62*>
63*> \param[in] N
64*> \verbatim
65*> N is INTEGER
66*> The order of the matrix H. N >= 0.
67*> \endverbatim
68*>
69*> \param[in] H
70*> \verbatim
71*> H is COMPLEX array, dimension (LDH,N)
72*> The upper Hessenberg matrix H.
73*> \endverbatim
74*>
75*> \param[in] LDH
76*> \verbatim
77*> LDH is INTEGER
78*> The leading dimension of the array H. LDH >= max(1,N).
79*> \endverbatim
80*>
81*> \param[in] W
82*> \verbatim
83*> W is COMPLEX
84*> The eigenvalue of H whose corresponding right or left
85*> eigenvector is to be computed.
86*> \endverbatim
87*>
88*> \param[in,out] V
89*> \verbatim
90*> V is COMPLEX array, dimension (N)
91*> On entry, if NOINIT = .FALSE., V must contain a starting
92*> vector for inverse iteration; otherwise V need not be set.
93*> On exit, V contains the computed eigenvector, normalized so
94*> that the component of largest magnitude has magnitude 1; here
95*> the magnitude of a complex number (x,y) is taken to be
96*> |x| + |y|.
97*> \endverbatim
98*>
99*> \param[out] B
100*> \verbatim
101*> B is COMPLEX array, dimension (LDB,N)
102*> \endverbatim
103*>
104*> \param[in] LDB
105*> \verbatim
106*> LDB is INTEGER
107*> The leading dimension of the array B. LDB >= max(1,N).
108*> \endverbatim
109*>
110*> \param[out] RWORK
111*> \verbatim
112*> RWORK is REAL array, dimension (N)
113*> \endverbatim
114*>
115*> \param[in] EPS3
116*> \verbatim
117*> EPS3 is REAL
118*> A small machine-dependent value which is used to perturb
119*> close eigenvalues, and to replace zero pivots.
120*> \endverbatim
121*>
122*> \param[in] SMLNUM
123*> \verbatim
124*> SMLNUM is REAL
125*> A machine-dependent value close to the underflow threshold.
126*> \endverbatim
127*>
128*> \param[out] INFO
129*> \verbatim
130*> INFO is INTEGER
131*> = 0: successful exit
132*> = 1: inverse iteration did not converge; V is set to the
133*> last iterate.
134*> \endverbatim
135*
136* Authors:
137* ========
138*
139*> \author Univ. of Tennessee
140*> \author Univ. of California Berkeley
141*> \author Univ. of Colorado Denver
142*> \author NAG Ltd.
143*
144*> \ingroup laein
145*
146* =====================================================================
147 SUBROUTINE claein( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
148 $ EPS3, SMLNUM, INFO )
149*
150* -- LAPACK auxiliary routine --
151* -- LAPACK is a software package provided by Univ. of Tennessee, --
152* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
153*
154* .. Scalar Arguments ..
155 LOGICAL NOINIT, RIGHTV
156 INTEGER INFO, LDB, LDH, N
157 REAL EPS3, SMLNUM
158 COMPLEX W
159* ..
160* .. Array Arguments ..
161 REAL RWORK( * )
162 COMPLEX B( LDB, * ), H( LDH, * ), V( * )
163* ..
164*
165* =====================================================================
166*
167* .. Parameters ..
168 REAL ONE, TENTH
169 parameter( one = 1.0e+0, tenth = 1.0e-1 )
170 COMPLEX ZERO
171 parameter( zero = ( 0.0e+0, 0.0e+0 ) )
172* ..
173* .. Local Scalars ..
174 CHARACTER NORMIN, TRANS
175 INTEGER I, IERR, ITS, J
176 REAL GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
177 COMPLEX CDUM, EI, EJ, TEMP, X
178* ..
179* .. External Functions ..
180 INTEGER ICAMAX
181 REAL SCASUM, SCNRM2
182 COMPLEX CLADIV
183 EXTERNAL icamax, scasum, scnrm2, cladiv
184* ..
185* .. External Subroutines ..
186 EXTERNAL clatrs, csscal
187* ..
188* .. Intrinsic Functions ..
189 INTRINSIC abs, aimag, max, real, sqrt
190* ..
191* .. Statement Functions ..
192 REAL CABS1
193* ..
194* .. Statement Function definitions ..
195 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
196* ..
197* .. Executable Statements ..
198*
199 info = 0
200*
201* GROWTO is the threshold used in the acceptance test for an
202* eigenvector.
203*
204 rootn = sqrt( real( n ) )
205 growto = tenth / rootn
206 nrmsml = max( one, eps3*rootn )*smlnum
207*
208* Form B = H - W*I (except that the subdiagonal elements are not
209* stored).
210*
211 DO 20 j = 1, n
212 DO 10 i = 1, j - 1
213 b( i, j ) = h( i, j )
214 10 CONTINUE
215 b( j, j ) = h( j, j ) - w
216 20 CONTINUE
217*
218 IF( noinit ) THEN
219*
220* Initialize V.
221*
222 DO 30 i = 1, n
223 v( i ) = eps3
224 30 CONTINUE
225 ELSE
226*
227* Scale supplied initial vector.
228*
229 vnorm = scnrm2( n, v, 1 )
230 CALL csscal( n, ( eps3*rootn ) / max( vnorm, nrmsml ), v, 1 )
231 END IF
232*
233 IF( rightv ) THEN
234*
235* LU decomposition with partial pivoting of B, replacing zero
236* pivots by EPS3.
237*
238 DO 60 i = 1, n - 1
239 ei = h( i+1, i )
240 IF( cabs1( b( i, i ) ).LT.cabs1( ei ) ) THEN
241*
242* Interchange rows and eliminate.
243*
244 x = cladiv( b( i, i ), ei )
245 b( i, i ) = ei
246 DO 40 j = i + 1, n
247 temp = b( i+1, j )
248 b( i+1, j ) = b( i, j ) - x*temp
249 b( i, j ) = temp
250 40 CONTINUE
251 ELSE
252*
253* Eliminate without interchange.
254*
255 IF( b( i, i ).EQ.zero )
256 $ b( i, i ) = eps3
257 x = cladiv( ei, b( i, i ) )
258 IF( x.NE.zero ) THEN
259 DO 50 j = i + 1, n
260 b( i+1, j ) = b( i+1, j ) - x*b( i, j )
261 50 CONTINUE
262 END IF
263 END IF
264 60 CONTINUE
265 IF( b( n, n ).EQ.zero )
266 $ b( n, n ) = eps3
267*
268 trans = 'N'
269*
270 ELSE
271*
272* UL decomposition with partial pivoting of B, replacing zero
273* pivots by EPS3.
274*
275 DO 90 j = n, 2, -1
276 ej = h( j, j-1 )
277 IF( cabs1( b( j, j ) ).LT.cabs1( ej ) ) THEN
278*
279* Interchange columns and eliminate.
280*
281 x = cladiv( b( j, j ), ej )
282 b( j, j ) = ej
283 DO 70 i = 1, j - 1
284 temp = b( i, j-1 )
285 b( i, j-1 ) = b( i, j ) - x*temp
286 b( i, j ) = temp
287 70 CONTINUE
288 ELSE
289*
290* Eliminate without interchange.
291*
292 IF( b( j, j ).EQ.zero )
293 $ b( j, j ) = eps3
294 x = cladiv( ej, b( j, j ) )
295 IF( x.NE.zero ) THEN
296 DO 80 i = 1, j - 1
297 b( i, j-1 ) = b( i, j-1 ) - x*b( i, j )
298 80 CONTINUE
299 END IF
300 END IF
301 90 CONTINUE
302 IF( b( 1, 1 ).EQ.zero )
303 $ b( 1, 1 ) = eps3
304*
305 trans = 'C'
306*
307 END IF
308*
309 normin = 'N'
310 DO 110 its = 1, n
311*
312* Solve U*x = scale*v for a right eigenvector
313* or U**H *x = scale*v for a left eigenvector,
314* overwriting x on v.
315*
316 CALL clatrs( 'Upper', trans, 'Nonunit', normin, n, b, ldb, v,
317 $ scale, rwork, ierr )
318 normin = 'Y'
319*
320* Test for sufficient growth in the norm of v.
321*
322 vnorm = scasum( n, v, 1 )
323 IF( vnorm.GE.growto*scale )
324 $ GO TO 120
325*
326* Choose new orthogonal starting vector and try again.
327*
328 rtemp = eps3 / ( rootn+one )
329 v( 1 ) = eps3
330 DO 100 i = 2, n
331 v( i ) = rtemp
332 100 CONTINUE
333 v( n-its+1 ) = v( n-its+1 ) - eps3*rootn
334 110 CONTINUE
335*
336* Failure to find eigenvector in N iterations.
337*
338 info = 1
339*
340 120 CONTINUE
341*
342* Normalize eigenvector.
343*
344 i = icamax( n, v, 1 )
345 CALL csscal( n, one / cabs1( v( i ) ), v, 1 )
346*
347 RETURN
348*
349* End of CLAEIN
350*
351 END
subroutine claein(rightv, noinit, n, h, ldh, w, v, b, ldb, rwork, eps3, smlnum, info)
CLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition claein.f:149
subroutine clatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
CLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition clatrs.f:239
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78