1 |
C dgbsl.f |
2 |
C is freely available from netlib. It is not subject to any GNU License |
3 |
C set by the authors of the ASCEND math programming system. |
4 |
C $Date: 1996/04/30 18:17:11 $ $Revision: 1.1.1.1 $ |
5 |
C |
6 |
subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job) |
7 |
integer lda,n,ml,mu,ipvt(1),job |
8 |
double precision abd(lda,1),b(1) |
9 |
c |
10 |
c dgbsl solves the double precision band system |
11 |
c a * x = b or trans(a) * x = b |
12 |
c using the factors computed by dgbco or dgbfa. |
13 |
c |
14 |
c on entry |
15 |
c |
16 |
c abd double precision(lda, n) |
17 |
c the output from dgbco or dgbfa. |
18 |
c |
19 |
c lda integer |
20 |
c the leading dimension of the array abd . |
21 |
c |
22 |
c n integer |
23 |
c the order of the original matrix. |
24 |
c |
25 |
c ml integer |
26 |
c number of diagonals below the main diagonal. |
27 |
c |
28 |
c mu integer |
29 |
c number of diagonals above the main diagonal. |
30 |
c |
31 |
c ipvt integer(n) |
32 |
c the pivot vector from dgbco or dgbfa. |
33 |
c |
34 |
c b double precision(n) |
35 |
c the right hand side vector. |
36 |
c |
37 |
c job integer |
38 |
c = 0 to solve a*x = b , |
39 |
c = nonzero to solve trans(a)*x = b , where |
40 |
c trans(a) is the transpose. |
41 |
c |
42 |
c on return |
43 |
c |
44 |
c b the solution vector x . |
45 |
c |
46 |
c error condition |
47 |
c |
48 |
c a division by zero will occur if the input factor contains a |
49 |
c zero on the diagonal. technically this indicates singularity |
50 |
c but it is often caused by improper arguments or improper |
51 |
c setting of lda . it will not occur if the subroutines are |
52 |
c called correctly and if dgbco has set rcond .gt. 0.0 |
53 |
c or dgbfa has set info .eq. 0 . |
54 |
c |
55 |
c to compute inverse(a) * c where c is a matrix |
56 |
c with p columns |
57 |
c call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z) |
58 |
c if (rcond is too small) go to ... |
59 |
c do 10 j = 1, p |
60 |
c call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) |
61 |
c 10 continue |
62 |
c |
63 |
c linpack. this version dated 08/14/78 . |
64 |
c cleve moler, university of new mexico, argonne national lab. |
65 |
c |
66 |
c subroutines and functions |
67 |
c |
68 |
c blas daxpy,ddot |
69 |
c fortran min0 |
70 |
c |
71 |
c internal variables |
72 |
c |
73 |
double precision ddot,t |
74 |
integer k,kb,l,la,lb,lm,m,nm1 |
75 |
c |
76 |
m = mu + ml + 1 |
77 |
nm1 = n - 1 |
78 |
if (job .ne. 0) go to 50 |
79 |
c |
80 |
c job = 0 , solve a * x = b |
81 |
c first solve l*y = b |
82 |
c |
83 |
if (ml .eq. 0) go to 30 |
84 |
if (nm1 .lt. 1) go to 30 |
85 |
do 20 k = 1, nm1 |
86 |
lm = min0(ml,n-k) |
87 |
l = ipvt(k) |
88 |
t = b(l) |
89 |
if (l .eq. k) go to 10 |
90 |
b(l) = b(k) |
91 |
b(k) = t |
92 |
10 continue |
93 |
call daxpy(lm,t,abd(m+1,k),1,b(k+1),1) |
94 |
20 continue |
95 |
30 continue |
96 |
c |
97 |
c now solve u*x = y |
98 |
c |
99 |
do 40 kb = 1, n |
100 |
k = n + 1 - kb |
101 |
b(k) = b(k)/abd(m,k) |
102 |
lm = min0(k,m) - 1 |
103 |
la = m - lm |
104 |
lb = k - lm |
105 |
t = -b(k) |
106 |
call daxpy(lm,t,abd(la,k),1,b(lb),1) |
107 |
40 continue |
108 |
go to 100 |
109 |
50 continue |
110 |
c |
111 |
c job = nonzero, solve trans(a) * x = b |
112 |
c first solve trans(u)*y = b |
113 |
c |
114 |
do 60 k = 1, n |
115 |
lm = min0(k,m) - 1 |
116 |
la = m - lm |
117 |
lb = k - lm |
118 |
t = ddot(lm,abd(la,k),1,b(lb),1) |
119 |
b(k) = (b(k) - t)/abd(m,k) |
120 |
60 continue |
121 |
c |
122 |
c now solve trans(l)*x = y |
123 |
c |
124 |
if (ml .eq. 0) go to 90 |
125 |
if (nm1 .lt. 1) go to 90 |
126 |
do 80 kb = 1, nm1 |
127 |
k = n - kb |
128 |
lm = min0(ml,n-k) |
129 |
b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1) |
130 |
l = ipvt(k) |
131 |
if (l .eq. k) go to 70 |
132 |
t = b(l) |
133 |
b(l) = b(k) |
134 |
b(k) = t |
135 |
70 continue |
136 |
80 continue |
137 |
90 continue |
138 |
100 continue |
139 |
return |
140 |
end |