1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1172.c,v 1.3 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1172
50
51 #include "sislP.h"
52
53 /*
54 * Forward declarations.
55 * ---------------------
56 */
57 #if defined(SISLNEEDPROTOTYPES)
58 static void s1172_s9corr(double *,double,double,double);
59 static void s1172_s9dir(double *,double []);
60 #else
61 static void s1172_s9corr();
62 static void s1172_s9dir();
63 #endif
64
65 #if defined(SISLNEEDPROTOTYPES)
66 void
s1172(SISLCurve * pcurve,double astart,double aend,double anext,double * cpos,int * jstat)67 s1172(SISLCurve *pcurve,double astart,
68 double aend, double anext, double *cpos,int *jstat)
69 #else
70 void s1172(pcurve,astart,aend,anext,cpos,jstat)
71 SISLCurve *pcurve;
72 double astart;
73 double aend;
74 double anext;
75 double *cpos;
76 int *jstat;
77 #endif
78 /*
79 *********************************************************************
80 *
81 *********************************************************************
82 *
83 * PURPOSE : Newton iteration on a onedimensional curve.
84 * The function finds a local extremum.
85 *
86 *
87 * INPUT : pcurve - Pointer to the curve.
88 * astart - Start values of parameter intervall.
89 * aend - End value of parameter intervall.
90 * anext - Parameter start value for iteration.
91 *
92 *
93 * OUTPUT : cpos - Parameter value of the found extremum.
94 * jstat - status messages
95 * = 1 : Extremum found.
96 * = 0 : Extremum NOT found.
97 * < 0 : error.
98 *
99 *
100 * METHOD : Newton iteration in one parameter direction.
101 *
102 *
103 * REFERENCES :
104 *
105 *
106 * WRITTEN BY : Ulf J. Krystad, SI, OCTOBER 1990
107 * CORRECTED BY: Ulf J. Krystad, SI, AUGUST 1991
108 *********************************************************************
109 */
110 {
111 int kstat = 0; /* Local status variable. */
112 int kpos = 0; /* Position of error. */
113 int kleft=0; /* Variables used in the evaluator. */
114 int kder=3; /* Order of derivatives to be calulated */
115 int knbit; /* Number of iterations */
116 int kdir; /* Changing direction. */
117 double tdelta; /* Parameter intervals of the Curve. */
118 double tdist; /* Euclidian norm of derivative vector */
119 double tprev; /* Previous Euclidian norm of derivative vector*/
120 double td,t1,tdn; /* Distances between old and new parameter
121 value in the two parameter directions. */
122 double sval[4]; /* Value ,first and second derivatiev of Curve.*/
123 double tnext; /* Parameter values */
124 double tol = (double)1000.0*REL_COMP_RES; /* Singularity tolerance */
125 /* --------------------------------------------------------------------- */
126
127 /* Test input. */
128 if (pcurve->idim != 1) goto err106;
129
130 /* Fetch endpoints and the interval of parameter interval of curves. */
131
132 tdelta = pcurve->et[pcurve->in] - pcurve->et[pcurve->ik - 1];
133
134 /* Evaluate 0-2.st derivatives of curve */
135 s1221(pcurve,kder,anext,&kleft,sval,&kstat);
136 if (kstat < 0) goto error;
137
138 /* Get Euclidian norm of derivative */
139 tprev = fabs(sval[1]);
140
141 /* Compute the Newton stepdistanse vector. */
142 s1172_s9dir(&td,sval);
143
144 /* Adjust if we are not inside the parameter intervall. */
145 t1 = td;
146 s1172_s9corr(&t1,anext,astart,aend);
147
148 /* Iterate to find the intersection point. */
149
150 for (knbit = 0; knbit < 50; knbit++)
151 {
152 /* Evaluate 0-3.st derivatives of curve */
153
154 tnext = anext + t1;
155
156 s1221(pcurve,kder,tnext,&kleft,sval,&kstat);
157 if (kstat < 0) goto error;
158
159 /* Get Euclidian norm of derivative */
160 tdist = fabs(sval[1]);
161
162 /* Compute the Newton stepdistanse vector. */
163 s1172_s9dir(&tdn,sval);
164
165 /* Check if the direction of the step have change. */
166
167 kdir = (td*tdn >= DZERO); /* 0 if changed. */
168
169 if (tdist <= tprev || kdir)
170 {
171 /* Ordinary converging. */
172
173 anext += t1;
174
175 td = t1 = tdn;
176
177 /* Adjust if we are not inside the parameter intervall. */
178 s1172_s9corr(&t1,anext,astart,aend);
179
180
181 if (fabs(t1/tdelta) <= REL_COMP_RES)
182 {
183 anext += t1;
184 break;
185 }
186
187 tprev = tdist;
188 }
189
190 else
191 {
192 /* Not converging, half step length try again. */
193
194 t1 /= (double)2;
195 /* knbit--; */
196 }
197 }
198
199 /* Iteration stopped, test if point is extremum */
200
201 if (tdist <= tol)
202 *jstat = 1;
203 else
204 *jstat = 0;
205
206
207 /* Test if the iteration is close to a knot */
208 if (fabs(anext - pcurve->et[kleft])/tdelta < tol)
209 anext = pcurve->et[kleft];
210 else if (fabs(anext - pcurve->et[kleft+1])/tdelta < tol)
211 anext = pcurve->et[kleft+1];
212
213 /* Uppdate output. */
214 *cpos = anext;
215
216 /* Iteration completed. */
217 goto out;
218
219 /* --------------------------------------------------------------------- */
220 /* Error in input. Dimension not equal to 1 */
221 err106: *jstat = -106;
222 s6err("s1172",*jstat,kpos);
223 goto out;
224
225 /* Error in lower level routine. */
226 error : *jstat = kstat;
227 s6err("s1172",*jstat,kpos);
228 goto out;
229
230 out:;
231 }
232
233 #if defined(SISLNEEDPROTOTYPES)
234 static void
s1172_s9corr(double * cd,double acoef,double astart,double aend)235 s1172_s9corr(double *cd, double acoef,double astart,double aend)
236 #else
237 static void s1172_s9corr(cd,acoef,astart,aend)
238 double *cd;
239 double acoef;
240 double astart;
241 double aend;
242 #endif
243 /*
244 *********************************************************************
245 *
246 *********************************************************************
247 *
248 * PURPOSE : To be sure that we are inside the boorder of the
249 * parameter plan. If we are outside clipping is used
250 * to corrigate the step value.
251 *
252 *
253 * INPUT : acoef - Parameter value to clipp
254 * astart - The lower boorder
255 * aend - The higher boorder
256 *
257 *
258 *
259 * INPUT/OUTPUT : cd - Old and new step value.
260 *
261 *
262 * METHOD : We are cutting a line.
263 * In this case we always know that the startpoint of
264 * the line is inside the rectangel, and we may therfor
265 * use a simple kind of clipping.
266 *
267 *
268 * REFERENCES :
269 *
270 *
271 * WRITTEN BY : Arne Laksaa, SI, Feb 1989
272 *
273 *********************************************************************
274 */
275 {
276 if (acoef + *cd < astart) *cd = astart - acoef;
277 else if (acoef + *cd > aend) *cd = aend - acoef;
278 }
279
280 #if defined(SISLNEEDPROTOTYPES)
281 static void
s1172_s9dir(double * cdiff,double evals[])282 s1172_s9dir(double *cdiff,double evals[])
283 #else
284 static void s1172_s9dir(cdiff,evals)
285 double *cdiff;
286 double evals[];
287 #endif
288 /*
289 *********************************************************************
290 *
291 *********************************************************************
292 *
293 * PURPOSE : To compute the step according to a Newton scheme
294 * for finding an extremal to a one dimensional curve.
295 *
296 *
297 * INPUT : evals - Value and derivatives in point on curve.
298 *
299 *
300 * OUTPUT : cdiff - Parameter increment in first direction.
301 *
302 *
303 *
304 * METHOD : This is a one dimensional case. We want to find x such that
305 *
306 * x : S'(x) = 0
307 *
308 * Using Taylor we get:
309 *
310 * x : S'+dxS''+0.5 dxdxS''' = 0
311 *
312 * The solution of this equation is the
313 * following function.
314 *
315 * REFERENCES :
316 *
317 *
318 * WRITTEN BY : Ulf J. Krystad, SI, OCTOBER 1990
319 *
320 *********************************************************************
321 */
322 {
323 double a,b,c,d,d1,d2;
324
325 a = evals[3];
326 b = evals[2];
327 c = b*b - 2.0*a*evals[1];
328
329 if (fabs(b) > DZERO) d = -evals[1]/b;
330 else d = 0.0;
331
332
333 if (c < DZERO) *cdiff = d;
334 else if (fabs(a) > DZERO)
335 {
336 c = sqrt(c);
337 d1 = (-b + c)/a;
338 d2 = (-b - c)/a;
339 if (DEQUAL(b,c)) *cdiff = d;
340 else
341 if (fabs(d1-d) < fabs(d2-d)) *cdiff = d1;
342 else *cdiff = d2;
343 }
344 else *cdiff = d;
345 }
346
347