179 SUBROUTINE cptrfs( UPLO, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
180 $ FERR, BERR, WORK, RWORK, INFO )
188 INTEGER INFO, LDB, LDX, N, NRHS
191 REAL BERR( * ), D( * ), DF( * ), FERR( * ),
193 COMPLEX B( LDB, * ), E( * ), EF( * ), WORK( * ),
201 parameter( itmax = 5 )
203 parameter( zero = 0.0e+0 )
205 parameter( one = 1.0e+0 )
207 parameter( two = 2.0e+0 )
209 parameter( three = 3.0e+0 )
213 INTEGER COUNT, I, IX, J, NZ
214 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
215 COMPLEX BI, CX, DX, EX, ZDUM
221 EXTERNAL lsame, isamax, slamch
227 INTRINSIC abs, aimag, cmplx, conjg, max, real
233 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
240 upper = lsame( uplo,
'U' )
241 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
243 ELSE IF( n.LT.0 )
THEN
245 ELSE IF( nrhs.LT.0 )
THEN
247 ELSE IF( ldb.LT.max( 1, n ) )
THEN
249 ELSE IF( ldx.LT.max( 1, n ) )
THEN
253 CALL xerbla(
'CPTRFS', -info )
259 IF( n.EQ.0 .OR. nrhs.EQ.0 )
THEN
270 eps = slamch(
'Epsilon' )
271 safmin = slamch(
'Safe minimum' )
272 safe1 = real( nz )*safmin
291 dx = d( 1 )*x( 1, j )
293 rwork( 1 ) = cabs1( bi ) + cabs1( dx )
296 dx = d( 1 )*x( 1, j )
297 ex = e( 1 )*x( 2, j )
298 work( 1 ) = bi - dx - ex
299 rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
300 $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
303 cx = conjg( e( i-1 ) )*x( i-1, j )
304 dx = d( i )*x( i, j )
305 ex = e( i )*x( i+1, j )
306 work( i ) = bi - cx - dx - ex
307 rwork( i ) = cabs1( bi ) +
308 $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
309 $ cabs1( dx ) + cabs1( e( i ) )*
310 $ cabs1( x( i+1, j ) )
313 cx = conjg( e( n-1 ) )*x( n-1, j )
314 dx = d( n )*x( n, j )
315 work( n ) = bi - cx - dx
316 rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
317 $ cabs1( x( n-1, j ) ) + cabs1( dx )
322 dx = d( 1 )*x( 1, j )
324 rwork( 1 ) = cabs1( bi ) + cabs1( dx )
327 dx = d( 1 )*x( 1, j )
328 ex = conjg( e( 1 ) )*x( 2, j )
329 work( 1 ) = bi - dx - ex
330 rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
331 $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
334 cx = e( i-1 )*x( i-1, j )
335 dx = d( i )*x( i, j )
336 ex = conjg( e( i ) )*x( i+1, j )
337 work( i ) = bi - cx - dx - ex
338 rwork( i ) = cabs1( bi ) +
339 $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
340 $ cabs1( dx ) + cabs1( e( i ) )*
341 $ cabs1( x( i+1, j ) )
344 cx = e( n-1 )*x( n-1, j )
345 dx = d( n )*x( n, j )
346 work( n ) = bi - cx - dx
347 rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
348 $ cabs1( x( n-1, j ) ) + cabs1( dx )
363 IF( rwork( i ).GT.safe2 )
THEN
364 s = max( s, cabs1( work( i ) ) / rwork( i ) )
366 s = max( s, ( cabs1( work( i ) )+safe1 ) /
367 $ ( rwork( i )+safe1 ) )
378 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
379 $ count.LE.itmax )
THEN
383 CALL cpttrs( uplo, n, 1, df, ef, work, n, info )
384 CALL caxpy( n, cmplx( one ), work, 1, x( 1, j ), 1 )
409 IF( rwork( i ).GT.safe2 )
THEN
410 rwork( i ) = cabs1( work( i ) ) + real( nz )*
413 rwork( i ) = cabs1( work( i ) ) + real( nz )*
414 $ eps*rwork( i ) + safe1
417 ix = isamax( n, rwork, 1 )
418 ferr( j ) = rwork( ix )
433 rwork( i ) = one + rwork( i-1 )*abs( ef( i-1 ) )
438 rwork( n ) = rwork( n ) / df( n )
439 DO 80 i = n - 1, 1, -1
440 rwork( i ) = rwork( i ) / df( i ) +
441 $ rwork( i+1 )*abs( ef( i ) )
446 ix = isamax( n, rwork, 1 )
447 ferr( j ) = ferr( j )*abs( rwork( ix ) )
453 lstres = max( lstres, abs( x( i, j ) ) )
456 $ ferr( j ) = ferr( j ) / lstres