174 SUBROUTINE dstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
175 $ iwork, ifail, info )
183 INTEGER INFO, LDZ, M, N
186 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
188 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( ldz, * )
194 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
195 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
196 $ odm3 = 1.0d-3, odm1 = 1.0d-1 )
197 INTEGER MAXITS, EXTRA
198 parameter( maxits = 5, extra = 2 )
201 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
202 $ indrv2, indrv3, indrv4, indrv5, its, j, j1,
203 $ jblk, jmax, nblk, nrmchk
204 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
205 $ scl, sep, tol, xj, xjm, ztr
212 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DNRM2
213 EXTERNAL idamax, dasum, ddot, dlamch, dnrm2
220 INTRINSIC abs, max, sqrt
233 ELSE IF( m.LT.0 .OR. m.GT.n )
THEN
235 ELSE IF( ldz.LT.max( 1, n ) )
THEN
239 IF( iblock( j ).LT.iblock( j-1 ) )
THEN
243 IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
253 CALL xerbla(
'DSTEIN', -info )
259 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
261 ELSE IF( n.EQ.1 )
THEN
268 eps = dlamch(
'Precision' )
287 DO 160 nblk = 1, iblock( m )
294 b1 = isplit( nblk-1 ) + 1
304 onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
305 onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
306 DO 50 i = b1 + 1, bn - 1
307 onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
312 dtpcrt = sqrt( odm1 / blksiz )
319 IF( iblock( j ).NE.nblk )
THEN
328 IF( blksiz.EQ.1 )
THEN
329 work( indrv1+1 ) = one
349 CALL dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
353 CALL dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
354 CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
355 CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
360 CALL dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
361 $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
373 scl = blksiz*onenrm*max( eps,
374 $ abs( work( indrv4+blksiz ) ) ) /
375 $ dasum( blksiz, work( indrv1+1 ), 1 )
376 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
380 CALL dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
381 $ work( indrv3+1 ), work( indrv5+1 ), iwork,
382 $ work( indrv1+1 ), tol, iinfo )
389 IF( abs( xj-xjm ).GT.ortol )
391 IF( gpind.NE.j )
THEN
392 DO 80 i = gpind, j - 1
393 ztr = -ddot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
395 CALL daxpy( blksiz, ztr, z( b1, i ), 1,
396 $ work( indrv1+1 ), 1 )
403 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
404 nrm = abs( work( indrv1+jmax ) )
412 IF( nrmchk.LT.extra+1 )
427 scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
428 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
429 IF( work( indrv1+jmax ).LT.zero )
431 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
437 z( b1+i-1, j ) = work( indrv1+i )
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
subroutine dlagtf(N, A, LAMBDA, B, C, TOL, D, IN, INFO)
DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix...
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dlagts(JOB, N, A, B, C, D, IN, Y, TOL, INFO)
DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal ma...