183 SUBROUTINE dporfs( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X,
184 $ ldx, ferr, berr, work, iwork, info )
193 INTEGER info, lda, ldaf, ldb, ldx, n, nrhs
197 DOUBLE PRECISION a( lda, * ), af( ldaf, * ), b( ldb, * ),
198 $ berr( * ), ferr( * ), work( * ), x( ldx, * )
205 parameter( itmax = 5 )
206 DOUBLE PRECISION zero
207 parameter( zero = 0.0d+0 )
209 parameter( one = 1.0d+0 )
211 parameter( two = 2.0d+0 )
212 DOUBLE PRECISION three
213 parameter( three = 3.0d+0 )
217 INTEGER count, i, j, k, kase, nz
218 DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
239 upper =
lsame( uplo,
'U' )
240 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
242 ELSE IF( n.LT.0 )
THEN
244 ELSE IF( nrhs.LT.0 )
THEN
246 ELSE IF( lda.LT.max( 1, n ) )
THEN
248 ELSE IF( ldaf.LT.max( 1, n ) )
THEN
250 ELSE IF( ldb.LT.max( 1, n ) )
THEN
252 ELSE IF( ldx.LT.max( 1, n ) )
THEN
256 CALL
xerbla(
'DPORFS', -info )
262 IF( n.EQ.0 .OR. nrhs.EQ.0 )
THEN
274 safmin =
dlamch(
'Safe minimum' )
290 CALL
dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
291 CALL
dsymv( uplo, n, -one, a, lda, x( 1, j ), 1, one,
304 work( i ) = abs( b( i, j ) )
312 xk = abs( x( k, j ) )
314 work( i ) = work( i ) + abs( a( i, k ) )*xk
315 s = s + abs( a( i, k ) )*abs( x( i, j ) )
317 work( k ) = work( k ) + abs( a( k, k ) )*xk + s
322 xk = abs( x( k, j ) )
323 work( k ) = work( k ) + abs( a( k, k ) )*xk
325 work( i ) = work( i ) + abs( a( i, k ) )*xk
326 s = s + abs( a( i, k ) )*abs( x( i, j ) )
328 work( k ) = work( k ) + s
333 IF( work( i ).GT.safe2 )
THEN
334 s = max( s, abs( work( n+i ) ) / work( i ) )
336 s = max( s, ( abs( work( n+i ) )+safe1 ) /
337 $ ( work( i )+safe1 ) )
348 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
349 $ count.LE.itmax )
THEN
353 CALL
dpotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
354 CALL
daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
383 IF( work( i ).GT.safe2 )
THEN
384 work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
386 work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
392 CALL
dlacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
399 CALL
dpotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
401 work( n+i ) = work( i )*work( n+i )
403 ELSE IF( kase.EQ.2 )
THEN
408 work( n+i ) = work( i )*work( n+i )
410 CALL
dpotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
419 lstres = max( lstres, abs( x( i, j ) ) )
422 $ ferr( j ) = ferr( j ) / lstres