182 SUBROUTINE cstein( 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 REAL d( * ), e( * ), w( * ), work( * )
204 parameter( czero = ( 0.0e+0, 0.0e+0 ),
205 $ cone = ( 1.0e+0, 0.0e+0 ) )
206 REAL zero, one, ten, odm3, odm1
207 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
208 $ odm3 = 1.0e-3, odm1 = 1.0e-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 REAL ctr, eps, eps1, nrm, onenrm, ortol, pertol,
217 $ scl, sep, stpcrt, tol, xj, xjm
231 INTRINSIC abs, cmplx, max,
REAL, 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(
'CSTEIN', -info )
270 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
272 ELSE IF( n.EQ.1 )
THEN
279 eps =
slamch(
'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 stpcrt = sqrt( odm1 / blksiz )
330 IF( iblock(
j ).NE.nblk )
THEN
339 IF( blksiz.EQ.1 )
THEN
340 work( indrv1+1 ) = one
360 CALL
slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
364 CALL
scopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
365 CALL
scopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
366 CALL
scopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
371 CALL
slagtf( 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 $
sasum( blksiz, work( indrv1+1 ), 1 )
387 CALL
sscal( blksiz, scl, work( indrv1+1 ), 1 )
391 CALL
slagts( -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 ctr = ctr + work( indrv1+jr )*
407 $
REAL( Z( B1-1+JR, I ) )
410 work( indrv1+jr ) = work( indrv1+jr ) -
411 $ ctr*
REAL( Z( B1-1+JR, I ) )
419 jmax =
isamax( blksiz, work( indrv1+1 ), 1 )
420 nrm = abs( work( indrv1+jmax ) )
428 IF( nrmchk.LT.extra+1 )
443 scl = one /
snrm2( blksiz, work( indrv1+1 ), 1 )
444 jmax =
isamax( blksiz, work( indrv1+1 ), 1 )
445 IF( work( indrv1+jmax ).LT.zero )
447 CALL
sscal( blksiz, scl, work( indrv1+1 ), 1 )
453 z( b1+i-1,
j ) = cmplx( work( indrv1+i ), zero )
INTEGER function isamax(N, SX, INCX)
ISAMAX
REAL function slamch(CMACH)
SLAMCH
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 cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
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
REAL function sasum(N, SX, INCX)
SASUM
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...
REAL function snrm2(N, X, INCX)
SNRM2
subroutine sscal(N, SA, SX, INCX)
SSCAL