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 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
385 scl = blksiz*onenrm*max( eps,
386 $ abs( work( indrv4+blksiz ) ) ) /
387 $ abs( work( indrv1+jmax ) )
388 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
392 CALL dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
393 $ work( indrv3+1 ), work( indrv5+1 ), iwork,
394 $ work( indrv1+1 ), tol, iinfo )
401 IF( abs( xj-xjm ).GT.ortol )
403 IF( gpind.NE.j )
THEN
404 DO 100 i = gpind, j - 1
407 ztr = ztr + work( indrv1+jr )*
408 $ dble( z( b1-1+jr, i ) )
411 work( indrv1+jr ) = work( indrv1+jr ) -
412 $ ztr*dble( z( b1-1+jr, i ) )
420 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
421 nrm = abs( work( indrv1+jmax ) )
429 IF( nrmchk.LT.extra+1 )
444 scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
445 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
446 IF( work( indrv1+jmax ).LT.zero )
448 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
454 z( b1+i-1, j ) = dcmplx( work( indrv1+i ), zero )
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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 zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
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...
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.