180 SUBROUTINE zstein( 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 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
194 COMPLEX*16 Z( LDZ, * )
200 COMPLEX*16 CZERO, CONE
201 parameter( czero = ( 0.0d+0, 0.0d+0 ),
202 $ cone = ( 1.0d+0, 0.0d+0 ) )
203 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
204 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
205 $ odm3 = 1.0d-3, odm1 = 1.0d-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 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
214 $ scl, sep, tol, xj, xjm, ztr
221 DOUBLE PRECISION DLAMCH, DNRM2
222 EXTERNAL idamax, dlamch, dnrm2
228 INTRINSIC abs, dble, dcmplx, max, 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(
'ZSTEIN', -info )
267 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
269 ELSE IF( n.EQ.1 )
THEN
276 eps = dlamch(
'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 dtpcrt = sqrt( odm1 / blksiz )
327 IF( iblock( j ).NE.nblk )
THEN
336 IF( blksiz.EQ.1 )
THEN
337 work( indrv1+1 ) = one
357 CALL dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
361 CALL dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
362 CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
363 CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
368 CALL dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
369 $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
381 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
382 scl = blksiz*onenrm*max( eps,
383 $ abs( work( indrv4+blksiz ) ) ) /
384 $ abs( work( indrv1+jmax ) )
385 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
389 CALL dlagts( -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 ztr = ztr + work( indrv1+jr )*
405 $ dble( z( b1-1+jr, i ) )
408 work( indrv1+jr ) = work( indrv1+jr ) -
409 $ ztr*dble( z( b1-1+jr, i ) )
417 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
418 nrm = abs( work( indrv1+jmax ) )
426 IF( nrmchk.LT.extra+1 )
441 scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
442 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
443 IF( work( indrv1+jmax ).LT.zero )
445 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
451 z( b1+i-1, j ) = dcmplx( work( indrv1+i ), zero )
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 dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine xerbla(SRNAME, INFO)
XERBLA
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 zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dscal(N, DA, DX, INCX)
DSCAL