LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zstein.f
Go to the documentation of this file.
1*> \brief \b ZSTEIN
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstein.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstein.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstein.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
22* IWORK, IFAIL, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDZ, M, N
26* ..
27* .. Array Arguments ..
28* INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
29* \$ IWORK( * )
30* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
31* COMPLEX*16 Z( LDZ, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> ZSTEIN computes the eigenvectors of a real symmetric tridiagonal
41*> matrix T corresponding to specified eigenvalues, using inverse
42*> iteration.
43*>
44*> The maximum number of iterations allowed for each eigenvector is
45*> specified by an internal parameter MAXITS (currently set to 5).
46*>
47*> Although the eigenvectors are real, they are stored in a complex
48*> array, which may be passed to ZUNMTR or ZUPMTR for back
49*> transformation to the eigenvectors of a complex Hermitian matrix
50*> which was reduced to tridiagonal form.
51*>
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] N
58*> \verbatim
59*> N is INTEGER
60*> The order of the matrix. N >= 0.
61*> \endverbatim
62*>
63*> \param[in] D
64*> \verbatim
65*> D is DOUBLE PRECISION array, dimension (N)
66*> The n diagonal elements of the tridiagonal matrix T.
67*> \endverbatim
68*>
69*> \param[in] E
70*> \verbatim
71*> E is DOUBLE PRECISION array, dimension (N-1)
72*> The (n-1) subdiagonal elements of the tridiagonal matrix
73*> T, stored in elements 1 to N-1.
74*> \endverbatim
75*>
76*> \param[in] M
77*> \verbatim
78*> M is INTEGER
79*> The number of eigenvectors to be found. 0 <= M <= N.
80*> \endverbatim
81*>
82*> \param[in] W
83*> \verbatim
84*> W is DOUBLE PRECISION array, dimension (N)
85*> The first M elements of W contain the eigenvalues for
86*> which eigenvectors are to be computed. The eigenvalues
87*> should be grouped by split-off block and ordered from
88*> smallest to largest within the block. ( The output array
89*> W from DSTEBZ with ORDER = 'B' is expected here. )
90*> \endverbatim
91*>
92*> \param[in] IBLOCK
93*> \verbatim
94*> IBLOCK is INTEGER array, dimension (N)
95*> The submatrix indices associated with the corresponding
96*> eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
97*> the first submatrix from the top, =2 if W(i) belongs to
98*> the second submatrix, etc. ( The output array IBLOCK
99*> from DSTEBZ is expected here. )
100*> \endverbatim
101*>
102*> \param[in] ISPLIT
103*> \verbatim
104*> ISPLIT is INTEGER array, dimension (N)
105*> The splitting points, at which T breaks up into submatrices.
106*> The first submatrix consists of rows/columns 1 to
107*> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
108*> through ISPLIT( 2 ), etc.
109*> ( The output array ISPLIT from DSTEBZ is expected here. )
110*> \endverbatim
111*>
112*> \param[out] Z
113*> \verbatim
114*> Z is COMPLEX*16 array, dimension (LDZ, M)
115*> The computed eigenvectors. The eigenvector associated
116*> with the eigenvalue W(i) is stored in the i-th column of
117*> Z. Any vector which fails to converge is set to its current
118*> iterate after MAXITS iterations.
119*> The imaginary parts of the eigenvectors are set to zero.
120*> \endverbatim
121*>
122*> \param[in] LDZ
123*> \verbatim
124*> LDZ is INTEGER
125*> The leading dimension of the array Z. LDZ >= max(1,N).
126*> \endverbatim
127*>
128*> \param[out] WORK
129*> \verbatim
130*> WORK is DOUBLE PRECISION array, dimension (5*N)
131*> \endverbatim
132*>
133*> \param[out] IWORK
134*> \verbatim
135*> IWORK is INTEGER array, dimension (N)
136*> \endverbatim
137*>
138*> \param[out] IFAIL
139*> \verbatim
140*> IFAIL is INTEGER array, dimension (M)
141*> On normal exit, all elements of IFAIL are zero.
142*> If one or more eigenvectors fail to converge after
143*> MAXITS iterations, then their indices are stored in
144*> array IFAIL.
145*> \endverbatim
146*>
147*> \param[out] INFO
148*> \verbatim
149*> INFO is INTEGER
150*> = 0: successful exit
151*> < 0: if INFO = -i, the i-th argument had an illegal value
152*> > 0: if INFO = i, then i eigenvectors failed to converge
153*> in MAXITS iterations. Their indices are stored in
154*> array IFAIL.
155*> \endverbatim
156*
157*> \par Internal Parameters:
158* =========================
159*>
160*> \verbatim
161*> MAXITS INTEGER, default = 5
162*> The maximum number of iterations performed.
163*>
164*> EXTRA INTEGER, default = 2
165*> The number of iterations performed after norm growth
166*> criterion is satisfied, should be at least 1.
167*> \endverbatim
168*
169* Authors:
170* ========
171*
172*> \author Univ. of Tennessee
173*> \author Univ. of California Berkeley
174*> \author Univ. of Colorado Denver
175*> \author NAG Ltd.
176*
177*> \ingroup stein
178*
179* =====================================================================
180 SUBROUTINE zstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
181 \$ IWORK, IFAIL, INFO )
182*
183* -- LAPACK computational routine --
184* -- LAPACK is a software package provided by Univ. of Tennessee, --
185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187* .. Scalar Arguments ..
188 INTEGER INFO, LDZ, M, N
189* ..
190* .. Array Arguments ..
191 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
192 \$ iwork( * )
193 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
194 COMPLEX*16 Z( LDZ, * )
195* ..
196*
197* =====================================================================
198*
199* .. Parameters ..
200 COMPLEX*16 CZERO, CONE
201 parameter( czero = ( 0.0d+0, 0.0d+0 ),
202 \$ cone = ( 1.0d+0, 0.0d+0 ) )
203 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
204 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
205 \$ odm3 = 1.0d-3, odm1 = 1.0d-1 )
206 INTEGER MAXITS, EXTRA
207 parameter( maxits = 5, extra = 2 )
208* ..
209* .. Local Scalars ..
210 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
211 \$ indrv2, indrv3, indrv4, indrv5, its, j, j1,
212 \$ jblk, jmax, jr, nblk, nrmchk
213 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
214 \$ scl, sep, tol, xj, xjm, ztr
215* ..
216* .. Local Arrays ..
217 INTEGER ISEED( 4 )
218* ..
219* .. External Functions ..
220 INTEGER IDAMAX
221 DOUBLE PRECISION DLAMCH, DNRM2
222 EXTERNAL idamax, dlamch, dnrm2
223* ..
224* .. External Subroutines ..
225 EXTERNAL dcopy, dlagtf, dlagts, dlarnv, dscal, xerbla
226* ..
227* .. Intrinsic Functions ..
228 INTRINSIC abs, dble, dcmplx, max, sqrt
229* ..
230* .. Executable Statements ..
231*
232* Test the input parameters.
233*
234 info = 0
235 DO 10 i = 1, m
236 ifail( i ) = 0
237 10 CONTINUE
238*
239 IF( n.LT.0 ) THEN
240 info = -1
241 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
242 info = -4
243 ELSE IF( ldz.LT.max( 1, n ) ) THEN
244 info = -9
245 ELSE
246 DO 20 j = 2, m
247 IF( iblock( j ).LT.iblock( j-1 ) ) THEN
248 info = -6
249 GO TO 30
250 END IF
251 IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
252 \$ THEN
253 info = -5
254 GO TO 30
255 END IF
256 20 CONTINUE
257 30 CONTINUE
258 END IF
259*
260 IF( info.NE.0 ) THEN
261 CALL xerbla( 'ZSTEIN', -info )
262 RETURN
263 END IF
264*
265* Quick return if possible
266*
267 IF( n.EQ.0 .OR. m.EQ.0 ) THEN
268 RETURN
269 ELSE IF( n.EQ.1 ) THEN
270 z( 1, 1 ) = cone
271 RETURN
272 END IF
273*
274* Get machine constants.
275*
276 eps = dlamch( 'Precision' )
277*
278* Initialize seed for random number generator DLARNV.
279*
280 DO 40 i = 1, 4
281 iseed( i ) = 1
282 40 CONTINUE
283*
284* Initialize pointers.
285*
286 indrv1 = 0
287 indrv2 = indrv1 + n
288 indrv3 = indrv2 + n
289 indrv4 = indrv3 + n
290 indrv5 = indrv4 + n
291*
292* Compute eigenvectors of matrix blocks.
293*
294 j1 = 1
295 DO 180 nblk = 1, iblock( m )
296*
297* Find starting and ending indices of block nblk.
298*
299 IF( nblk.EQ.1 ) THEN
300 b1 = 1
301 ELSE
302 b1 = isplit( nblk-1 ) + 1
303 END IF
304 bn = isplit( nblk )
305 blksiz = bn - b1 + 1
306 IF( blksiz.EQ.1 )
307 \$ GO TO 60
308 gpind = j1
309*
310* Compute reorthogonalization criterion and stopping criterion.
311*
312 onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
313 onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
314 DO 50 i = b1 + 1, bn - 1
315 onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
316 \$ abs( e( i ) ) )
317 50 CONTINUE
318 ortol = odm3*onenrm
319*
320 dtpcrt = sqrt( odm1 / blksiz )
321*
322* Loop through eigenvalues of block nblk.
323*
324 60 CONTINUE
325 jblk = 0
326 DO 170 j = j1, m
327 IF( iblock( j ).NE.nblk ) THEN
328 j1 = j
329 GO TO 180
330 END IF
331 jblk = jblk + 1
332 xj = w( j )
333*
334* Skip all the work if the block size is one.
335*
336 IF( blksiz.EQ.1 ) THEN
337 work( indrv1+1 ) = one
338 GO TO 140
339 END IF
340*
341* If eigenvalues j and j-1 are too close, add a relatively
342* small perturbation.
343*
344 IF( jblk.GT.1 ) THEN
345 eps1 = abs( eps*xj )
346 pertol = ten*eps1
347 sep = xj - xjm
348 IF( sep.LT.pertol )
349 \$ xj = xjm + pertol
350 END IF
351*
352 its = 0
353 nrmchk = 0
354*
355* Get random starting vector.
356*
357 CALL dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
358*
359* Copy the matrix T so it won't be destroyed in factorization.
360*
361 CALL dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
362 CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
363 CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
364*
365* Compute LU factors with partial pivoting ( PT = LU )
366*
367 tol = zero
368 CALL dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
369 \$ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
370 \$ iinfo )
371*
372* Update iteration count.
373*
374 70 CONTINUE
375 its = its + 1
376 IF( its.GT.maxits )
377 \$ GO TO 120
378*
379* Normalize and scale the righthand side vector Pb.
380*
381 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
382 scl = blksiz*onenrm*max( eps,
383 \$ abs( work( indrv4+blksiz ) ) ) /
384 \$ abs( work( indrv1+jmax ) )
385 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
386*
387* Solve the system LU = Pb.
388*
389 CALL dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
390 \$ work( indrv3+1 ), work( indrv5+1 ), iwork,
391 \$ work( indrv1+1 ), tol, iinfo )
392*
393* Reorthogonalize by modified Gram-Schmidt if eigenvalues are
394* close enough.
395*
396 IF( jblk.EQ.1 )
397 \$ GO TO 110
398 IF( abs( xj-xjm ).GT.ortol )
399 \$ gpind = j
400 IF( gpind.NE.j ) THEN
401 DO 100 i = gpind, j - 1
402 ztr = zero
403 DO 80 jr = 1, blksiz
404 ztr = ztr + work( indrv1+jr )*
405 \$ dble( z( b1-1+jr, i ) )
406 80 CONTINUE
407 DO 90 jr = 1, blksiz
408 work( indrv1+jr ) = work( indrv1+jr ) -
409 \$ ztr*dble( z( b1-1+jr, i ) )
410 90 CONTINUE
411 100 CONTINUE
412 END IF
413*
414* Check the infinity norm of the iterate.
415*
416 110 CONTINUE
417 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
418 nrm = abs( work( indrv1+jmax ) )
419*
420* Continue for additional iterations after norm reaches
421* stopping criterion.
422*
423 IF( nrm.LT.dtpcrt )
424 \$ GO TO 70
425 nrmchk = nrmchk + 1
426 IF( nrmchk.LT.extra+1 )
427 \$ GO TO 70
428*
429 GO TO 130
430*
431* If stopping criterion was not satisfied, update info and
432* store eigenvector number in array ifail.
433*
434 120 CONTINUE
435 info = info + 1
436 ifail( info ) = j
437*
438* Accept iterate as jth eigenvector.
439*
440 130 CONTINUE
441 scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
442 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
443 IF( work( indrv1+jmax ).LT.zero )
444 \$ scl = -scl
445 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
446 140 CONTINUE
447 DO 150 i = 1, n
448 z( i, j ) = czero
449 150 CONTINUE
450 DO 160 i = 1, blksiz
451 z( b1+i-1, j ) = dcmplx( work( indrv1+i ), zero )
452 160 CONTINUE
453*
454* Save the shift to check eigenvalue spacing at next
455* iteration.
456*
457 xjm = xj
458*
459 170 CONTINUE
460 180 CONTINUE
461*
462 RETURN
463*
464* End of ZSTEIN
465*
466 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlagtf(n, a, lambda, b, c, tol, d, in, info)
DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix,...
Definition dlagtf.f:156
subroutine dlagts(job, n, a, b, c, d, in, y, tol, info)
DLAGTS solves the system of equations (T-λI)x = y or (T-λI)^Tx = y, where T is a general tridiagonal ...
Definition dlagts.f:163
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
ZSTEIN
Definition zstein.f:182