subroutine solsy (wm, iwm, x, tem) clll. optimize integer iwm integer iownd, iowns, 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu integer i, meband, ml, mu double precision wm, x, tem double precision rownd, rowns, 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround double precision di, hl0, phl0, r dimension wm(1), iwm(1), x(1), tem(1) common /ls0001/ rownd, rowns(209), 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, 3 iownd(14), iowns(6), 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu c----------------------------------------------------------------------- c this routine manages the solution of the linear system arising from c a chord iteration. it is called if miter .ne. 0. c if miter is 1 or 2, it calls dgesl to accomplish this. c if miter = 3 it updates the coefficient h*el0 in the diagonal c matrix, and then computes the solution. c if miter is 4 or 5, it calls dgbsl. c communication with solsy uses the following variables.. c wm = real work space containing the inverse diagonal matrix if c miter = 3 and the lu decomposition of the matrix otherwise. c storage of matrix elements starts at wm(3). c wm also contains the following matrix-related data.. c wm(1) = sqrt(uround) (not used here), c wm(2) = hl0, the previous value of h*el0, used if miter = 3. c iwm = integer work space containing pivot information, starting at c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5. c x = the right-hand side vector on input, and the solution vector c on output, of length n. c tem = vector of work space of length n, not used in this version. c iersl = output flag (in common). iersl = 0 if no trouble occurred. c iersl = 1 if a singular matrix arose with miter = 3. c this routine also uses the common variables el0, h, miter, and n. c----------------------------------------------------------------------- iersl = 0 go to (100, 100, 300, 400, 400), miter 100 call dgesl (wm(3), n, n, iwm(21), x, 0) return c 300 phl0 = wm(2) hl0 = h*el0 wm(2) = hl0 if (hl0 .eq. phl0) go to 330 r = hl0/phl0 do 320 i = 1,n di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2)) if (dabs(di) .eq. 0.0d0) go to 390 320 wm(i+2) = 1.0d0/di 330 do 340 i = 1,n 340 x(i) = wm(i+2)*x(i) return 390 iersl = 1 return c 400 ml = iwm(1) mu = iwm(2) meband = 2*ml + mu + 1 call dgbsl (wm(3), meband, n, ml, mu, iwm(21), x, 0) return c----------------------- end of subroutine solsy ----------------------- end