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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
|
subroutine dpofa(a,lda,n,info)
integer lda,n,info
double precision a(lda,*)
c
c dpofa factors a double precision symmetric positive definite
c matrix.
c
c dpofa is usually called by dpoco, but it can be called
c directly with a saving in time if rcond is not needed.
c (time for dpoco) = (1 + 18/n)*(time for dpofa) .
c
c on entry
c
c a double precision(lda, n)
c the symmetric matrix to be factored. only the
c diagonal and upper triangle are used.
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 on return
c
c a an upper triangular matrix r so that a = trans(r)*r
c where trans(r) is the transpose.
c the strict lower triangle is unaltered.
c if info .ne. 0 , the factorization is not complete.
c
c info integer
c = 0 for normal return.
c = k signals an error condition. the leading minor
c of order k is not positive definite.
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 ddot
c fortran sqrt
c
c internal variables
c
double precision ddot,t
double precision s
integer j,jm1,k
c begin block with ...exits to 40
c
c
do 30 j = 1, n
info = j
s = 0.0d0
jm1 = j - 1
if (jm1 .lt. 1) go to 20
do 10 k = 1, jm1
t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1)
t = t/a(k,k)
a(k,j) = t
s = s + t*t
10 continue
20 continue
s = a(j,j) - s
c ......exit
if (s .le. 0.0d0) go to 40
a(j,j) = sqrt(s)
30 continue
info = 0
40 continue
return
end
c====================== The end of dpofa ===============================
subroutine dtrsl(t,ldt,n,b,job,info)
integer ldt,n,job,info
double precision t(ldt,*),b(*)
c
c
c dtrsl solves systems of the form
c
c t * x = b
c or
c trans(t) * x = b
c
c where t is a triangular matrix of order n. here trans(t)
c denotes the transpose of the matrix t.
c
c on entry
c
c t double precision(ldt,n)
c t contains the matrix of the system. the zero
c elements of the matrix are not referenced, and
c the corresponding elements of the array can be
c used to store other information.
c
c ldt integer
c ldt is the leading dimension of the array t.
c
c n integer
c n is the order of the system.
c
c b double precision(n).
c b contains the right hand side of the system.
c
c job integer
c job specifies what kind of system is to be solved.
c if job is
c
c 00 solve t*x=b, t lower triangular,
c 01 solve t*x=b, t upper triangular,
c 10 solve trans(t)*x=b, t lower triangular,
c 11 solve trans(t)*x=b, t upper triangular.
c
c on return
c
c b b contains the solution, if info .eq. 0.
c otherwise b is unaltered.
c
c info integer
c info contains zero if the system is nonsingular.
c otherwise info contains the index of
c the first zero diagonal element of t.
c
c linpack. this version dated 08/14/78 .
c g. w. stewart, university of maryland, argonne national lab.
c
c subroutines and functions
c
c blas daxpy,ddot
c fortran mod
c
c internal variables
c
double precision ddot,temp
integer case,j,jj
c
c begin block permitting ...exits to 150
c
c check for zero diagonal elements.
c
do 10 info = 1, n
c ......exit
if (t(info,info) .eq. 0.0d0) go to 150
10 continue
info = 0
c
c determine the task and go to it.
c
case = 1
if (mod(job,10) .ne. 0) case = 2
if (mod(job,100)/10 .ne. 0) case = case + 2
go to (20,50,80,110), case
c
c solve t*x=b for t lower triangular
c
20 continue
b(1) = b(1)/t(1,1)
if (n .lt. 2) go to 40
do 30 j = 2, n
temp = -b(j-1)
call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1)
b(j) = b(j)/t(j,j)
30 continue
40 continue
go to 140
c
c solve t*x=b for t upper triangular.
c
50 continue
b(n) = b(n)/t(n,n)
if (n .lt. 2) go to 70
do 60 jj = 2, n
j = n - jj + 1
temp = -b(j+1)
call daxpy(j,temp,t(1,j+1),1,b(1),1)
b(j) = b(j)/t(j,j)
60 continue
70 continue
go to 140
c
c solve trans(t)*x=b for t lower triangular.
c
80 continue
b(n) = b(n)/t(n,n)
if (n .lt. 2) go to 100
do 90 jj = 2, n
j = n - jj + 1
b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1)
b(j) = b(j)/t(j,j)
90 continue
100 continue
go to 140
c
c solve trans(t)*x=b for t upper triangular.
c
110 continue
b(1) = b(1)/t(1,1)
if (n .lt. 2) go to 130
do 120 j = 2, n
b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1)
b(j) = b(j)/t(j,j)
120 continue
130 continue
140 continue
150 continue
return
end
c====================== The end of dtrsl ===============================
|
