GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ssolsy.f
Go to the documentation of this file.
1  SUBROUTINE ssolsy (WM, IWM, X, TEM)
2 C***BEGIN PROLOGUE SSOLSY
3 C***SUBSIDIARY
4 C***PURPOSE ODEPACK linear system solver.
5 C***TYPE SINGLE PRECISION (SSOLSY-S, DSOLSY-D)
6 C***AUTHOR Hindmarsh, Alan C., (LLNL)
7 C***DESCRIPTION
8 C
9 C This routine manages the solution of the linear system arising from
10 C a chord iteration. It is called if MITER .ne. 0.
11 C If MITER is 1 or 2, it calls SGETRF to accomplish this.
12 C If MITER = 3 it updates the coefficient h*EL0 in the diagonal
13 C matrix, and then computes the solution.
14 C If MITER is 4 or 5, it calls SGBTRS.
15 C Communication with SSOLSY uses the following variables:
16 C WM = real work space containing the inverse diagonal matrix if
17 C MITER = 3 and the LU decomposition of the matrix otherwise.
18 C Storage of matrix elements starts at WM(3).
19 C WM also contains the following matrix-related data:
20 C WM(1) = SQRT(UROUND) (not used here),
21 C WM(2) = HL0, the previous value of h*EL0, used if MITER = 3.
22 C IWM = integer work space containing pivot information, starting at
23 C IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
24 C parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
25 C X = the right-hand side vector on input, and the solution vector
26 C on output, of length N.
27 C TEM = vector of work space of length N, not used in this version.
28 C IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred.
29 C IERSL = 1 if a singular matrix arose with MITER = 3.
30 C This routine also uses the COMMON variables EL0, H, MITER, and N.
31 C
32 C***SEE ALSO SLSODE
33 C***ROUTINES CALLED SGBTRS, SGETRS
34 C***COMMON BLOCKS SLS001
35 C***REVISION HISTORY (YYMMDD)
36 C 791129 DATE WRITTEN
37 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
38 C 890503 Minor cosmetic changes. (FNF)
39 C 930809 Renamed to allow single/double precision versions. (ACH)
40 C 010412 Reduced size of Common block /SLS001/. (ACH)
41 C 031105 Restored 'own' variables to Common block /SLS001/, to
42 C enable interrupt/restart feature. (ACH)
43 C***END PROLOGUE SSOLSY
44 C**End
45  INTEGER iwm
46  REAL wm, x, tem
47  dimension wm(*), iwm(*), x(*), tem(*)
48  INTEGER iownd, iowns,
49  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
50  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
51  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
52  REAL rowns,
53  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
54  COMMON /sls001/ rowns(209),
55  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
56  2 iownd(6), iowns(6),
57  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
58  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
59  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
60  INTEGER i, meband, ml, mu
61  REAL di, hl0, phl0, r
62 C
63 C***FIRST EXECUTABLE STATEMENT SSOLSY
64  iersl = 0
65  go to(100, 100, 300, 400, 400), miter
66  100 CALL sgetrs( 'N', n, 1, wm(3), n, iwm(21), x, n, inlpck)
67  RETURN
68 C
69  300 phl0 = wm(2)
70  hl0 = h*el0
71  wm(2) = hl0
72  IF (hl0 .EQ. phl0) go to 330
73  r = hl0/phl0
74  DO 320 i = 1,n
75  di = 1.0e0 - r*(1.0e0 - 1.0e0/wm(i+2))
76  IF (abs(di) .EQ. 0.0e0) go to 390
77  320 wm(i+2) = 1.0e0/di
78  330 DO 340 i = 1,n
79  340 x(i) = wm(i+2)*x(i)
80  RETURN
81  390 iersl = 1
82  RETURN
83 C
84  400 ml = iwm(1)
85  mu = iwm(2)
86  meband = 2*ml + mu + 1
87  CALL sgbtrs( 'N', n, ml, mu, 1, wm(3), meband, iwm(21), x, n,
88  * inlpck)
89  RETURN
90 C----------------------- END OF SUBROUTINE SSOLSY ----------------------
91  END