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
ddaini.f
Go to the documentation of this file.
1  SUBROUTINE ddaini (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR,
2  + ipar, phi, delta, e, wm, iwm, hmin, uround, nonneg, ntemp)
3 C***BEGIN PROLOGUE DDAINI
4 C***SUBSIDIARY
5 C***PURPOSE Initialization routine for DDASSL.
6 C***LIBRARY SLATEC (DASSL)
7 C***TYPE DOUBLE PRECISION (SDAINI-S, DDAINI-D)
8 C***AUTHOR PETZOLD, LINDA R., (LLNL)
9 C***DESCRIPTION
10 C-----------------------------------------------------------------
11 C DDAINI TAKES ONE STEP OF SIZE H OR SMALLER
12 C WITH THE BACKWARD EULER METHOD, TO
13 C FIND YPRIME. X AND Y ARE UPDATED TO BE CONSISTENT WITH THE
14 C NEW STEP. A MODIFIED DAMPED NEWTON ITERATION IS USED TO
15 C SOLVE THE CORRECTOR ITERATION.
16 C
17 C THE INITIAL GUESS FOR YPRIME IS USED IN THE
18 C PREDICTION, AND IN FORMING THE ITERATION
19 C MATRIX, BUT IS NOT INVOLVED IN THE
20 C ERROR TEST. THIS MAY HAVE TROUBLE
21 C CONVERGING IF THE INITIAL GUESS IS NO
22 C GOOD, OR IF G(X,Y,YPRIME) DEPENDS
23 C NONLINEARLY ON YPRIME.
24 C
25 C THE PARAMETERS REPRESENT:
26 C X -- INDEPENDENT VARIABLE
27 C Y -- SOLUTION VECTOR AT X
28 C YPRIME -- DERIVATIVE OF SOLUTION VECTOR
29 C NEQ -- NUMBER OF EQUATIONS
30 C H -- STEPSIZE. IMDER MAY USE A STEPSIZE
31 C SMALLER THAN H.
32 C WT -- VECTOR OF WEIGHTS FOR ERROR
33 C CRITERION
34 C IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS
35 C IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
36 C IDID=-12 -- DDAINI FAILED TO FIND YPRIME
37 C RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS
38 C THAT ARE NOT ALTERED BY DDAINI
39 C PHI -- WORK SPACE FOR DDAINI
40 C DELTA,E -- WORK SPACE FOR DDAINI
41 C WM,IWM -- REAL AND INTEGER ARRAYS STORING
42 C MATRIX INFORMATION
43 C
44 C-----------------------------------------------------------------
45 C***ROUTINES CALLED DDAJAC, DDANRM, DDASLV
46 C***REVISION HISTORY (YYMMDD)
47 C 830315 DATE WRITTEN
48 C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
49 C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
50 C 901026 Added explicit declarations for all variables and minor
51 C cosmetic changes to prologue. (FNF)
52 C 901030 Minor corrections to declarations. (FNF)
53 C***END PROLOGUE DDAINI
54 C
55  INTEGER neq, idid, ipar(*), iwm(*), nonneg, ntemp
56  double precision
57  * x, y(*), yprime(*), h, wt(*), rpar(*), phi(neq,*), delta(*),
58  * e(*), wm(*), hmin, uround
59  EXTERNAL res, jac
60 C
61  EXTERNAL ddajac, ddanrm, ddaslv
62  DOUBLE PRECISION ddanrm
63 C
64  INTEGER i, ier, ires, jcalc, lnje, lnre, m, maxit, mjac, ncf,
65  * nef, nsf
66  double precision
67  * cj, damp, delnrm, err, oldnrm, r, rate, s, xold, ynorm
68  LOGICAL convgd
69 C
70  parameter(lnre=12)
71  parameter(lnje=13)
72 C
73  DATA maxit/10/,mjac/5/
74  DATA damp/0.75d0/
75 C
76 C
77 C---------------------------------------------------
78 C BLOCK 1.
79 C INITIALIZATIONS.
80 C---------------------------------------------------
81 C
82 C***FIRST EXECUTABLE STATEMENT DDAINI
83  idid=1
84  nef=0
85  ncf=0
86  nsf=0
87  xold=x
88  ynorm=ddanrm(neq,y,wt,rpar,ipar)
89 C
90 C SAVE Y AND YPRIME IN PHI
91  DO 100 i=1,neq
92  phi(i,1)=y(i)
93 100 phi(i,2)=yprime(i)
94 C
95 C
96 C----------------------------------------------------
97 C BLOCK 2.
98 C DO ONE BACKWARD EULER STEP.
99 C----------------------------------------------------
100 C
101 C SET UP FOR START OF CORRECTOR ITERATION
102 200 cj=1.0d0/h
103  x=x+h
104 C
105 C PREDICT SOLUTION AND DERIVATIVE
106  DO 250 i=1,neq
107 250 y(i)=y(i)+h*yprime(i)
108 C
109  jcalc=-1
110  m=0
111  convgd=.true.
112 C
113 C
114 C CORRECTOR LOOP.
115 300 iwm(lnre)=iwm(lnre)+1
116  ires=0
117 C
118  CALL res(x,y,yprime,delta,ires,rpar,ipar)
119  IF (ires.LT.0) go to 430
120 C
121 C
122 C EVALUATE THE ITERATION MATRIX
123  IF (jcalc.NE.-1) go to 310
124  iwm(lnje)=iwm(lnje)+1
125  jcalc=0
126  CALL ddajac(neq,x,y,yprime,delta,cj,h,
127  * ier,wt,e,wm,iwm,res,ires,
128  * uround,jac,rpar,ipar,ntemp)
129 C
130  s=1000000.d0
131  IF (ires.LT.0) go to 430
132  IF (ier.NE.0) go to 430
133  nsf=0
134 C
135 C
136 C
137 C MULTIPLY RESIDUAL BY DAMPING FACTOR
138 310 CONTINUE
139  DO 320 i=1,neq
140 320 delta(i)=delta(i)*damp
141 C
142 C COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
143 C STORE THE CORRECTION IN DELTA
144 C
145  CALL ddaslv(neq,delta,wm,iwm)
146 C
147 C UPDATE Y AND YPRIME
148  DO 330 i=1,neq
149  y(i)=y(i)-delta(i)
150 330 yprime(i)=yprime(i)-cj*delta(i)
151 C
152 C TEST FOR CONVERGENCE OF THE ITERATION.
153 C
154  delnrm=ddanrm(neq,delta,wt,rpar,ipar)
155  IF (delnrm.LE.100.d0*uround*ynorm)
156  * go to 400
157 C
158  IF (m.GT.0) go to 340
159  oldnrm=delnrm
160  go to 350
161 C
162 340 rate=(delnrm/oldnrm)**(1.0d0/m)
163  IF (rate.GT.0.90d0) go to 430
164  s=rate/(1.0d0-rate)
165 C
166 350 IF (s*delnrm .LE. 0.33d0) go to 400
167 C
168 C
169 C THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
170 C M AND AND TEST WHETHER THE MAXIMUM
171 C NUMBER OF ITERATIONS HAVE BEEN TRIED.
172 C EVERY MJAC ITERATIONS, GET A NEW
173 C ITERATION MATRIX.
174 C
175  m=m+1
176  IF (m.GE.maxit) go to 430
177 C
178  IF ((m/mjac)*mjac.EQ.m) jcalc=-1
179  go to 300
180 C
181 C
182 C THE ITERATION HAS CONVERGED.
183 C CHECK NONNEGATIVITY CONSTRAINTS
184 400 IF (nonneg.EQ.0) go to 450
185  DO 410 i=1,neq
186 410 delta(i)=min(y(i),0.0d0)
187 C
188  delnrm=ddanrm(neq,delta,wt,rpar,ipar)
189  IF (delnrm.GT.0.33d0) go to 430
190 C
191  DO 420 i=1,neq
192  y(i)=y(i)-delta(i)
193 420 yprime(i)=yprime(i)-cj*delta(i)
194  go to 450
195 C
196 C
197 C EXITS FROM CORRECTOR LOOP.
198 430 convgd=.false.
199 450 IF (.NOT.convgd) go to 600
200 C
201 C
202 C
203 C-----------------------------------------------------
204 C BLOCK 3.
205 C THE CORRECTOR ITERATION CONVERGED.
206 C DO ERROR TEST.
207 C-----------------------------------------------------
208 C
209  DO 510 i=1,neq
210 510 e(i)=y(i)-phi(i,1)
211  err=ddanrm(neq,e,wt,rpar,ipar)
212 C
213  IF (err.LE.1.0d0) RETURN
214 C
215 C
216 C
217 C--------------------------------------------------------
218 C BLOCK 4.
219 C THE BACKWARD EULER STEP FAILED. RESTORE X, Y
220 C AND YPRIME TO THEIR ORIGINAL VALUES.
221 C REDUCE STEPSIZE AND TRY AGAIN, IF
222 C POSSIBLE.
223 C---------------------------------------------------------
224 C
225 600 CONTINUE
226  x = xold
227  DO 610 i=1,neq
228  y(i)=phi(i,1)
229 610 yprime(i)=phi(i,2)
230 C
231  IF (convgd) go to 640
232  IF (ier.EQ.0) go to 620
233  nsf=nsf+1
234  h=h*0.25d0
235  IF (nsf.LT.3.AND.abs(h).GE.hmin) go to 690
236  idid=-12
237  RETURN
238 620 IF (ires.GT.-2) go to 630
239  idid=-12
240  RETURN
241 630 ncf=ncf+1
242  h=h*0.25d0
243  IF (ncf.LT.10.AND.abs(h).GE.hmin) go to 690
244  idid=-12
245  RETURN
246 C
247 640 nef=nef+1
248  r=0.90d0/(2.0d0*err+0.0001d0)
249  r=max(0.1d0,min(0.5d0,r))
250  h=h*r
251  IF (abs(h).GE.hmin.AND.nef.LT.10) go to 690
252  idid=-12
253  RETURN
254 690 go to 200
255 C
256 C-------------END OF SUBROUTINE DDAINI----------------------
257  END