180 SUBROUTINE cstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
181 $ IWORK, IFAIL, INFO )
188 INTEGER INFO, LDZ, M, N
191 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
193 REAL D( * ), E( * ), W( * ), WORK( * )
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 )
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
222 EXTERNAL isamax, slamch, snrm2
228 INTRINSIC abs, cmplx, max, real, sqrt
241 ELSE IF( m.LT.0 .OR. m.GT.n )
THEN
243 ELSE IF( ldz.LT.max( 1, n ) )
THEN
247 IF( iblock( j ).LT.iblock( j-1 ) )
THEN
251 IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
261 CALL xerbla(
'CSTEIN', -info )
267 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
269 ELSE IF( n.EQ.1 )
THEN
276 eps = slamch(
'Precision' )
295 DO 180 nblk = 1, iblock( m )
302 b1 = isplit( nblk-1 ) + 1
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 ) )+
320 stpcrt = sqrt( odm1 / blksiz )
327 IF( iblock( j ).NE.nblk )
THEN
336 IF( blksiz.EQ.1 )
THEN
337 work( indrv1+1 ) = one
357 CALL slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
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 )
368 CALL slagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
369 $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
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 )
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 )
398 IF( abs( xj-xjm ).GT.ortol )
400 IF( gpind.NE.j )
THEN
401 DO 100 i = gpind, j - 1
404 ctr = ctr + work( indrv1+jr )*
405 $ real( z( b1-1+jr, i ) )
408 work( indrv1+jr ) = work( indrv1+jr ) -
409 $ ctr*real( z( b1-1+jr, i ) )
417 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
418 nrm = abs( work( indrv1+jmax ) )
426 IF( nrmchk.LT.extra+1 )
441 scl = one / snrm2( blksiz, work( indrv1+1 ), 1 )
442 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
443 IF( work( indrv1+jmax ).LT.zero )
445 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
451 z( b1+i-1, j ) = cmplx( work( indrv1+i ), zero )
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
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...
subroutine xerbla(SRNAME, INFO)
XERBLA
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,...
subroutine cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sscal(N, SA, SX, INCX)
SSCAL