GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
solsy.f
Go to the documentation of this file.
1  SUBROUTINE solsy (WM, IWM, X, TEM)
2 CLLL. OPTIMIZE
3  INTEGER IWM
4  INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
5  1 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, NYH,
6  2 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP
7  INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
8  2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
9  INTEGER I, MEBAND, ML, MU
10  DOUBLE PRECISION WM, X, TEM
11  DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
12  1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
13  DOUBLE PRECISION DI, HL0, PHL0, R
14  dimension wm(*), iwm(*), x(*), tem(*)
15  COMMON /ls0001/ conit, crate, el(13), elco(13,12),
16  1 hold, rmax, tesco(3,12),
17  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
18  2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
19  3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh,
20  3 ialth, ipup, lmax, meo, nqnyh, nslp,
21  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
22  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
23 C-----------------------------------------------------------------------
24 C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM
25 C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0.
26 C IF MITER IS 1 OR 2, IT CALLS DGETRS TO ACCOMPLISH THIS.
27 C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
28 C MATRIX, AND THEN COMPUTES THE SOLUTION.
29 C IF MITER IS 4 OR 5, IT CALLS DGBTRS.
30 C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES..
31 C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF
32 C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE.
33 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
34 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
35 C WM(1) = SQRT(UROUND) (NOT USED HERE),
36 C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3.
37 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
38 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
39 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
40 C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR
41 C ON OUTPUT, OF LENGTH N.
42 C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION.
43 C IERSL = OUTPUT FLAG (IN COMMON). IERSL = 0 IF NO TROUBLE OCCURRED.
44 C IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
45 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
46 C-----------------------------------------------------------------------
47  iersl = 0
48  GO TO (100, 100, 300, 400, 400), miter
49  100 CALL dgetrs ( 'N', n, 1, wm(3), n, iwm(21), x, n, inlpck)
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.0d0 - r*(1.0d0 - 1.0d0/wm(i+2))
59  IF (dabs(di) .EQ. 0.0d0) GO TO 390
60  320 wm(i+2) = 1.0d0/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 dgbtrs ( 'N', n, ml, mu, 1, wm(3), meband, iwm(21), x, n,
71  * inlpck)
72  RETURN
73 C----------------------- END OF SUBROUTINE SOLSY -----------------------
74  END
subroutine solsy(WM, IWM, X, TEM)
Definition: solsy.f:2