![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to Band.out CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 41BSD / cmd / efl / efltest|
Return to Band.out CVS log | Up to [CSRG BSD Unix] / 41BSD / cmd / efl / efltest |
1.1 root 1: subroutine bnds(n, ml, m, g, nb, b)
2: integer m, n, nb
3: integer ml
4: real g(n, m), b(n, nb)
5: common /cstak/ ds
6: double precision ds(500)
7: integer istkgt, nerror, max0, iint, nerr, il
8: integer is(1000)
9: real rs(1000), ws(500)
10: logical ls(1000)
11: complex cs(500)
12: equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
13: c to solve a*x = b, where a is a banded matrix, using gaussian
14: c elimination with partial pivoting.
15: c mnemonic - double precision band solution of a system of
16: c linear algebraic equations.
17: c input -
18: c n - the order of the system.
19: c ml - the number of nonzero elements of a on and below the diagonal.
20: c m - the total number of nonzero elements in each row of a.
21: c g - the matrix a, with g(i,j) = a(i,i+j-ml).
22: c nb - the number of right-hand-sides b.
23: c b - the right-hand-sides.
24: c output -
25: c g - has been clobbered.
26: c b - the solution vectors, x.
27: c scratch space allocated - n*( (ml-1)*mu + 1 ) words.
28: c error states -
29: c 1 - n.lt.1.
30: c 2 - ml.lt.1.
31: c 3 - ml.gt.m.
32: c 4 - nb.lt.1.
33: c 5 - singular matrix. (recoverable)
34: c check the input for errors.
35: if (n .lt. 1) call seterr(14h bnds - n.lt.1, 14, 1, 2)
36: if (ml .lt. 1) call seterr(15h bnds - ml.lt.1, 16, 2, 2)
37: if (ml .gt. m) call seterr(15h bnds - ml.gt.m, 15, 3, 2)
38: if (nb .lt. 1) call seterr(15h bnds - nb.lt.1, 15, 4, 2)
39: call enter(1)
40: il = istkgt(max0(n*(ml-1), 1), 3)
41: iint = istkgt(n, 2)
42: call bndlu(n, ml, m, g, ws(il), is(iint))
43: if (nerror(nerr) .ne. 0) goto 1
44: call bndfb(n, ml, m, ws(il), g, is(iint), nb, b)
45: goto 2
46: 1 call erroff
47: call seterr(23h bnds - singular matrix, 23, 5, 1)
48: 2 call leave
49: return
50: end
51: subroutine bndlu(n, ml, m, g, l, int)
52: integer m, n, ml
53: integer int(n)
54: real g(n, m), l(n, ml)
55: integer min0, i, j, k, m1, m2
56: integer ll
57: real abs, eps, x, norm, amax1, r1mach
58: logical sing
59: integer temp, temp1
60: c to obtain the lu decomposition of a banded matrix,
61: c using gaussian elimination with partial pivoting.
62: c mnemonic - double precision band lu decomposition.
63: c input -
64: c n - the order of the matrix.
65: c ml - the number of nonzero elements of a on and below the diagonal.
66: c m - the number of nonzero elements in each row of a.
67: c g - the matrix a, with g(i,j) = a(i,i+j-ml).
68: c output -
69: c l - the lower triangular banded factor of a.
70: c g - the upper triangular banded factor of a.
71: c int - the row pivoting used.
72: c scratch storage allocated - none.
73: c error states -
74: c 1 - n.lt.1.
75: c 2 - ml.lt.1.
76: c 3 - m.lt.ml.
77: c 4 - singular matrix. (recoverable)
78: c l(n,ml-1).
79: c check the input for errors.
80: if (n .lt. 1) call seterr(15h bndlu - n.lt.1, 15, 1, 2)
81: if (ml .lt. 1) call seterr(16h bndlu - ml.lt.1, 16, 2, 2)
82: if (m .lt. ml) call seterr(16h bndlu - m.lt.ml, 16, 3, 2)
83: c protect against an existing error state.
84: call entsrc(i, 0)
85: sing = .false.
86: eps = r1mach(4)
87: m1 = ml-1
88: m2 = m-ml
89: ll = m1
90: i = 1
91: goto 2
92: 1 i = i+1
93: 2 if (i .gt. min0(m1, n)) goto 5
94: c set to 0 those elements
95: c of g which are undefined.
96: temp = ml+1-i
97: do 3 j = temp, m
98: temp1 = j-ll
99: g(i, temp1) = g(i, j)
100: 3 continue
101: ll = ll-1
102: temp = m-ll
103: do 4 j = temp, m
104: g(i, j) = 0.0e0
105: 4 continue
106: goto 1
107: 5 i = 1
108: goto 7
109: 6 i = i+1
110: 7 if (i .gt. min0(m2, n)) goto 9
111: c zero out lower rhs wart.
112: temp = ml+i
113: do 8 j = temp, m
114: temp1 = n+1-i
115: g(temp1, j) = 0.0e0
116: 8 continue
117: goto 6
118: c get || a || sub infinity.
119: 9 norm = 0.0e0
120: do 11 i = 1, n
121: int(i) = i
122: x = 0.0e0
123: do 10 j = 1, m
124: x = x+abs(g(i, j))
125: 10 continue
126: norm = amax1(norm, x)
127: 11 continue
128: do 20 k = 1, n
129: x = g(k, 1)
130: i = k
131: ll = min0(m1+k, n)
132: if (k .ge. ll) goto 14
133: temp = k+1
134: do 13 j = temp, ll
135: c get the pivot row.
136: if (abs(g(j, 1)) .le. abs(x)) goto 12
137: x = g(j, 1)
138: i = j
139: 12 continue
140: 13 continue
141: 14 int(k) = i
142: if (x .ne. 0.0e0) goto 15
143: sing = .true.
144: g(k, 1) = norm*eps
145: 15 if (ml .eq. 1 .or. k .eq. n) goto 20
146: if (i .eq. k) goto 17
147: do 16 j = 1, m
148: c need to interchange the rows.
149: x = g(k, j)
150: g(k, j) = g(i, j)
151: g(i, j) = x
152: 16 continue
153: 17 if (k .ge. ll) goto 20
154: temp = k+1
155: do 19 i = temp, ll
156: x = g(i, 1)/g(k, 1)
157: temp1 = i-k
158: l(k, temp1) = x
159: do 18 j = 2, m
160: g(i, j-1) = g(i, j)-x*g(k, j)
161: 18 continue
162: g(i, m) = 0.0e0
163: 19 continue
164: 20 continue
165: if (sing) call seterr(24h bndlu - singular matrix, 24, 4, 1)
166: return
167: end
168: subroutine bndfb(n, ml, m, l, u, int, nb, b)
169: integer m, n, nb, ml
170: integer int(n)
171: real l(n, ml), u(n, m), b(n, nb)
172: integer nerror, nerr
173: c to solve l*u*x = b, where l and u result from a call to bnds.
174: c mnemonic - double precision band forward elimination and
175: c back-solve.
176: c input -
177: c n - the order of the system.
178: c ml - the number of nonzero entries of l on and below
179: c the diagonal.
180: c m - the number of nonzero elements of u on and above
181: c the diagonal.
182: c l - the lower triangular banded factor.
183: c u - the upper triangular banded factor.
184: c int - the ordering of the rows of the system, due to pivoting.
185: c nb - the number of right-hand-sides.
186: c b - the right-hand-sides.
187: c output -
188: c b - the solution vectors.
189: c scratch space allocated - none.
190: c error states -
191: c 1 - n.lt.1.
192: c 2 - ml.lt.1.
193: c 3 - m.lt.ml.
194: c 4 - nb.lt.1.
195: c l(n,ml-1).
196: c check the input for errors.
197: if (n .lt. 1) call seterr(15h bndfb - n.lt.1, 15, 1, 2)
198: if (ml .lt. 1) call seterr(16h bndfb - ml.lt.1, 16, 2, 2)
199: if (m .lt. ml) call seterr(16h bndfb - m.lt.ml, 16, 3, 2)
200: if (nb .lt. 1) call seterr(16h bndfb - nb.lt.1, 16, 4, 2)
201: c protect against an existing error state.
202: call entsrc(nerr, 0)
203: call bndfe(n, ml, l, int, nb, b)
204: c do the forward-elimination.
205: call bndbs(n, m, u, nb, b)
206: c do the back-substitution.
207: return
208: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.