LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zstein()

subroutine zstein ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
integer  M,
double precision, dimension( * )  W,
integer, dimension( * )  IBLOCK,
integer, dimension( * )  ISPLIT,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

ZSTEIN

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

Purpose:
 ZSTEIN computes the eigenvectors of a real symmetric tridiagonal
 matrix T corresponding to specified eigenvalues, using inverse
 iteration.

 The maximum number of iterations allowed for each eigenvector is
 specified by an internal parameter MAXITS (currently set to 5).

 Although the eigenvectors are real, they are stored in a complex
 array, which may be passed to ZUNMTR or ZUPMTR for back
 transformation to the eigenvectors of a complex Hermitian matrix
 which was reduced to tridiagonal form.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of the tridiagonal matrix
          T, stored in elements 1 to N-1.
[in]M
          M is INTEGER
          The number of eigenvectors to be found.  0 <= M <= N.
[in]W
          W is DOUBLE PRECISION array, dimension (N)
          The first M elements of W contain the eigenvalues for
          which eigenvectors are to be computed.  The eigenvalues
          should be grouped by split-off block and ordered from
          smallest to largest within the block.  ( The output array
          W from DSTEBZ with ORDER = 'B' is expected here. )
[in]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          The submatrix indices associated with the corresponding
          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
          the first submatrix from the top, =2 if W(i) belongs to
          the second submatrix, etc.  ( The output array IBLOCK
          from DSTEBZ is expected here. )
[in]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into submatrices.
          The first submatrix consists of rows/columns 1 to
          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
          through ISPLIT( 2 ), etc.
          ( The output array ISPLIT from DSTEBZ is expected here. )
[out]Z
          Z is COMPLEX*16 array, dimension (LDZ, M)
          The computed eigenvectors.  The eigenvector associated
          with the eigenvalue W(i) is stored in the i-th column of
          Z.  Any vector which fails to converge is set to its current
          iterate after MAXITS iterations.
          The imaginary parts of the eigenvectors are set to zero.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= max(1,N).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (5*N)
[out]IWORK
          IWORK is INTEGER array, dimension (N)
[out]IFAIL
          IFAIL is INTEGER array, dimension (M)
          On normal exit, all elements of IFAIL are zero.
          If one or more eigenvectors fail to converge after
          MAXITS iterations, then their indices are stored in
          array IFAIL.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, then i eigenvectors failed to converge
               in MAXITS iterations.  Their indices are stored in
               array IFAIL.
Internal Parameters:
  MAXITS  INTEGER, default = 5
          The maximum number of iterations performed.

  EXTRA   INTEGER, default = 2
          The number of iterations performed after norm growth
          criterion is satisfied, should be at least 1.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 184 of file zstein.f.

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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:73
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
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
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
Here is the call graph for this function:
Here is the caller graph for this function: