-
Notifications
You must be signed in to change notification settings - Fork 315
/
dtrsl.f
141 lines (141 loc) · 3.64 KB
/
dtrsl.f
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
c Triangular Solve dtrsl()
c ----------------
c 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
subroutine dtrsl(t,ldt,n,b,job,info)
integer ldt,n,job,info
double precision t(ldt,*),b(*)
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
if (t(info,info) .eq. 0.0d0) go to 150
c ......exit
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 Case 1 (job = 00):
c solve t*x=b for t lower triangular
c
20 continue
b(1) = b(1)/t(1,1)
if (n .ge. 2) then
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
endif
go to 140
c
C Case 2 (job = 01):
c solve t*x=b for t upper triangular.
c
50 continue
b(n) = b(n)/t(n,n)
if (n .ge. 2) then
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
endif
go to 140
c
C Case 3 (job = 10):
c solve trans(t)*x=b for t lower triangular.
c
80 continue
b(n) = b(n)/t(n,n)
if (n .ge. 2) then
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
endif
go to 140
c
C Case 4 (job = 11):
c solve trans(t)*x=b for t upper triangular.
c
110 continue
b(1) = b(1)/t(1,1)
if (n .ge. 2) then
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
endif
C
140 continue
c EXIT:
150 continue
return
end