182 SUBROUTINE zstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
183 $ iwork, ifail, info )
191 INTEGER INFO, LDZ, M, N
194 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
196 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
197 COMPLEX*16 Z( ldz, * )
203 COMPLEX*16 CZERO, CONE
204 parameter( czero = ( 0.0d+0, 0.0d+0 ),
205 $ cone = ( 1.0d+0, 0.0d+0 ) )
206 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
207 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
208 $ odm3 = 1.0d-3, odm1 = 1.0d-1 )
209 INTEGER MAXITS, EXTRA
210 parameter( maxits = 5, extra = 2 )
213 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
214 $ indrv2, indrv3, indrv4, indrv5, its, j, j1,
215 $ jblk, jmax, jr, nblk, nrmchk
216 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
217 $ scl, sep, tol, xj, xjm, ztr
224 DOUBLE PRECISION DASUM, DLAMCH, DNRM2
225 EXTERNAL idamax, dasum, dlamch, dnrm2
231 INTRINSIC abs, dble, dcmplx, max, sqrt
244 ELSE IF( m.LT.0 .OR. m.GT.n )
THEN
246 ELSE IF( ldz.LT.max( 1, n ) )
THEN
250 IF( iblock( j ).LT.iblock( j-1 ) )
THEN
254 IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
264 CALL
xerbla(
'ZSTEIN', -info )
270 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
272 ELSE IF( n.EQ.1 )
THEN
279 eps = dlamch(
'Precision' )
298 DO 180 nblk = 1, iblock( m )
305 b1 = isplit( nblk-1 ) + 1
315 onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
316 onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
317 DO 50 i = b1 + 1, bn - 1
318 onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
323 dtpcrt = sqrt( odm1 / blksiz )
330 IF( iblock( j ).NE.nblk )
THEN
339 IF( blksiz.EQ.1 )
THEN
340 work( indrv1+1 ) = one
360 CALL
dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
364 CALL
dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
365 CALL
dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
366 CALL
dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
371 CALL
dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
372 $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
384 scl = blksiz*onenrm*max( eps,
385 $ abs( work( indrv4+blksiz ) ) ) /
386 $ dasum( blksiz, work( indrv1+1 ), 1 )
387 CALL
dscal( blksiz, scl, work( indrv1+1 ), 1 )
391 CALL
dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
392 $ work( indrv3+1 ), work( indrv5+1 ), iwork,
393 $ work( indrv1+1 ), tol, iinfo )
400 IF( abs( xj-xjm ).GT.ortol )
402 IF( gpind.NE.j )
THEN
403 DO 100 i = gpind, j - 1
406 ztr = ztr + work( indrv1+jr )*
407 $ dble( z( b1-1+jr, i ) )
410 work( indrv1+jr ) = work( indrv1+jr ) -
411 $ ztr*dble( z( b1-1+jr, i ) )
419 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
420 nrm = abs( work( indrv1+jmax ) )
428 IF( nrmchk.LT.extra+1 )
443 scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
444 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
445 IF( work( indrv1+jmax ).LT.zero )
447 CALL
dscal( blksiz, scl, work( indrv1+1 ), 1 )
453 z( b1+i-1, j ) = dcmplx( work( indrv1+i ), zero )
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 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 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...