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 DDOT, DLAMCH, DNRM2
213 EXTERNAL idamax, 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 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
374 scl = blksiz*onenrm*max( eps,
375 $ abs( work( indrv4+blksiz ) ) ) /
376 $ abs( work( indrv1+jmax ) )
377 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
381 CALL dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
382 $ work( indrv3+1 ), work( indrv5+1 ), iwork,
383 $ work( indrv1+1 ), tol, iinfo )
390 IF( abs( xj-xjm ).GT.ortol )
392 IF( gpind.NE.j )
THEN
393 DO 80 i = gpind, j - 1
394 ztr = -ddot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
396 CALL daxpy( blksiz, ztr, z( b1, i ), 1,
397 $ work( indrv1+1 ), 1 )
404 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
405 nrm = abs( work( indrv1+jmax ) )
413 IF( nrmchk.LT.extra+1 )
428 scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
429 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
430 IF( work( indrv1+jmax ).LT.zero )
432 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
438 z( b1+i-1, j ) = work( indrv1+i )
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...
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
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 xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL