1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
subroutine dpodi(a,lda,n,det,job)
integer lda,n,job
double precision a(lda,1)
double precision det(2)
c
c dpodi computes the determinant and inverse of a certain
c double precision symmetric positive definite matrix (see below)
c using the factors computed by dpoco, dpofa or dqrdc.
c
c on entry
c
c a double precision(lda, n)
c the output a from dpoco or dpofa
c or the output x from dqrdc.
c
c lda integer
c the leading dimension of the array a .
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 a if dpoco or dpofa was used to factor a then
c dpodi produces the upper half of inverse(a) .
c if dqrdc was used to decompose x then
c dpodi produces the upper half of inverse(trans(x)*x)
c where trans(x) is the transpose.
c elements of a below the diagonal are unchanged.
c if the units digit of job is zero, a is unchanged.
c
c det double precision(2)
c determinant of a or of trans(x)*x 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 dpoco or dpofa 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 daxpy,dscal
c fortran mod
c
c internal variables
c
double precision t
double precision s
integer i,j,jm1,k,kp1
c
c compute determinant
c
if (job/10 .eq. 0) go to 70
det(1) = 1.0d0
det(2) = 0.0d0
s = 10.0d0
do 50 i = 1, n
det(1) = a(i,i)**2*det(1)
c ...exit
if (det(1) .eq. 0.0d0) go to 60
10 if (det(1) .ge. 1.0d0) go to 20
det(1) = s*det(1)
det(2) = det(2) - 1.0d0
go to 10
20 continue
30 if (det(1) .lt. s) go to 40
det(1) = det(1)/s
det(2) = det(2) + 1.0d0
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
do 100 k = 1, n
a(k,k) = 1.0d0/a(k,k)
t = -a(k,k)
call dscal(k-1,t,a(1,k),1)
kp1 = k + 1
if (n .lt. kp1) go to 90
do 80 j = kp1, n
t = a(k,j)
a(k,j) = 0.0d0
call daxpy(k,t,a(1,k),1,a(1,j),1)
80 continue
90 continue
100 continue
c
c form inverse(r) * trans(inverse(r))
c
do 130 j = 1, n
jm1 = j - 1
if (jm1 .lt. 1) go to 120
do 110 k = 1, jm1
t = a(k,j)
call daxpy(k,t,a(1,j),1,a(1,k),1)
110 continue
120 continue
t = a(j,j)
call dscal(j,t,a(1,j),1)
130 continue
140 continue
return
end
|
