110 SUBROUTINE dsptri( UPLO, N, AP, IPIV, WORK, INFO )
123 DOUBLE PRECISION ap( * ), work( * )
129 DOUBLE PRECISION one, zero
130 parameter( one = 1.0d+0, zero = 0.0d+0 )
134 INTEGER j, k, kc, kcnext, kp, kpc, kstep, kx, npp
135 DOUBLE PRECISION ak, akkp1, akp1, d, t, temp
139 DOUBLE PRECISION ddot
153 upper =
lsame( uplo,
'U' )
154 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
156 ELSE IF( n.LT.0 )
THEN
160 CALL
xerbla(
'DSPTRI', -info )
176 DO 10 info = n, 1, -1
177 IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
187 IF( ipiv( info ).GT.0 .AND. ap( kp ).EQ.zero )
189 kp = kp + n - info + 1
211 IF( ipiv( k ).GT.0 )
THEN
217 ap( kc+k-1 ) = one / ap( kc+k-1 )
222 CALL
dcopy( k-1, ap( kc ), 1, work, 1 )
223 CALL
dspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
225 ap( kc+k-1 ) = ap( kc+k-1 ) -
226 $
ddot( k-1, work, 1, ap( kc ), 1 )
235 t = abs( ap( kcnext+k-1 ) )
236 ak = ap( kc+k-1 ) / t
237 akp1 = ap( kcnext+k ) / t
238 akkp1 = ap( kcnext+k-1 ) / t
239 d = t*( ak*akp1-one )
240 ap( kc+k-1 ) = akp1 / d
241 ap( kcnext+k ) = ak / d
242 ap( kcnext+k-1 ) = -akkp1 / d
247 CALL
dcopy( k-1, ap( kc ), 1, work, 1 )
248 CALL
dspmv( uplo, k-1, -one, ap, work, 1, zero, ap( kc ),
250 ap( kc+k-1 ) = ap( kc+k-1 ) -
251 $
ddot( k-1, work, 1, ap( kc ), 1 )
252 ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -
253 $
ddot( k-1, ap( kc ), 1, ap( kcnext ),
255 CALL
dcopy( k-1, ap( kcnext ), 1, work, 1 )
256 CALL
dspmv( uplo, k-1, -one, ap, work, 1, zero,
258 ap( kcnext+k ) = ap( kcnext+k ) -
259 $
ddot( k-1, work, 1, ap( kcnext ), 1 )
262 kcnext = kcnext + k + 1
265 kp = abs( ipiv( k ) )
271 kpc = ( kp-1 )*kp / 2 + 1
272 CALL
dswap( kp-1, ap( kc ), 1, ap( kpc ), 1 )
274 DO 40
j = kp + 1, k - 1
277 ap( kc+
j-1 ) = ap( kx )
281 ap( kc+k-1 ) = ap( kpc+kp-1 )
282 ap( kpc+kp-1 ) = temp
283 IF( kstep.EQ.2 )
THEN
284 temp = ap( kc+k+k-1 )
285 ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
286 ap( kc+k+kp-1 ) = temp
312 kcnext = kc - ( n-k+2 )
313 IF( ipiv( k ).GT.0 )
THEN
319 ap( kc ) = one / ap( kc )
324 CALL
dcopy( n-k, ap( kc+1 ), 1, work, 1 )
325 CALL
dspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1,
326 $ zero, ap( kc+1 ), 1 )
327 ap( kc ) = ap( kc ) -
ddot( n-k, work, 1, ap( kc+1 ), 1 )
336 t = abs( ap( kcnext+1 ) )
337 ak = ap( kcnext ) / t
339 akkp1 = ap( kcnext+1 ) / t
340 d = t*( ak*akp1-one )
341 ap( kcnext ) = akp1 / d
343 ap( kcnext+1 ) = -akkp1 / d
348 CALL
dcopy( n-k, ap( kc+1 ), 1, work, 1 )
349 CALL
dspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
350 $ zero, ap( kc+1 ), 1 )
351 ap( kc ) = ap( kc ) -
ddot( n-k, work, 1, ap( kc+1 ), 1 )
352 ap( kcnext+1 ) = ap( kcnext+1 ) -
353 $
ddot( n-k, ap( kc+1 ), 1,
354 $ ap( kcnext+2 ), 1 )
355 CALL
dcopy( n-k, ap( kcnext+2 ), 1, work, 1 )
356 CALL
dspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1,
357 $ zero, ap( kcnext+2 ), 1 )
358 ap( kcnext ) = ap( kcnext ) -
359 $
ddot( n-k, work, 1, ap( kcnext+2 ), 1 )
362 kcnext = kcnext - ( n-k+3 )
365 kp = abs( ipiv( k ) )
371 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2 + 1
373 $ CALL
dswap( n-kp, ap( kc+kp-k+1 ), 1, ap( kpc+1 ), 1 )
375 DO 70
j = k + 1, kp - 1
378 ap( kc+
j-k ) = ap( kx )
384 IF( kstep.EQ.2 )
THEN
385 temp = ap( kc-n+k-1 )
386 ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
387 ap( kc-n+kp-1 ) = temp
LOGICAL function lsame(CA, CB)
LSAME
subroutine dspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
DSPMV
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dsptri(UPLO, N, AP, IPIV, WORK, INFO)
DSPTRI
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP