|
|
subroutine sppdi(ap,n,det,job)
|
|
|
integer n,job
|
|
|
real ap(1)
|
|
|
real det(2)
|
|
|
c
|
|
|
c sppdi computes the determinant and inverse
|
|
|
c of a real symmetric positive definite matrix
|
|
|
c using the factors computed by sppco or sppfa .
|
|
|
c
|
|
|
c on entry
|
|
|
c
|
|
|
c ap real (n*(n+1)/2)
|
|
|
c the output from sppco or sppfa.
|
|
|
c
|
|
|
c n integer
|
|
|
c the order of the matrix a .
|
|
|
c
|
|
|
c job integer
|
|
|
c = 11 both determinant and inverse.
|
|
|
c = 01 inverse only.
|
|
|
c = 10 determinant only.
|
|
|
c
|
|
|
c on return
|
|
|
c
|
|
|
c ap the upper triangular half of the inverse .
|
|
|
c the strict lower triangle is unaltered.
|
|
|
c
|
|
|
c det real(2)
|
|
|
c determinant of original matrix if requested.
|
|
|
c otherwise not referenced.
|
|
|
c determinant = det(1) * 10.0**det(2)
|
|
|
c with 1.0 .le. det(1) .lt. 10.0
|
|
|
c or det(1) .eq. 0.0 .
|
|
|
c
|
|
|
c error condition
|
|
|
c
|
|
|
c a division by zero will occur if the input factor contains
|
|
|
c a zero on the diagonal and the inverse is requested.
|
|
|
c it will not occur if the subroutines are called correctly
|
|
|
c and if spoco or spofa has set info .eq. 0 .
|
|
|
c
|
|
|
c linpack. this version dated 08/14/78 .
|
|
|
c cleve moler, university of new mexico, argonne national lab.
|
|
|
c
|
|
|
c subroutines and functions
|
|
|
c
|
|
|
c blas saxpy,sscal
|
|
|
c fortran mod
|
|
|
c
|
|
|
c internal variables
|
|
|
c
|
|
|
real t
|
|
|
real s
|
|
|
integer i,ii,j,jj,jm1,j1,k,kj,kk,kp1,k1
|
|
|
c
|
|
|
c compute determinant
|
|
|
c
|
|
|
if (job/10 .eq. 0) go to 70
|
|
|
det(1) = 1.0e0
|
|
|
det(2) = 0.0e0
|
|
|
s = 10.0e0
|
|
|
ii = 0
|
|
|
do 50 i = 1, n
|
|
|
ii = ii + i
|
|
|
det(1) = ap(ii)**2*det(1)
|
|
|
c ...exit
|
|
|
if (det(1) .eq. 0.0e0) go to 60
|
|
|
10 if (det(1) .ge. 1.0e0) go to 20
|
|
|
det(1) = s*det(1)
|
|
|
det(2) = det(2) - 1.0e0
|
|
|
go to 10
|
|
|
20 continue
|
|
|
30 if (det(1) .lt. s) go to 40
|
|
|
det(1) = det(1)/s
|
|
|
det(2) = det(2) + 1.0e0
|
|
|
go to 30
|
|
|
40 continue
|
|
|
50 continue
|
|
|
60 continue
|
|
|
70 continue
|
|
|
c
|
|
|
c compute inverse(r)
|
|
|
c
|
|
|
if (mod(job,10) .eq. 0) go to 140
|
|
|
kk = 0
|
|
|
do 100 k = 1, n
|
|
|
k1 = kk + 1
|
|
|
kk = kk + k
|
|
|
ap(kk) = 1.0e0/ap(kk)
|
|
|
t = -ap(kk)
|
|
|
call sscal(k-1,t,ap(k1),1)
|
|
|
kp1 = k + 1
|
|
|
j1 = kk + 1
|
|
|
kj = kk + k
|
|
|
if (n .lt. kp1) go to 90
|
|
|
do 80 j = kp1, n
|
|
|
t = ap(kj)
|
|
|
ap(kj) = 0.0e0
|
|
|
call saxpy(k,t,ap(k1),1,ap(j1),1)
|
|
|
j1 = j1 + j
|
|
|
kj = kj + j
|
|
|
80 continue
|
|
|
90 continue
|
|
|
100 continue
|
|
|
c
|
|
|
c form inverse(r) * trans(inverse(r))
|
|
|
c
|
|
|
jj = 0
|
|
|
do 130 j = 1, n
|
|
|
j1 = jj + 1
|
|
|
jj = jj + j
|
|
|
jm1 = j - 1
|
|
|
k1 = 1
|
|
|
kj = j1
|
|
|
if (jm1 .lt. 1) go to 120
|
|
|
do 110 k = 1, jm1
|
|
|
t = ap(kj)
|
|
|
call saxpy(k,t,ap(j1),1,ap(k1),1)
|
|
|
k1 = k1 + k
|
|
|
kj = kj + 1
|
|
|
110 continue
|
|
|
120 continue
|
|
|
t = ap(jj)
|
|
|
call sscal(j,t,ap(j1),1)
|
|
|
130 continue
|
|
|
140 continue
|
|
|
return
|
|
|
end
|
|
|
|
