LAPACK  3.10.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.

Definition at line 180 of file cstein.f.

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