41BSD/cmd/efl/efltest/Band.out - view

File: [CSRG BSD Unix] / 41BSD / cmd / efl / efltest / Band.out
Revision 1.1: download - view: text, annotated - select for diffs
Tue Apr 24 16:12:53 2018 UTC (8 years, 1 month ago) by root
CVS tags: MAIN, HEAD

Initial revision

subroutine bnds(n, ml, m, g, nb, b) integer m, n, nb integer ml real g(n, m), b(n, nb) common /cstak/ ds double precision ds(500) integer istkgt, nerror, max0, iint, nerr, il integer is(1000) real rs(1000), ws(500) logical ls(1000) complex cs(500) equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1)) c to solve a*x = b, where a is a banded matrix, using gaussian c elimination with partial pivoting. c mnemonic - double precision band solution of a system of c linear algebraic equations. c input - c n - the order of the system. c ml - the number of nonzero elements of a on and below the diagonal. c m - the total number of nonzero elements in each row of a. c g - the matrix a, with g(i,j) = a(i,i+j-ml). c nb - the number of right-hand-sides b. c b - the right-hand-sides. c output - c g - has been clobbered. c b - the solution vectors, x. c scratch space allocated - n*( (ml-1)*mu + 1 ) words. c error states - c 1 - n.lt.1. c 2 - ml.lt.1. c 3 - ml.gt.m. c 4 - nb.lt.1. c 5 - singular matrix. (recoverable) c check the input for errors. if (n .lt. 1) call seterr(14h bnds - n.lt.1, 14, 1, 2) if (ml .lt. 1) call seterr(15h bnds - ml.lt.1, 16, 2, 2) if (ml .gt. m) call seterr(15h bnds - ml.gt.m, 15, 3, 2) if (nb .lt. 1) call seterr(15h bnds - nb.lt.1, 15, 4, 2) call enter(1) il = istkgt(max0(n*(ml-1), 1), 3) iint = istkgt(n, 2) call bndlu(n, ml, m, g, ws(il), is(iint)) if (nerror(nerr) .ne. 0) goto 1 call bndfb(n, ml, m, ws(il), g, is(iint), nb, b) goto 2 1 call erroff call seterr(23h bnds - singular matrix, 23, 5, 1) 2 call leave return end subroutine bndlu(n, ml, m, g, l, int) integer m, n, ml integer int(n) real g(n, m), l(n, ml) integer min0, i, j, k, m1, m2 integer ll real abs, eps, x, norm, amax1, r1mach logical sing integer temp, temp1 c to obtain the lu decomposition of a banded matrix, c using gaussian elimination with partial pivoting. c mnemonic - double precision band lu decomposition. c input - c n - the order of the matrix. c ml - the number of nonzero elements of a on and below the diagonal. c m - the number of nonzero elements in each row of a. c g - the matrix a, with g(i,j) = a(i,i+j-ml). c output - c l - the lower triangular banded factor of a. c g - the upper triangular banded factor of a. c int - the row pivoting used. c scratch storage allocated - none. c error states - c 1 - n.lt.1. c 2 - ml.lt.1. c 3 - m.lt.ml. c 4 - singular matrix. (recoverable) c l(n,ml-1). c check the input for errors. if (n .lt. 1) call seterr(15h bndlu - n.lt.1, 15, 1, 2) if (ml .lt. 1) call seterr(16h bndlu - ml.lt.1, 16, 2, 2) if (m .lt. ml) call seterr(16h bndlu - m.lt.ml, 16, 3, 2) c protect against an existing error state. call entsrc(i, 0) sing = .false. eps = r1mach(4) m1 = ml-1 m2 = m-ml ll = m1 i = 1 goto 2 1 i = i+1 2 if (i .gt. min0(m1, n)) goto 5 c set to 0 those elements c of g which are undefined. temp = ml+1-i do 3 j = temp, m temp1 = j-ll g(i, temp1) = g(i, j) 3 continue ll = ll-1 temp = m-ll do 4 j = temp, m g(i, j) = 0.0e0 4 continue goto 1 5 i = 1 goto 7 6 i = i+1 7 if (i .gt. min0(m2, n)) goto 9 c zero out lower rhs wart. temp = ml+i do 8 j = temp, m temp1 = n+1-i g(temp1, j) = 0.0e0 8 continue goto 6 c get || a || sub infinity. 9 norm = 0.0e0 do 11 i = 1, n int(i) = i x = 0.0e0 do 10 j = 1, m x = x+abs(g(i, j)) 10 continue norm = amax1(norm, x) 11 continue do 20 k = 1, n x = g(k, 1) i = k ll = min0(m1+k, n) if (k .ge. ll) goto 14 temp = k+1 do 13 j = temp, ll c get the pivot row. if (abs(g(j, 1)) .le. abs(x)) goto 12 x = g(j, 1) i = j 12 continue 13 continue 14 int(k) = i if (x .ne. 0.0e0) goto 15 sing = .true. g(k, 1) = norm*eps 15 if (ml .eq. 1 .or. k .eq. n) goto 20 if (i .eq. k) goto 17 do 16 j = 1, m c need to interchange the rows. x = g(k, j) g(k, j) = g(i, j) g(i, j) = x 16 continue 17 if (k .ge. ll) goto 20 temp = k+1 do 19 i = temp, ll x = g(i, 1)/g(k, 1) temp1 = i-k l(k, temp1) = x do 18 j = 2, m g(i, j-1) = g(i, j)-x*g(k, j) 18 continue g(i, m) = 0.0e0 19 continue 20 continue if (sing) call seterr(24h bndlu - singular matrix, 24, 4, 1) return end subroutine bndfb(n, ml, m, l, u, int, nb, b) integer m, n, nb, ml integer int(n) real l(n, ml), u(n, m), b(n, nb) integer nerror, nerr c to solve l*u*x = b, where l and u result from a call to bnds. c mnemonic - double precision band forward elimination and c back-solve. c input - c n - the order of the system. c ml - the number of nonzero entries of l on and below c the diagonal. c m - the number of nonzero elements of u on and above c the diagonal. c l - the lower triangular banded factor. c u - the upper triangular banded factor. c int - the ordering of the rows of the system, due to pivoting. c nb - the number of right-hand-sides. c b - the right-hand-sides. c output - c b - the solution vectors. c scratch space allocated - none. c error states - c 1 - n.lt.1. c 2 - ml.lt.1. c 3 - m.lt.ml. c 4 - nb.lt.1. c l(n,ml-1). c check the input for errors. if (n .lt. 1) call seterr(15h bndfb - n.lt.1, 15, 1, 2) if (ml .lt. 1) call seterr(16h bndfb - ml.lt.1, 16, 2, 2) if (m .lt. ml) call seterr(16h bndfb - m.lt.ml, 16, 3, 2) if (nb .lt. 1) call seterr(16h bndfb - nb.lt.1, 16, 4, 2) c protect against an existing error state. call entsrc(nerr, 0) call bndfe(n, ml, l, int, nb, b) c do the forward-elimination. call bndbs(n, m, u, nb, b) c do the back-substitution. return end

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.