← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/integ/solsy.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/integ/solsy.f
 
1
C/MEMBR ADD NAME=SOLSY,SSI=0
 
2
      subroutine solsy (wm, iwm, x, tem)
 
3
clll. optimize
 
4
      integer iwm
 
5
      integer iownd, iowns,
 
6
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
 
7
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
 
8
      integer i, meband, ml, mu
 
9
      double precision wm, x, tem
 
10
      double precision rownd, rowns,
 
11
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
 
12
      double precision di, hl0, phl0, r
 
13
      dimension wm(*), iwm(*), x(*), tem(*)
 
14
      common /ls0001/ rownd, rowns(209),
 
15
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
 
16
     3   iownd(14), iowns(6),
 
17
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
 
18
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
 
19
c-----------------------------------------------------------------------
 
20
c%purpose
 
21
c this routine manages the solution of the linear system arising from
 
22
c a chord iteration.  it is called if miter .ne. 0.
 
23
c if miter is 1 or 2, it calls dgesl to accomplish this.
 
24
c if miter = 3 it updates the coefficient h*el0 in the diagonal
 
25
c matrix, and then computes the solution.
 
26
c if miter is 4 or 5, it calls dgbsl.
 
27
c%calling sequence
 
28
c communication with solsy uses the following variables..
 
29
c wm    = real work space containing the inverse diagonal matrix if
 
30
c         miter = 3 and the lu decomposition of the matrix otherwise.
 
31
c         storage of matrix elements starts at wm(3).
 
32
c         wm also contains the following matrix-related data..
 
33
c         wm(1) = sqrt(uround) (not used here),
 
34
c         wm(2) = hl0, the previous value of h*el0, used if miter = 3.
 
35
c iwm   = integer work space containing pivot information, starting at
 
36
c         iwm(21), if miter is 1, 2, 4, or 5.  iwm also contains band
 
37
c         parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
 
38
c x     = the right-hand side vector on input, and the solution vector
 
39
c         on output, of length n.
 
40
c tem   = vector of work space of length n, not used in this version.
 
41
c iersl = output flag (in common).  iersl = 0 if no trouble occurred.
 
42
c         iersl = 1 if a singular matrix arose with miter = 3.
 
43
c this routine also uses the common variables el0, h, miter, and n.
 
44
c!
 
45
c-----------------------------------------------------------------------
 
46
      save/ls0001/
 
47
      iersl = 0
 
48
      go to (100, 100, 300, 400, 400), miter
 
49
 100  call dgesl (wm(3), n, n, iwm(21), x, 0)
 
50
      return
 
51
c
 
52
 300  phl0 = wm(2)
 
53
      hl0 = h*el0
 
54
      wm(2) = hl0
 
55
      if (hl0 .eq. phl0) go to 330
 
56
      r = hl0/phl0
 
57
      do 320 i = 1,n
 
58
        di = 1.0d+0 - r*(1.0d+0 - 1.0d+0/wm(i+2))
 
59
        if (abs(di) .eq. 0.0d+0) go to 390
 
60
 320    wm(i+2) = 1.0d+0/di
 
61
 330  do 340 i = 1,n
 
62
 340    x(i) = wm(i+2)*x(i)
 
63
      return
 
64
 390  iersl = 1
 
65
      return
 
66
c
 
67
 400  ml = iwm(1)
 
68
      mu = iwm(2)
 
69
      meband = 2*ml + mu + 1
 
70
      call dgbsl (wm(3), meband, n, ml, mu, iwm(21), x, 0)
 
71
      return
 
72
c----------------------- end of subroutine solsy -----------------------
 
73
      end