LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ cstein()

subroutine cstein ( integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
integer  M,
real, dimension( * )  W,
integer, dimension( * )  IBLOCK,
integer, dimension( * )  ISPLIT,
complex, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

CSTEIN

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

Purpose:
 CSTEIN 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 CUNMTR or CUPMTR 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 REAL array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E
          E is REAL 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 REAL 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 SSTEBZ 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 SSTEBZ 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 SSTEBZ is expected here. )
[out]Z
          Z is COMPLEX 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 REAL 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 cstein.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  REAL d( * ), e( * ), w( * ), work( * )
197  COMPLEX z( ldz, * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  COMPLEX czero, cone
204  parameter( czero = ( 0.0e+0, 0.0e+0 ),
205  $ cone = ( 1.0e+0, 0.0e+0 ) )
206  REAL zero, one, ten, odm3, odm1
207  parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
208  $ odm3 = 1.0e-3, odm1 = 1.0e-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  REAL ctr, eps, eps1, nrm, onenrm, ortol, pertol,
217  $ scl, sep, stpcrt, tol, xj, xjm
218 * ..
219 * .. Local Arrays ..
220  INTEGER iseed( 4 )
221 * ..
222 * .. External Functions ..
223  INTEGER isamax
224  REAL slamch, snrm2
225  EXTERNAL isamax, slamch, snrm2
226 * ..
227 * .. External Subroutines ..
228  EXTERNAL scopy, slagtf, slagts, slarnv, sscal, xerbla
229 * ..
230 * .. Intrinsic Functions ..
231  INTRINSIC abs, cmplx, max, REAL, 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( 'CSTEIN', -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 = slamch( 'Precision' )
280 *
281 * Initialize seed for random number generator SLARNV.
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  stpcrt = 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 slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
361 *
362 * Copy the matrix T so it won't be destroyed in factorization.
363 *
364  CALL scopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
365  CALL scopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
366  CALL scopy( 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 slagtf( 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 = isamax( blksiz, work( indrv1+1 ), 1 )
385  scl = blksiz*onenrm*max( eps,
386  $ abs( work( indrv4+blksiz ) ) ) /
387  $ abs( work( indrv1+jmax ) )
388  CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
389 *
390 * Solve the system LU = Pb.
391 *
392  CALL slagts( -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  ctr = zero
406  DO 80 jr = 1, blksiz
407  ctr = ctr + work( indrv1+jr )*
408  $ REAL( Z( B1-1+JR, I ) )
409  80 CONTINUE
410  DO 90 jr = 1, blksiz
411  work( indrv1+jr ) = work( indrv1+jr ) -
412  $ ctr*REAL( 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 = isamax( 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.stpcrt )
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 / snrm2( blksiz, work( indrv1+1 ), 1 )
445  jmax = isamax( blksiz, work( indrv1+1 ), 1 )
446  IF( work( indrv1+jmax ).LT.zero )
447  $ scl = -scl
448  CALL sscal( 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 ) = cmplx( 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 CSTEIN
468 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:73
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slagts(JOB, N, A, B, C, D, IN, Y, TOL, INFO)
SLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal ma...
Definition: slagts.f:163
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slagtf(N, A, LAMBDA, B, C, TOL, D, IN, INFO)
SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix...
Definition: slagtf.f:158
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:99
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
Here is the call graph for this function:
Here is the caller graph for this function: