V3FIT
solsy.f
1  SUBROUTINE solsy(wm, iwm, x, tem)
2  USE stel_kinds
3  USE liprec, ONLY: li_gbsl, li_gesl
4  IMPLICIT NONE
5 C-----------------------------------------------
6 C D u m m y A r g u m e n t s
7 C-----------------------------------------------
8  INTEGER, DIMENSION(*) :: iwm
9  REAL(rprec), DIMENSION(*) :: wm, x, tem
10 C-----------------------------------------------
11 C C o m m o n B l o c k s
12 C-----------------------------------------------
13 C... /ls0001/
14  COMMON /ls0001/ rowns(209), ccmax, el0, h, hmin, hmxi, hu, rc, tn
15  1 , uround, iownd(14), iowns(6), icf, ierpj, iersl, jcur, jstart
16  2 , kflag, l, meth, miter, maxord, maxcor, msbp, mxncf, n, nq,
17  3 nst, nfe, nje, nqu
18  INTEGER iownd, iowns, icf, ierpj, iersl, jcur, jstart, kflag, l
19  1 , meth, miter, maxord, maxcor, msbp, mxncf, n, nq, nst, nfe,
20  2 nje, nqu
21  REAL(rprec) ::rowns, ccmax, el0, h, hmin, hmxi, hu, rc, tn,
22  1 uround
23 C-----------------------------------------------
24 C L o c a l V a r i a b l e s
25 C-----------------------------------------------
26  INTEGER :: i, meband, ml, mu
27  REAL(rprec) :: di, hl0, phl0, r
28 C-----------------------------------------------
29 clll. optimize
30 c-----------------------------------------------------------------------
31 c this routine manages the solution of the linear system arising from
32 c a chord iteration. it is called if miter .ne. 0.
33 c if miter is 1 or 2, it calls sgesl to accomplish this.
34 c if miter = 3 it updates the coefficient h*el0 in the diagonal
35 c matrix, and then computes the solution.
36 c if miter is 4 or 5, it calls sgbsl.
37 c communication with solsy uses the following variables..
38 c wm = real work space containing the inverse diagonal matrix if
39 c miter = 3 and the lu decomposition of the matrix otherwise.
40 c storage of matrix elements starts at wm(3).
41 c wm also contains the following matrix-related data..
42 c wm(1) = sqrt(uround) (not used here),
43 c wm(2) = hl0, the previous value of h*el0, used if miter = 3.
44 c iwm = integer work space containing pivot information, starting at
45 c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band
46 c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
47 c x = the right-hand side vector on input, and the solution vector
48 c on output, of length n.
49 c tem = vector of work space of length n, not used in this version.
50 c iersl = output flag (in common). iersl = 0 if no trouble occurred.
51 c iersl = 1 if a singular matrix arose with miter = 3.
52 c this routine also uses the common variables el0, h, miter, and n.
53 c-----------------------------------------------------------------------
54  iersl = 0
55  SELECT CASE (miter)
56  CASE DEFAULT
57  CALL li_gesl (wm(3:3+n*n), n, n, iwm(21:21+n), x, 0)
58  RETURN
59 c
60  CASE (3)
61  phl0 = wm(2)
62  hl0 = h*el0
63  wm(2) = hl0
64  IF (hl0 /= phl0) THEN
65  r = hl0/phl0
66  DO i = 1, n
67  di = 1.0_dp - r*(1.0_dp - 1.0_dp/wm(i+2))
68  IF (abs(di) == 0.0_dp) go to 390
69  wm(i+2) = 1.0_dp/di
70  END DO
71  END IF
72  x(:n) = wm(3:n+2)*x(:n)
73  RETURN
74  390 CONTINUE
75  iersl = 1
76  RETURN
77 c
78  CASE (4:5)
79  ml = iwm(1)
80  mu = iwm(2)
81  meband = 2*ml + mu + 1
82  CALL li_gbsl (wm(3:3+meband*n), meband, n, ml, mu,
83  1 iwm(21:21+n), x, 0)
84  RETURN
85  END SELECT
86 
87  END SUBROUTINE solsy
liprec::li_gbsl
Definition: liprec.f:37
liprec::li_gesl
Definition: liprec.f:105