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