174 SUBROUTINE dstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
175 $ iwork, ifail, info )
183 INTEGER info, ldz, m, n
186 INTEGER iblock( * ), ifail( * ), isplit( * ),
188 DOUBLE PRECISION d( * ), e( * ), w( * ), work( * ), z( ldz, * )
194 DOUBLE PRECISION zero, one, ten, odm3, odm1
195 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
196 $ odm3 = 1.0d-3, odm1 = 1.0d-1 )
197 INTEGER maxits, extra
198 parameter( maxits = 5, extra = 2 )
201 INTEGER b1, blksiz, bn, gpind, i, iinfo, indrv1,
202 $ indrv2, indrv3, indrv4, indrv5, its,
j, j1,
203 $ jblk, jmax, nblk, nrmchk
204 DOUBLE PRECISION dtpcrt, eps, eps1, nrm, onenrm, ortol, pertol,
205 $ scl, sep, tol, xj, xjm, ztr
220 INTRINSIC abs, max, sqrt
233 ELSE IF( m.LT.0 .OR. m.GT.n )
THEN
235 ELSE IF( ldz.LT.max( 1, n ) )
THEN
239 IF( iblock(
j ).LT.iblock(
j-1 ) )
THEN
243 IF( iblock(
j ).EQ.iblock(
j-1 ) .AND. w(
j ).LT.w(
j-1 ) )
253 CALL
xerbla(
'DSTEIN', -info )
259 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
261 ELSE IF( n.EQ.1 )
THEN
268 eps =
dlamch(
'Precision' )
287 DO 160 nblk = 1, iblock( m )
294 b1 = isplit( nblk-1 ) + 1
304 onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
305 onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
306 DO 50 i = b1 + 1, bn - 1
307 onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
312 dtpcrt = sqrt( odm1 / blksiz )
319 IF( iblock(
j ).NE.nblk )
THEN
328 IF( blksiz.EQ.1 )
THEN
329 work( indrv1+1 ) = one
349 CALL
dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
353 CALL
dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
354 CALL
dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
355 CALL
dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
360 CALL
dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
361 $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
373 scl = blksiz*onenrm*max( eps,
374 $ abs( work( indrv4+blksiz ) ) ) /
375 $
dasum( blksiz, work( indrv1+1 ), 1 )
376 CALL
dscal( blksiz, scl, work( indrv1+1 ), 1 )
380 CALL
dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
381 $ work( indrv3+1 ), work( indrv5+1 ), iwork,
382 $ work( indrv1+1 ), tol, iinfo )
389 IF( abs( xj-xjm ).GT.ortol )
391 IF( gpind.NE.
j )
THEN
392 DO 80 i = gpind,
j - 1
393 ztr = -
ddot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
395 CALL
daxpy( blksiz, ztr, z( b1, i ), 1,
396 $ work( indrv1+1 ), 1 )
403 jmax =
idamax( blksiz, work( indrv1+1 ), 1 )
404 nrm = abs( work( indrv1+jmax ) )
412 IF( nrmchk.LT.extra+1 )
427 scl = one /
dnrm2( blksiz, work( indrv1+1 ), 1 )
428 jmax =
idamax( blksiz, work( indrv1+1 ), 1 )
429 IF( work( indrv1+jmax ).LT.zero )
431 CALL
dscal( blksiz, scl, work( indrv1+1 ), 1 )
437 z( b1+i-1,
j ) = work( indrv1+i )
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
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 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
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
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...
integer function idamax(N, DX, INCX)
IDAMAX
double precision function dnrm2(N, X, INCX)
DNRM2
double precision function dasum(N, DX, INCX)
DASUM