LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
cstein.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine cstein (N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
 CSTEIN More...
 

Function/Subroutine Documentation

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
November 2011

Definition at line 184 of file cstein.f.

184 *
185 * -- LAPACK computational routine (version 3.4.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 * November 2011
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 sasum, slamch, snrm2
225  EXTERNAL isamax, sasum, 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 = b1
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  scl = blksiz*onenrm*max( eps,
385  $ abs( work( indrv4+blksiz ) ) ) /
386  $ sasum( blksiz, work( indrv1+1 ), 1 )
387  CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
388 *
389 * Solve the system LU = Pb.
390 *
391  CALL slagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
392  $ work( indrv3+1 ), work( indrv5+1 ), iwork,
393  $ work( indrv1+1 ), tol, iinfo )
394 *
395 * Reorthogonalize by modified Gram-Schmidt if eigenvalues are
396 * close enough.
397 *
398  IF( jblk.EQ.1 )
399  $ GO TO 110
400  IF( abs( xj-xjm ).GT.ortol )
401  $ gpind = j
402  IF( gpind.NE.j ) THEN
403  DO 100 i = gpind, j - 1
404  ctr = zero
405  DO 80 jr = 1, blksiz
406  ctr = ctr + work( indrv1+jr )*
407  $ REAL( Z( B1-1+JR, I ) )
408  80 CONTINUE
409  DO 90 jr = 1, blksiz
410  work( indrv1+jr ) = work( indrv1+jr ) -
411  $ ctr*REAL( Z( B1-1+JR, I ) )
412  90 CONTINUE
413  100 CONTINUE
414  END IF
415 *
416 * Check the infinity norm of the iterate.
417 *
418  110 CONTINUE
419  jmax = isamax( blksiz, work( indrv1+1 ), 1 )
420  nrm = abs( work( indrv1+jmax ) )
421 *
422 * Continue for additional iterations after norm reaches
423 * stopping criterion.
424 *
425  IF( nrm.LT.stpcrt )
426  $ GO TO 70
427  nrmchk = nrmchk + 1
428  IF( nrmchk.LT.extra+1 )
429  $ GO TO 70
430 *
431  GO TO 130
432 *
433 * If stopping criterion was not satisfied, update info and
434 * store eigenvector number in array ifail.
435 *
436  120 CONTINUE
437  info = info + 1
438  ifail( info ) = j
439 *
440 * Accept iterate as jth eigenvector.
441 *
442  130 CONTINUE
443  scl = one / snrm2( blksiz, work( indrv1+1 ), 1 )
444  jmax = isamax( blksiz, work( indrv1+1 ), 1 )
445  IF( work( indrv1+jmax ).LT.zero )
446  $ scl = -scl
447  CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
448  140 CONTINUE
449  DO 150 i = 1, n
450  z( i, j ) = czero
451  150 CONTINUE
452  DO 160 i = 1, blksiz
453  z( b1+i-1, j ) = cmplx( work( indrv1+i ), zero )
454  160 CONTINUE
455 *
456 * Save the shift to check eigenvalue spacing at next
457 * iteration.
458 *
459  xjm = xj
460 *
461  170 CONTINUE
462  180 CONTINUE
463 *
464  RETURN
465 *
466 * End of CSTEIN
467 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
real function sasum(N, SX, INCX)
SASUM
Definition: sasum.f:54
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:99
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
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: