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...