1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
subroutine solsy (wm, iwm, x, tem)
clll. optimize
integer iwm
integer iownd, iowns,
1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
integer i, meband, ml, mu
double precision wm, x, tem
double precision rowns,
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
double precision di, hl0, phl0, r
dimension wm(*), iwm(*), x(1), tem(1)
common /ls0001/ rowns(209),
2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
3 iownd(14), iowns(6),
4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
c-----------------------------------------------------------------------
c this routine manages the solution of the linear system arising from
c a chord iteration. it is called if miter .ne. 0.
c if miter is 1 or 2, it calls dgetrs to accomplish this.
c if miter = 3 it updates the coefficient h*el0 in the diagonal
c matrix, and then computes the solution.
c if miter is 4 or 5, it calls dgbtrs.
c communication with solsy uses the following variables..
c wm = real work space containing the inverse diagonal matrix if
c miter = 3 and the lu decomposition of the matrix otherwise.
c storage of matrix elements starts at wm(3).
c wm also contains the following matrix-related data..
c wm(1) = sqrt(uround) (not used here),
c wm(2) = hl0, the previous value of h*el0, used if miter = 3.
c iwm = integer work space containing pivot information, starting at
c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band
c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
c x = the right-hand side vector on input, and the solution vector
c on output, of length n.
c tem = vector of work space of length n, not used in this version.
c iersl = output flag (in common). iersl = 0 if no trouble occurred.
c iersl = 1 if a singular matrix arose with miter = 3.
c this routine also uses the common variables el0, h, miter, and n.
c-----------------------------------------------------------------------
iersl = 0
go to (100, 100, 300, 400, 400), miter
c Replaced LINPACK dgesl with LAPACK dgetrs
c 100 call dgesl (wm(3), n, n, iwm(21), x, 0)
100 call dgetrs ('N', n, 1, wm(3), n, iwm(21), x, n, ier)
return
c
300 phl0 = wm(2)
hl0 = h*el0
wm(2) = hl0
if (hl0 .eq. phl0) go to 330
r = hl0/phl0
do 320 i = 1,n
di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2))
if (dabs(di) .eq. 0.0d0) go to 390
320 wm(i+2) = 1.0d0/di
330 do 340 i = 1,n
340 x(i) = wm(i+2)*x(i)
return
390 iersl = 1
return
c
400 ml = iwm(1)
mu = iwm(2)
meband = 2*ml + mu + 1
c Replaced LINPACK dgbsl with LAPACK dgbtrs
c call dgbsl (wm(3), meband, n, ml, mu, iwm(21), x, 0)
call dgbtrs ('N', n, ml, mu, 1, wm(3), meband, iwm(21), x, n, ier)
return
c----------------------- end of subroutine solsy -----------------------
end
|
