![[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.