LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ claein()

subroutine claein ( logical  rightv,
logical  noinit,
integer  n,
complex, dimension( ldh, * )  h,
integer  ldh,
complex  w,
complex, dimension( * )  v,
complex, dimension( ldb, * )  b,
integer  ldb,
real, dimension( * )  rwork,
real  eps3,
real  smlnum,
integer  info 
)

CLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.

Download CLAEIN + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CLAEIN uses inverse iteration to find a right or left eigenvector
 corresponding to the eigenvalue W of a complex upper Hessenberg
 matrix H.
Parameters
[in]RIGHTV
          RIGHTV is LOGICAL
          = .TRUE. : compute right eigenvector;
          = .FALSE.: compute left eigenvector.
[in]NOINIT
          NOINIT is LOGICAL
          = .TRUE. : no initial vector supplied in V
          = .FALSE.: initial vector supplied in V.
[in]N
          N is INTEGER
          The order of the matrix H.  N >= 0.
[in]H
          H is COMPLEX array, dimension (LDH,N)
          The upper Hessenberg matrix H.
[in]LDH
          LDH is INTEGER
          The leading dimension of the array H.  LDH >= max(1,N).
[in]W
          W is COMPLEX
          The eigenvalue of H whose corresponding right or left
          eigenvector is to be computed.
[in,out]V
          V is COMPLEX array, dimension (N)
          On entry, if NOINIT = .FALSE., V must contain a starting
          vector for inverse iteration; otherwise V need not be set.
          On exit, V contains the computed eigenvector, normalized so
          that the component of largest magnitude has magnitude 1; here
          the magnitude of a complex number (x,y) is taken to be
          |x| + |y|.
[out]B
          B is COMPLEX array, dimension (LDB,N)
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[out]RWORK
          RWORK is REAL array, dimension (N)
[in]EPS3
          EPS3 is REAL
          A small machine-dependent value which is used to perturb
          close eigenvalues, and to replace zero pivots.
[in]SMLNUM
          SMLNUM is REAL
          A machine-dependent value close to the underflow threshold.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          = 1:  inverse iteration did not converge; V is set to the
                last iterate.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 147 of file claein.f.

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*
scasum
real function scasum(n, cx, incx)
SCASUM
Definition scasum.f:72
icamax
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
cladiv
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition cladiv.f:64
clatrs
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
scnrm2
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition scnrm2.f90:90
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: