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 DLAMCH, DNRM2
225 EXTERNAL idamax, 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 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 dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
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 xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL