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