GNU Octave  4.2.1
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
ddasik.f
Go to the documentation of this file.
1 C Work perfored 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 ddasik(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACK,PSOL,H,WT,
6  * jskip,rpar,ipar,savr,delta,r,yic,ypic,pwk,wm,iwm,cj,uround,
7  * epli,sqrtn,rsqrtn,epcon,ratemx,stptol,jflg,
8  * icnflg,icnstr,iernls)
9 C
10 C***BEGIN PROLOGUE DDASIK
11 C***REFER TO DDASPK
12 C***DATE WRITTEN 941026 (YYMMDD)
13 C***REVISION DATE 950808 (YYMMDD)
14 C***REVISION DATE 951110 Removed unreachable block 390.
15 C
16 C
17 C-----------------------------------------------------------------------
18 C***DESCRIPTION
19 C
20 C
21 C DDASIK solves a nonlinear system of algebraic equations of the
22 C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
23 C the initial conditions.
24 C
25 C An initial value for Y and initial guess for YPRIME are input.
26 C
27 C The method used is a Newton scheme with Krylov iteration and a
28 C linesearch algorithm.
29 C
30 C The parameters represent
31 C
32 C X -- Independent variable.
33 C Y -- Solution vector at x.
34 C YPRIME -- Derivative of solution vector.
35 C NEQ -- Number of equations to be integrated.
36 C ICOPT -- Initial condition option chosen (1 or 2).
37 C ID -- Array of dimension NEQ, which must be initialized
38 C if ICOPT = 1. See DDASIC.
39 C RES -- External user-supplied subroutine
40 C to evaluate the residual. See RES description
41 C in DDASPK prologue.
42 C JACK -- External user-supplied routine to update
43 C the preconditioner. (This is optional).
44 C See JAC description for the case
45 C INFO(12) = 1 in the DDASPK prologue.
46 C PSOL -- External user-supplied routine to solve
47 C a linear system using preconditioning.
48 C (This is optional). See explanation inside DDASPK.
49 C H -- Scaling factor for this initial condition calc.
50 C WT -- Vector of weights for error criterion.
51 C JSKIP -- input flag to signal if initial JAC call is to be
52 C skipped. 1 => skip the call, 0 => do not skip call.
53 C RPAR,IPAR -- Real and integer arrays used for communication
54 C between the calling program and external user
55 C routines. They are not altered within DASPK.
56 C SAVR -- Work vector for DDASIK of length NEQ.
57 C DELTA -- Work vector for DDASIK of length NEQ.
58 C R -- Work vector for DDASIK of length NEQ.
59 C YIC,YPIC -- Work vectors for DDASIK, each of length NEQ.
60 C PWK -- Work vector for DDASIK of length NEQ.
61 C WM,IWM -- Real and integer arrays storing
62 C matrix information for linear system
63 C solvers, and various other information.
64 C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
65 C UROUND -- Unit roundoff.
66 C EPLI -- convergence test constant.
67 C See DDASPK prologue for more details.
68 C SQRTN -- Square root of NEQ.
69 C RSQRTN -- reciprical of square root of NEQ.
70 C EPCON -- Tolerance to test for convergence of the Newton
71 C iteration.
72 C RATEMX -- Maximum convergence rate for which Newton iteration
73 C is considered converging.
74 C JFLG -- Flag showing whether a Jacobian routine is supplied.
75 C ICNFLG -- Integer scalar. If nonzero, then constraint
76 C violations in the proposed new approximate solution
77 C will be checked for, and the maximum step length
78 C will be adjusted accordingly.
79 C ICNSTR -- Integer array of length NEQ containing flags for
80 C checking constraints.
81 C IERNLS -- Error flag for nonlinear solver.
82 C 0 ==> nonlinear solver converged.
83 C 1,2 ==> recoverable error inside nonlinear solver.
84 C 1 => retry with current Y, YPRIME
85 C 2 => retry with original Y, YPRIME
86 C -1 ==> unrecoverable error in nonlinear solver.
87 C
88 C-----------------------------------------------------------------------
89 C
90 C***ROUTINES CALLED
91 C RES, JACK, DNSIK, DCOPY
92 C
93 C***END PROLOGUE DDASIK
94 C
95 C
96  IMPLICIT DOUBLE PRECISION(a-h,o-z)
97  dimension y(*),yprime(*),id(*),wt(*),icnstr(*)
98  dimension savr(*),delta(*),r(*),yic(*),ypic(*),pwk(*)
99  dimension wm(*),iwm(*), rpar(*),ipar(*)
100  EXTERNAL res, jack, psol
101 C
102  parameter(lnre=12, lnje=13, llocwp=29, llciwp=30)
103  parameter(lmxnit=32, lmxnj=33)
104 C
105 C
106 C Perform initializations.
107 C
108  lwp = iwm(llocwp)
109  liwp = iwm(llciwp)
110  mxnit = iwm(lmxnit)
111  mxnj = iwm(lmxnj)
112  iernls = 0
113  nj = 0
114  eplin = epli*epcon
115 C
116 C Call RES to initialize DELTA.
117 C
118  ires = 0
119  iwm(lnre) = iwm(lnre) + 1
120  CALL res(x,y,yprime,cj,delta,ires,rpar,ipar)
121  IF (ires .LT. 0) go to 370
122 C
123 C Looping point for updating the preconditioner.
124 C
125  300 CONTINUE
126 C
127 C Initialize all error flags to zero.
128 C
129  ierpj = 0
130  ires = 0
131  iernew = 0
132 C
133 C If a Jacobian routine was supplied, call it.
134 C
135  IF (jflg .EQ. 1 .AND. jskip .EQ. 0) THEN
136  nj = nj + 1
137  iwm(lnje)=iwm(lnje)+1
138  CALL jack(res, ires, neq, x, y, yprime, wt, delta, r, h, cj,
139  * wm(lwp), iwm(liwp), ierpj, rpar, ipar)
140  IF (ires .LT. 0 .OR. ierpj .NE. 0) go to 370
141  ENDIF
142  jskip = 0
143 C
144 C Call the nonlinear Newton solver for up to MXNIT iterations.
145 C
146  CALL dnsik(x,y,yprime,neq,icopt,id,res,psol,wt,rpar,ipar,
147  * savr,delta,r,yic,ypic,pwk,wm,iwm,cj,sqrtn,rsqrtn,
148  * eplin,epcon,ratemx,mxnit,stptol,icnflg,icnstr,iernew)
149 C
150  IF (iernew .EQ. 1 .AND. nj .LT. mxnj .AND. jflg .EQ. 1) THEN
151 C
152 C Up to MXNIT iterations were done, the convergence rate is < 1,
153 C a Jacobian routine is supplied, and the number of JACK calls
154 C is less than MXNJ.
155 C Copy the residual SAVR to DELTA, call JACK, and try again.
156 C
157  CALL dcopy(neq, savr, 1, delta, 1)
158  go to 300
159  ENDIF
160 C
161  IF (iernew .NE. 0) go to 380
162  RETURN
163 C
164 C
165 C Unsuccessful exits from nonlinear solver.
166 C Set IERNLS accordingly.
167 C
168  370 iernls = 2
169  IF (ires .LE. -2) iernls = -1
170  RETURN
171 C
172  380 iernls = min(iernew,2)
173  RETURN
174 C
175 C----------------------- END OF SUBROUTINE DDASIK-----------------------
176  END
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
double h
Definition: graphics.cc:11205
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
subroutine ddasik(X, Y, YPRIME, NEQ, ICOPT, ID, RES, JACK, PSOL, H, WT, JSKIP, RPAR, IPAR, SAVR, DELTA, R, YIC, YPIC, PWK, WM, IWM, CJ, UROUND, EPLI, SQRTN, RSQRTN, EPCON, RATEMX, STPTOL, JFLG, ICNFLG, ICNSTR, IERNLS)
Definition: ddasik.f:5
subroutine dnsik(X, Y, YPRIME, NEQ, ICOPT, ID, RES, PSOL, WT, RPAR, IPAR, SAVR, DELTA, R, YIC, YPIC, PWK, WM, IWM, CJ, SQRTN, RSQRTN, EPLIN, EPCON, RATEMX, MAXIT, STPTOL, ICNFLG, ICNSTR, IERNEW)
Definition: dnsik.f:5
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:205