SUBROUTINE DSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
CHARACTER UPLO
INTEGER INFO, N
INTEGER IPIV( * )
DOUBLE PRECISION AP( * ), WORK( * )
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
LOGICAL UPPER
INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
LOGICAL LSAME
DOUBLE PRECISION DDOT
EXTERNAL LSAME, DDOT
EXTERNAL DCOPY, DSPMV, DSWAP, XERBLA
INTRINSIC ABS
INFO = 0
UPPER = LSAME( UPLO, 'U' )
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DSPTRI', -INFO )
RETURN
END IF
IF( N.EQ.0 )
$ RETURN
IF( UPPER ) THEN
KP = N*( N+1 ) / 2
DO 10 INFO = N, 1, -1
IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
$ RETURN
KP = KP - INFO
10 CONTINUE
ELSE
KP = 1
DO 20 INFO = 1, N
IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
$ RETURN
KP = KP + N - INFO + 1
20 CONTINUE
END IF
INFO = 0
IF( UPPER ) THEN
K = 1
KC = 1
30 CONTINUE
IF( K.GT.N )
$ GO TO 50
KCNEXT = KC + K
IF( IPIV( K ).GT.0 ) THEN
AP( KC+K-1 ) = ONE / AP( KC+K-1 )
IF( K.GT.1 ) THEN
CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
$ 1 )
AP( KC+K-1 ) = AP( KC+K-1 ) -
$ DDOT( K-1, WORK, 1, AP( KC ), 1 )
END IF
KSTEP = 1
ELSE
T = ABS( AP( KCNEXT+K-1 ) )
AK = AP( KC+K-1 ) / T
AKP1 = AP( KCNEXT+K ) / T
AKKP1 = AP( KCNEXT+K-1 ) / T
D = T*( AK*AKP1-ONE )
AP( KC+K-1 ) = AKP1 / D
AP( KCNEXT+K ) = AK / D
AP( KCNEXT+K-1 ) = -AKKP1 / D
IF( K.GT.1 ) THEN
CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
$ 1 )
AP( KC+K-1 ) = AP( KC+K-1 ) -
$ DDOT( K-1, WORK, 1, AP( KC ), 1 )
AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
$ DDOT( K-1, AP( KC ), 1, AP( KCNEXT ),
$ 1 )
CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
$ AP( KCNEXT ), 1 )
AP( KCNEXT+K ) = AP( KCNEXT+K ) -
$ DDOT( K-1, WORK, 1, AP( KCNEXT ), 1 )
END IF
KSTEP = 2
KCNEXT = KCNEXT + K + 1
END IF
KP = ABS( IPIV( K ) )
IF( KP.NE.K ) THEN
KPC = ( KP-1 )*KP / 2 + 1
CALL DSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
KX = KPC + KP - 1
DO 40 J = KP + 1, K - 1
KX = KX + J - 1
TEMP = AP( KC+J-1 )
AP( KC+J-1 ) = AP( KX )
AP( KX ) = TEMP
40 CONTINUE
TEMP = AP( KC+K-1 )
AP( KC+K-1 ) = AP( KPC+KP-1 )
AP( KPC+KP-1 ) = TEMP
IF( KSTEP.EQ.2 ) THEN
TEMP = AP( KC+K+K-1 )
AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
AP( KC+K+KP-1 ) = TEMP
END IF
END IF
K = K + KSTEP
KC = KCNEXT
GO TO 30
50 CONTINUE
ELSE
NPP = N*( N+1 ) / 2
K = N
KC = NPP
60 CONTINUE
IF( K.LT.1 )
$ GO TO 80
KCNEXT = KC - ( N-K+2 )
IF( IPIV( K ).GT.0 ) THEN
AP( KC ) = ONE / AP( KC )
IF( K.LT.N ) THEN
CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
$ ZERO, AP( KC+1 ), 1 )
AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
END IF
KSTEP = 1
ELSE
T = ABS( AP( KCNEXT+1 ) )
AK = AP( KCNEXT ) / T
AKP1 = AP( KC ) / T
AKKP1 = AP( KCNEXT+1 ) / T
D = T*( AK*AKP1-ONE )
AP( KCNEXT ) = AKP1 / D
AP( KC ) = AK / D
AP( KCNEXT+1 ) = -AKKP1 / D
IF( K.LT.N ) THEN
CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
$ ZERO, AP( KC+1 ), 1 )
AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
$ DDOT( N-K, AP( KC+1 ), 1,
$ AP( KCNEXT+2 ), 1 )
CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
$ ZERO, AP( KCNEXT+2 ), 1 )
AP( KCNEXT ) = AP( KCNEXT ) -
$ DDOT( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
END IF
KSTEP = 2
KCNEXT = KCNEXT - ( N-K+3 )
END IF
KP = ABS( IPIV( K ) )
IF( KP.NE.K ) THEN
KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
IF( KP.LT.N )
$ CALL DSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
KX = KC + KP - K
DO 70 J = K + 1, KP - 1
KX = KX + N - J + 1
TEMP = AP( KC+J-K )
AP( KC+J-K ) = AP( KX )
AP( KX ) = TEMP
70 CONTINUE
TEMP = AP( KC )
AP( KC ) = AP( KPC )
AP( KPC ) = TEMP
IF( KSTEP.EQ.2 ) THEN
TEMP = AP( KC-N+K-1 )
AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
AP( KC-N+KP-1 ) = TEMP
END IF
END IF
K = K - KSTEP
KC = KCNEXT
GO TO 60
80 CONTINUE
END IF
RETURN
END