LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 *> \date December 2016
170 *
171 *> \ingroup doubleOTHERcomputational
172 *
173 * =====================================================================
174  SUBROUTINE dstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
175  $ IWORK, IFAIL, INFO )
176 *
177 * -- LAPACK computational routine (version 3.7.0) --
178 * -- LAPACK is a software package provided by Univ. of Tennessee, --
179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 * December 2016
181 *
182 * .. Scalar Arguments ..
183  INTEGER INFO, LDZ, M, N
184 * ..
185 * .. Array Arguments ..
186  INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
187  $ iwork( * )
188  DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
189 * ..
190 *
191 * =====================================================================
192 *
193 * .. Parameters ..
194  DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
195  parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
196  $ odm3 = 1.0d-3, odm1 = 1.0d-1 )
197  INTEGER MAXITS, EXTRA
198  parameter( maxits = 5, extra = 2 )
199 * ..
200 * .. Local Scalars ..
201  INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
202  $ indrv2, indrv3, indrv4, indrv5, its, j, j1,
203  $ jblk, jmax, nblk, nrmchk
204  DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
205  $ scl, sep, tol, xj, xjm, ztr
206 * ..
207 * .. Local Arrays ..
208  INTEGER ISEED( 4 )
209 * ..
210 * .. External Functions ..
211  INTEGER IDAMAX
212  DOUBLE PRECISION DDOT, DLAMCH, DNRM2
213  EXTERNAL idamax, ddot, dlamch, dnrm2
214 * ..
215 * .. External Subroutines ..
216  EXTERNAL daxpy, dcopy, dlagtf, dlagts, dlarnv, dscal,
217  $ xerbla
218 * ..
219 * .. Intrinsic Functions ..
220  INTRINSIC abs, max, sqrt
221 * ..
222 * .. Executable Statements ..
223 *
224 * Test the input parameters.
225 *
226  info = 0
227  DO 10 i = 1, m
228  ifail( i ) = 0
229  10 CONTINUE
230 *
231  IF( n.LT.0 ) THEN
232  info = -1
233  ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
234  info = -4
235  ELSE IF( ldz.LT.max( 1, n ) ) THEN
236  info = -9
237  ELSE
238  DO 20 j = 2, m
239  IF( iblock( j ).LT.iblock( j-1 ) ) THEN
240  info = -6
241  GO TO 30
242  END IF
243  IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
244  $ THEN
245  info = -5
246  GO TO 30
247  END IF
248  20 CONTINUE
249  30 CONTINUE
250  END IF
251 *
252  IF( info.NE.0 ) THEN
253  CALL xerbla( 'DSTEIN', -info )
254  RETURN
255  END IF
256 *
257 * Quick return if possible
258 *
259  IF( n.EQ.0 .OR. m.EQ.0 ) THEN
260  RETURN
261  ELSE IF( n.EQ.1 ) THEN
262  z( 1, 1 ) = one
263  RETURN
264  END IF
265 *
266 * Get machine constants.
267 *
268  eps = dlamch( 'Precision' )
269 *
270 * Initialize seed for random number generator DLARNV.
271 *
272  DO 40 i = 1, 4
273  iseed( i ) = 1
274  40 CONTINUE
275 *
276 * Initialize pointers.
277 *
278  indrv1 = 0
279  indrv2 = indrv1 + n
280  indrv3 = indrv2 + n
281  indrv4 = indrv3 + n
282  indrv5 = indrv4 + n
283 *
284 * Compute eigenvectors of matrix blocks.
285 *
286  j1 = 1
287  DO 160 nblk = 1, iblock( m )
288 *
289 * Find starting and ending indices of block nblk.
290 *
291  IF( nblk.EQ.1 ) THEN
292  b1 = 1
293  ELSE
294  b1 = isplit( nblk-1 ) + 1
295  END IF
296  bn = isplit( nblk )
297  blksiz = bn - b1 + 1
298  IF( blksiz.EQ.1 )
299  $ GO TO 60
300  gpind = j1
301 *
302 * Compute reorthogonalization criterion and stopping criterion.
303 *
304  onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
305  onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
306  DO 50 i = b1 + 1, bn - 1
307  onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
308  $ abs( e( i ) ) )
309  50 CONTINUE
310  ortol = odm3*onenrm
311 *
312  dtpcrt = sqrt( odm1 / blksiz )
313 *
314 * Loop through eigenvalues of block nblk.
315 *
316  60 CONTINUE
317  jblk = 0
318  DO 150 j = j1, m
319  IF( iblock( j ).NE.nblk ) THEN
320  j1 = j
321  GO TO 160
322  END IF
323  jblk = jblk + 1
324  xj = w( j )
325 *
326 * Skip all the work if the block size is one.
327 *
328  IF( blksiz.EQ.1 ) THEN
329  work( indrv1+1 ) = one
330  GO TO 120
331  END IF
332 *
333 * If eigenvalues j and j-1 are too close, add a relatively
334 * small perturbation.
335 *
336  IF( jblk.GT.1 ) THEN
337  eps1 = abs( eps*xj )
338  pertol = ten*eps1
339  sep = xj - xjm
340  IF( sep.LT.pertol )
341  $ xj = xjm + pertol
342  END IF
343 *
344  its = 0
345  nrmchk = 0
346 *
347 * Get random starting vector.
348 *
349  CALL dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
350 *
351 * Copy the matrix T so it won't be destroyed in factorization.
352 *
353  CALL dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
354  CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
355  CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
356 *
357 * Compute LU factors with partial pivoting ( PT = LU )
358 *
359  tol = zero
360  CALL dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
361  $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
362  $ iinfo )
363 *
364 * Update iteration count.
365 *
366  70 CONTINUE
367  its = its + 1
368  IF( its.GT.maxits )
369  $ GO TO 100
370 *
371 * Normalize and scale the righthand side vector Pb.
372 *
373  jmax = idamax( blksiz, work( indrv1+1 ), 1 )
374  scl = blksiz*onenrm*max( eps,
375  $ abs( work( indrv4+blksiz ) ) ) /
376  $ abs( work( indrv1+jmax ) )
377  CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
378 *
379 * Solve the system LU = Pb.
380 *
381  CALL dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
382  $ work( indrv3+1 ), work( indrv5+1 ), iwork,
383  $ work( indrv1+1 ), tol, iinfo )
384 *
385 * Reorthogonalize by modified Gram-Schmidt if eigenvalues are
386 * close enough.
387 *
388  IF( jblk.EQ.1 )
389  $ GO TO 90
390  IF( abs( xj-xjm ).GT.ortol )
391  $ gpind = j
392  IF( gpind.NE.j ) THEN
393  DO 80 i = gpind, j - 1
394  ztr = -ddot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
395  $ 1 )
396  CALL daxpy( blksiz, ztr, z( b1, i ), 1,
397  $ work( indrv1+1 ), 1 )
398  80 CONTINUE
399  END IF
400 *
401 * Check the infinity norm of the iterate.
402 *
403  90 CONTINUE
404  jmax = idamax( blksiz, work( indrv1+1 ), 1 )
405  nrm = abs( work( indrv1+jmax ) )
406 *
407 * Continue for additional iterations after norm reaches
408 * stopping criterion.
409 *
410  IF( nrm.LT.dtpcrt )
411  $ GO TO 70
412  nrmchk = nrmchk + 1
413  IF( nrmchk.LT.extra+1 )
414  $ GO TO 70
415 *
416  GO TO 110
417 *
418 * If stopping criterion was not satisfied, update info and
419 * store eigenvector number in array ifail.
420 *
421  100 CONTINUE
422  info = info + 1
423  ifail( info ) = j
424 *
425 * Accept iterate as jth eigenvector.
426 *
427  110 CONTINUE
428  scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
429  jmax = idamax( blksiz, work( indrv1+1 ), 1 )
430  IF( work( indrv1+jmax ).LT.zero )
431  $ scl = -scl
432  CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
433  120 CONTINUE
434  DO 130 i = 1, n
435  z( i, j ) = zero
436  130 CONTINUE
437  DO 140 i = 1, blksiz
438  z( b1+i-1, j ) = work( indrv1+i )
439  140 CONTINUE
440 *
441 * Save the shift to check eigenvalue spacing at next
442 * iteration.
443 *
444  xjm = xj
445 *
446  150 CONTINUE
447  160 CONTINUE
448 *
449  RETURN
450 *
451 * End of DSTEIN
452 *
453  END
dstein
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:176
dlagts
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 ma...
Definition: dlagts.f:163
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
dlarnv
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlagtf
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:158
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91