174 SUBROUTINE sstein( 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 REAL d( * ), e( * ), w( * ), work( * ), z( ldz, * )
194 REAL zero, one, ten, odm3, odm1
195 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
196 $ odm3 = 1.0e-3, odm1 = 1.0e-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 REAL ctr, eps, eps1, nrm, onenrm, ortol, pertol,
205 $ scl, sep, stpcrt, tol, xj, xjm
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(
'SSTEIN', -info )
259 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
261 ELSE IF( n.EQ.1 )
THEN
268 eps =
slamch(
'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 stpcrt = sqrt( odm1 / blksiz )
319 IF( iblock(
j ).NE.nblk )
THEN
328 IF( blksiz.EQ.1 )
THEN
329 work( indrv1+1 ) = one
349 CALL
slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
353 CALL
scopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
354 CALL
scopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
355 CALL
scopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
360 CALL
slagtf( 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 $
sasum( blksiz, work( indrv1+1 ), 1 )
376 CALL
sscal( blksiz, scl, work( indrv1+1 ), 1 )
380 CALL
slagts( -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 ctr = -
sdot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
395 CALL
saxpy( blksiz, ctr, z( b1, i ), 1,
396 $ work( indrv1+1 ), 1 )
403 jmax =
isamax( blksiz, work( indrv1+1 ), 1 )
404 nrm = abs( work( indrv1+jmax ) )
412 IF( nrmchk.LT.extra+1 )
427 scl = one /
snrm2( blksiz, work( indrv1+1 ), 1 )
428 jmax =
isamax( blksiz, work( indrv1+1 ), 1 )
429 IF( work( indrv1+jmax ).LT.zero )
431 CALL
sscal( blksiz, scl, work( indrv1+1 ), 1 )
437 z( b1+i-1,
j ) = work( indrv1+i )
integer function isamax(N, SX, INCX)
ISAMAX
real function sasum(N, SX, INCX)
SASUM
real function sdot(N, SX, INCX, SY, INCY)
SDOT
subroutine slagtf(N, A, LAMBDA, B, C, TOL, D, IN, INFO)
SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
real function slamch(CMACH)
SLAMCH
real function snrm2(N, X, INCX)
SNRM2
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slagts(JOB, N, A, B, C, D, IN, Y, TOL, INFO)
SLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal ma...
subroutine sscal(N, SA, SX, INCX)
SSCAL