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
solsy.f
Go to the documentation of this file.
1  SUBROUTINE solsy (WM, IWM, X, TEM)
2 CLLL. OPTIMIZE
3  INTEGER iwm
4  INTEGER iownd, iowns,
5  1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
6  2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
7  INTEGER i, meband, ml, mu
8  DOUBLE PRECISION wm, x, tem
9  DOUBLE PRECISION rowns,
10  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
11  DOUBLE PRECISION di, hl0, phl0, r
12  dimension wm(*), iwm(*), x(*), tem(*)
13  COMMON /ls0001/ rowns(209),
14  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
15  3 iownd(14), iowns(6),
16  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
17  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
18 C-----------------------------------------------------------------------
19 C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM
20 C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0.
21 C IF MITER IS 1 OR 2, IT CALLS DGETRS TO ACCOMPLISH THIS.
22 C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
23 C MATRIX, AND THEN COMPUTES THE SOLUTION.
24 C IF MITER IS 4 OR 5, IT CALLS DGBTRS.
25 C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES..
26 C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF
27 C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE.
28 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
29 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
30 C WM(1) = SQRT(UROUND) (NOT USED HERE),
31 C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3.
32 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
33 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
34 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
35 C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR
36 C ON OUTPUT, OF LENGTH N.
37 C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION.
38 C IERSL = OUTPUT FLAG (IN COMMON). IERSL = 0 IF NO TROUBLE OCCURRED.
39 C IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
40 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
41 C-----------------------------------------------------------------------
42  iersl = 0
43  go to(100, 100, 300, 400, 400), miter
44  100 CALL dgetrs( 'N', n, 1, wm(3), n, iwm(21), x, n, inlpck)
45  RETURN
46 C
47  300 phl0 = wm(2)
48  hl0 = h*el0
49  wm(2) = hl0
50  IF (hl0 .EQ. phl0) go to 330
51  r = hl0/phl0
52  DO 320 i = 1,n
53  di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2))
54  IF (dabs(di) .EQ. 0.0d0) go to 390
55  320 wm(i+2) = 1.0d0/di
56  330 DO 340 i = 1,n
57  340 x(i) = wm(i+2)*x(i)
58  RETURN
59  390 iersl = 1
60  RETURN
61 C
62  400 ml = iwm(1)
63  mu = iwm(2)
64  meband = 2*ml + mu + 1
65  CALL dgbtrs( 'N', n, ml, mu, 1, wm(3), meband, iwm(21), x, n,
66  * inlpck)
67  RETURN
68 C----------------------- END OF SUBROUTINE SOLSY -----------------------
69  END