LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slagtf
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
slagts
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 ...
Definition slagts.f:163
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slarnv
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: