114 SUBROUTINE dsptrs( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
122 INTEGER INFO, LDB, N, NRHS
126 DOUBLE PRECISION AP( * ), B( LDB, * )
133 parameter( one = 1.0d+0 )
138 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
153 upper = lsame( uplo,
'U' )
154 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
156 ELSE IF( n.LT.0 )
THEN
158 ELSE IF( nrhs.LT.0 )
THEN
160 ELSE IF( ldb.LT.max( 1, n ) )
THEN
164 CALL xerbla(
'DSPTRS', -info )
170 IF( n.EQ.0 .OR. nrhs.EQ.0 )
183 kc = n*( n+1 ) / 2 + 1
192 IF( ipiv( k ).GT.0 )
THEN
200 $
CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
205 CALL dger( k-1, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
210 CALL dscal( nrhs, one / ap( kc+k-1 ), b( k, 1 ), ldb )
220 $
CALL dswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
225 CALL dger( k-2, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
227 CALL dger( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1,
228 $ b( k-1, 1 ), ldb, b( 1, 1 ), ldb )
233 akm1 = ap( kc-1 ) / akm1k
234 ak = ap( kc+k-1 ) / akm1k
235 denom = akm1*ak - one
237 bkm1 = b( k-1, j ) / akm1k
238 bk = b( k, j ) / akm1k
239 b( k-1, j ) = ( ak*bkm1-bk ) / denom
240 b( k, j ) = ( akm1*bk-bkm1 ) / denom
263 IF( ipiv( k ).GT.0 )
THEN
270 CALL dgemv(
'Transpose', k-1, nrhs, -one, b, ldb, ap( kc ),
271 $ 1, one, b( k, 1 ), ldb )
277 $
CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
287 CALL dgemv(
'Transpose', k-1, nrhs, -one, b, ldb, ap( kc ),
288 $ 1, one, b( k, 1 ), ldb )
289 CALL dgemv(
'Transpose', k-1, nrhs, -one, b, ldb,
290 $ ap( kc+k ), 1, one, b( k+1, 1 ), ldb )
296 $
CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
322 IF( ipiv( k ).GT.0 )
THEN
330 $
CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
336 $
CALL dger( n-k, nrhs, -one, ap( kc+1 ), 1, b( k, 1 ),
337 $ ldb, b( k+1, 1 ), ldb )
341 CALL dscal( nrhs, one / ap( kc ), b( k, 1 ), ldb )
352 $
CALL dswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
358 CALL dger( n-k-1, nrhs, -one, ap( kc+2 ), 1, b( k, 1 ),
359 $ ldb, b( k+2, 1 ), ldb )
360 CALL dger( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1,
361 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
367 akm1 = ap( kc ) / akm1k
368 ak = ap( kc+n-k+1 ) / akm1k
369 denom = akm1*ak - one
371 bkm1 = b( k, j ) / akm1k
372 bk = b( k+1, j ) / akm1k
373 b( k, j ) = ( ak*bkm1-bk ) / denom
374 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
376 kc = kc + 2*( n-k ) + 1
389 kc = n*( n+1 ) / 2 + 1
398 IF( ipiv( k ).GT.0 )
THEN
406 $
CALL dgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
407 $ ldb, ap( kc+1 ), 1, one, b( k, 1 ), ldb )
413 $
CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
423 CALL dgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
424 $ ldb, ap( kc+1 ), 1, one, b( k, 1 ), ldb )
425 CALL dgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
426 $ ldb, ap( kc-( n-k ) ), 1, one, b( k-1, 1 ),
434 $
CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
subroutine xerbla(srname, info)
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
subroutine dsptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
DSPTRS
subroutine dscal(n, da, dx, incx)
DSCAL
subroutine dswap(n, dx, incx, dy, incy)
DSWAP