GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dslvk.f
Go to the documentation of this file.
1 C Work performed under the auspices of the U.S. Department of Energy
2 C by Lawrence Livermore National Laboratory under contract number
3 C W-7405-Eng-48.
4 C
5  SUBROUTINE dslvk (NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM,
6  * RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
7  * RPAR, IPAR)
8 C
9 C***BEGIN PROLOGUE DSLVK
10 C***REFER TO DDASPK
11 C***DATE WRITTEN 890101 (YYMMDD)
12 C***REVISION DATE 900926 (YYMMDD)
13 C***REVISION DATE 940928 Removed MNEWT and added RHOK in call list.
14 C
15 C
16 C-----------------------------------------------------------------------
17 C***DESCRIPTION
18 C
19 C DSLVK uses a restart algorithm and interfaces to DSPIGM for
20 C the solution of the linear system arising from a Newton iteration.
21 C
22 C In addition to variables described elsewhere,
23 C communication with DSLVK uses the following variables..
24 C WM = Real work space containing data for the algorithm
25 C (Krylov basis vectors, Hessenberg matrix, etc.).
26 C IWM = Integer work space containing data for the algorithm.
27 C X = The right-hand side vector on input, and the solution vector
28 C on output, of length NEQ.
29 C IRES = Error flag from RES.
30 C IERSL = Output flag ..
31 C IERSL = 0 means no trouble occurred (or user RES routine
32 C returned IRES < 0)
33 C IERSL = 1 means the iterative method failed to converge
34 C (DSPIGM returned IFLAG > 0.)
35 C IERSL = -1 means there was a nonrecoverable error in the
36 C iterative solver, and an error exit will occur.
37 C-----------------------------------------------------------------------
38 C***ROUTINES CALLED
39 C DSCAL, DCOPY, DSPIGM
40 C
41 C***END PROLOGUE DSLVK
42 C
43  INTEGER NEQ, IWM, IRES, IERSL, IPAR
44  DOUBLE PRECISION Y, TN, YPRIME, SAVR, X, EWT, WM, CJ, EPLIN,
45  1 sqrtn, rsqrtn, rhok, rpar
46  dimension y(*), yprime(*), savr(*), x(*), ewt(*),
47  1 wm(*), iwm(*), rpar(*), ipar(*)
48 C
49  INTEGER IFLAG, IRST, NRSTS, NRMAX, LR, LDL, LHES, LGMR, LQ, LV,
50  1 LWK, LZ, MAXLP1, NPSL
51  INTEGER NLI, NPS, NCFL, NRE, MAXL, KMP, MITER
52  EXTERNAL RES, PSOL
53 C
54  parameter(lnre=12, lncfl=16, lnli=20, lnps=21)
55  parameter(llocwp=29, llciwp=30)
56  parameter(lmiter=23, lmaxl=24, lkmp=25, lnrmax=26)
57 C
58 C-----------------------------------------------------------------------
59 C IRST is set to 1, to indicate restarting is in effect.
60 C NRMAX is the maximum number of restarts.
61 C-----------------------------------------------------------------------
62  DATA irst/1/
63 C
64  liwp = iwm(llciwp)
65  nli = iwm(lnli)
66  nps = iwm(lnps)
67  ncfl = iwm(lncfl)
68  nre = iwm(lnre)
69  lwp = iwm(llocwp)
70  maxl = iwm(lmaxl)
71  kmp = iwm(lkmp)
72  nrmax = iwm(lnrmax)
73  miter = iwm(lmiter)
74  iersl = 0
75  ires = 0
76 C-----------------------------------------------------------------------
77 C Use a restarting strategy to solve the linear system
78 C P*X = -F. Parse the work vector, and perform initializations.
79 C Note that zero is the initial guess for X.
80 C-----------------------------------------------------------------------
81  maxlp1 = maxl + 1
82  lv = 1
83  lr = lv + neq*maxl
84  lhes = lr + neq + 1
85  lq = lhes + maxl*maxlp1
86  lwk = lq + 2*maxl
87  ldl = lwk + min0(1,maxl-kmp)*neq
88  lz = ldl + neq
89  CALL dscal (neq, rsqrtn, ewt, 1)
90  CALL dcopy (neq, x, 1, wm(lr), 1)
91  DO 110 i = 1,neq
92  110 x(i) = 0.d0
93 C-----------------------------------------------------------------------
94 C Top of loop for the restart algorithm. Initial pass approximates
95 C X and sets up a transformed system to perform subsequent restarts
96 C to update X. NRSTS is initialized to -1, because restarting
97 C does not occur until after the first pass.
98 C Update NRSTS; conditionally copy DL to R; call the DSPIGM
99 C algorithm to solve A*Z = R; updated counters; update X with
100 C the residual solution.
101 C Note: if convergence is not achieved after NRMAX restarts,
102 C then the linear solver is considered to have failed.
103 C-----------------------------------------------------------------------
104  nrsts = -1
105  115 CONTINUE
106  nrsts = nrsts + 1
107  IF (nrsts .GT. 0) CALL dcopy (neq, wm(ldl), 1, wm(lr),1)
108  CALL dspigm (neq, tn, y, yprime, savr, wm(lr), ewt, maxl, maxlp1,
109  1 kmp, eplin, cj, res, ires, nres, psol, npsl, wm(lz), wm(lv),
110  2 wm(lhes), wm(lq), lgmr, wm(lwp), iwm(liwp), wm(lwk),
111  3 wm(ldl), rhok, iflag, irst, nrsts, rpar, ipar)
112  nli = nli + lgmr
113  nps = nps + npsl
114  nre = nre + nres
115  DO 120 i = 1,neq
116  120 x(i) = x(i) + wm(lz+i-1)
117  IF ((iflag .EQ. 1) .AND. (nrsts .LT. nrmax) .AND. (ires .EQ. 0))
118  1 GO TO 115
119 C-----------------------------------------------------------------------
120 C The restart scheme is finished. Test IRES and IFLAG to see if
121 C convergence was not achieved, and set flags accordingly.
122 C-----------------------------------------------------------------------
123  IF (ires .LT. 0) THEN
124  ncfl = ncfl + 1
125  ELSE IF (iflag .NE. 0) THEN
126  ncfl = ncfl + 1
127  IF (iflag .GT. 0) iersl = 1
128  IF (iflag .LT. 0) iersl = -1
129  ENDIF
130 C-----------------------------------------------------------------------
131 C Update IWM with counters, rescale EWT, and return.
132 C-----------------------------------------------------------------------
133  iwm(lnli) = nli
134  iwm(lnps) = nps
135  iwm(lncfl) = ncfl
136  iwm(lnre) = nre
137  CALL dscal (neq, sqrtn, ewt, 1)
138  RETURN
139 C
140 C------END OF SUBROUTINE DSLVK------------------------------------------
141  END
subroutine dslvk(NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM, RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK, RPAR, IPAR)
Definition: dslvk.f:8
subroutine dspigm(NEQ, TN, Y, YPRIME, SAVR, R, WGHT, MAXL, MAXLP1, KMP, EPLIN, CJ, RES, IRES, NRE, PSOL, NPSL, Z, V, HES, Q, LGMR, WP, IWP, WK, DL, RHOK, IFLAG, IRST, NRSTS, RPAR, IPAR)
Definition: dspigm.f:9