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