GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
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 INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH,
49  1 ialth, ipup, lmax, meo, nqnyh, nslp,
50  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
51  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
52  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
53  REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
54  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
55  COMMON /sls001/ conit, crate, el(13), elco(13,12),
56  1 hold, rmax, tesco(3,12),
57  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
58  2 init, mxstep, mxhnil, nhnil, nslast, nyh,
59  3 ialth, ipup, lmax, meo, nqnyh, nslp,
60  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
61  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
62  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
63  INTEGER I, MEBAND, ML, MU
64  REAL DI, HL0, PHL0, R
65 C
66 C***FIRST EXECUTABLE STATEMENT SSOLSY
67  iersl = 0
68  go to(100, 100, 300, 400, 400), miter
69  100 CALL sgetrs( 'N', n, 1, wm(3), n, iwm(21), x, n, inlpck)
70  RETURN
71 C
72  300 phl0 = wm(2)
73  hl0 = h*el0
74  wm(2) = hl0
75  IF (hl0 .EQ. phl0) go to 330
76  r = hl0/phl0
77  DO 320 i = 1,n
78  di = 1.0e0 - r*(1.0e0 - 1.0e0/wm(i+2))
79  IF (abs(di) .EQ. 0.0e0) go to 390
80  320 wm(i+2) = 1.0e0/di
81  330 DO 340 i = 1,n
82  340 x(i) = wm(i+2)*x(i)
83  RETURN
84  390 iersl = 1
85  RETURN
86 C
87  400 ml = iwm(1)
88  mu = iwm(2)
89  meband = 2*ml + mu + 1
90  CALL sgbtrs( 'N', n, ml, mu, 1, wm(3), meband, iwm(21), x, n,
91  * inlpck)
92  RETURN
93 C----------------------- END OF SUBROUTINE SSOLSY ----------------------
94  END
subroutine ssolsy(WM, IWM, X, TEM)
Definition: ssolsy.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to