/* Copyright (C) 2019-2020, Dirk Krause SPDX-License-Identifier: BSD-3-Clause */ /* WARNING: This file was generated by the dkct program (see http://dktools.sourceforge.net/ for details). Changes you make here will be lost if dkct is run again! You should modify the original source and run dkct on it. Original source: dk4itmin.ctr */ /** @file dk4itmin.c The dk4itmin module. */ #line 68 "dk4itmin.ctr" #include "dk4conf.h" #include #if DK4_HAVE_ASSERT_H #ifndef ASSERT_H_INCLUDED #include #define ASSERT_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDLIB_H #ifndef STDLIB_H_INCLUDED #include #define STDLIB_H_INCLUDED 1 #endif #endif #if DK4_HAVE_LIMITS_H #ifndef LIMITS_H_INCLUDED #include #define LIMITS_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDINT_H #ifndef STDINT_H_INCLUDED #include #define STDINT_H_INCLUDED 1 #endif #endif #ifndef DK4MEM_H_INCLUDED #include "dk4mem.h" #endif #ifndef DK4ITER_H_INCLUDED #include "dk4iter.h" #endif #ifndef DK4MIN_H_INCLUDED #include "dk4itmin.h" #endif #ifndef DK4MATH_H_INCLUDED #include "dk4math.h" #endif static int dk4itmin_check_tolerances( double eps_y, double y, double yold, double eps_x, double x, double xold ) { int back = 0; if ((!(isgreater(eps_y, 0.0))) || (isgreater(eps_y, fabs(y - yold)))) { if ((!(isgreater(eps_x, 0.0))) || (isgreater(eps_x, fabs(x - xold)))) { back = 1; } } return back; } static int dk4itmin_golden_section( double *pd, unsigned long *pp, dk4_iter_fct_t *fct, void const *ps, double a, double fa, double b, double fb, dk4_iter_ctx_t const *ctx ) { double c = 0.0; /* Probe position c */ double fc = 0.0; /* Function value at c */ double d = 0.0; /* Probe position d */ double fd = 0.0; /* Function value at d */ double f1 = 0.0; /* Factor for probe calculation */ double f2 = 0.0; /* Factor for probe calculation */ double x = 0.0; /* Current probe point */ double xold = 0.0; /* Previous probe point */ double y = 0.0; /* Current probe value */ double yold = 0.0; /* Previous probe value */ double max = 0.0; /* Maximum of fa and fb */ unsigned long passno = 0UL; /* Current pass number */ int res = 0; /* Operations result */ int cc = 1; /* 1=go on, 0=finished, -1=abort */ int back = DK4_ITER_RESULT_E_ARGS; #line 171 "dk4itmin.ctr" /* Initialize allowed maximum function value */ /* Initialize factors for probe calculation */ f1 = sqrt(5.0) - 1.0; #line 178 "dk4itmin.ctr" f2 = 3.0 - sqrt(5.0); #line 179 "dk4itmin.ctr" max = ((fa > fb) ? (fa) : (fb)); /* Initial probe points */ c = (f1 * a + f2 * b) / 2.0; d = (f2 * a + f1 * b) / 2.0; res = (*fct)(&fc, c, ps); if ((0 == res) || (!(dk4ma_is_finite(fc)))) { #line 187 "dk4itmin.ctr" back = DK4_ITER_RESULT_E_FCT; goto finished; } res = (*fct)(&fd, d, ps); if ((0 == res) || (!(dk4ma_is_finite(fd)))) { #line 192 "dk4itmin.ctr" back = DK4_ITER_RESULT_E_FCT; goto finished; } #line 196 "dk4itmin.ctr" #line 197 "dk4itmin.ctr" #line 198 "dk4itmin.ctr" #line 199 "dk4itmin.ctr" if (fc > max) { #line 200 "dk4itmin.ctr" back = DK4_ITER_RESULT_E_CONV; goto finished; } if (fd > max) { #line 204 "dk4itmin.ctr" back = DK4_ITER_RESULT_E_CONV; goto finished; } /* Iteration loop */ do { /* Increase pass number */ if (ULONG_MAX > passno) { passno++; } #line 215 "dk4itmin.ctr" #line 216 "dk4itmin.ctr" #line 217 "dk4itmin.ctr" #line 218 "dk4itmin.ctr" #line 219 "dk4itmin.ctr" /* Save probe point and value */ xold = x; yold = y; /* Process the probe points */ if (fc < fd) { #line 228 "dk4itmin.ctr" b = d; d = c; c = (f1 * a + f2 * b) / 2.0; fb = fd; fd = fc; res = (*fct)(&fc, c, ps); if ((0 == res) || (!(dk4ma_is_finite(fc)))) { back = DK4_ITER_RESULT_E_FCT; #line 236 "dk4itmin.ctr" cc = -1; } else { max = ((fa > fb) ? (fa) : (fb)); if (fc > max) { #line 241 "dk4itmin.ctr" back = DK4_ITER_RESULT_E_CONV; cc = -1; } else { x = c; #line 246 "dk4itmin.ctr" y = fc; #line 247 "dk4itmin.ctr" } } } else { #line 251 "dk4itmin.ctr" a = c; c = d; d = (f2 * a + f1 * b) / 2.0; fa = fc; fc = fd; res = (*fct)(&fd, d, ps); if ((0 == res) || (!(dk4ma_is_finite(fd)))) { back = DK4_ITER_RESULT_E_FCT; #line 259 "dk4itmin.ctr" cc = -1; } else { max = ((fa > fb) ? (fa) : (fb)); if (fd > max) { #line 264 "dk4itmin.ctr" back = DK4_ITER_RESULT_E_CONV; cc = -1; } else { x = d; #line 269 "dk4itmin.ctr" y = fd; #line 270 "dk4itmin.ctr" } } } /* Check for reached tolerance */ if (1 == cc) { if (1UL < passno) { if (x == xold) { cc = 0; #line 280 "dk4itmin.ctr" } else { if (0 == ctx->exact) { res = dk4itmin_check_tolerances( ctx->eps_y, y, yold, ctx->eps_x, x, xold ); if (0 != res) { cc = 0; #line 288 "dk4itmin.ctr" } #if TRACE_DEBUG else { #line 291 "dk4itmin.ctr" } #endif } } } #if TRACE_DEBUG else { #line 298 "dk4itmin.ctr" } #endif } /* Abort if number of passes exceeded */ if (1 == cc) { if ((0UL < ctx->maxpass) && (passno >= ctx->maxpass)) { back = DK4_ITER_RESULT_E_PASSES; cc = -1; #line 308 "dk4itmin.ctr" } #if TRACE_DEBUG else { #line 311 "dk4itmin.ctr" } #endif } } while (1 == cc); /* Transfer results on success */ if (0 == cc) { #line 320 "dk4itmin.ctr" pd[0] = x; pd[1] = y; back = DK4_ITER_RESULT_SUCCESS; if (NULL != pp) { *pp = passno; } } finished: #line 328 "dk4itmin.ctr" return back; } int dk4itmin_interval( double *d, unsigned long *pp, dk4_iter_fct_t *fct, void const *ps, double a, double b, dk4_iter_ctx_t const *ctx ) { dk4_iter_ctx_t mctx; /* Iteration context copy */ double fa = 0.0; /* Value at left border */ double fb = 0.0; /* Value at right border */ int res = 0; /* Operation result */ int back = DK4_ITER_RESULT_E_ARGS; #line 350 "dk4itmin.ctr" /* Check function cal arguments */ #if DK4_USE_ASSERT assert(NULL != d) assert(NULL != fct) assert(a < b) #endif if ((NULL == d) || (NULL == fct)) { #line 359 "dk4itmin.ctr" goto finished; } if (a >= b) { #line 362 "dk4itmin.ctr" goto finished; } /* Calculate values for the three given points */ res = (*fct)(&fa, a, ps); if ((0 == res) || (!(dk4ma_is_finite(fa)))) { back = DK4_ITER_RESULT_E_FCT; #line 370 "dk4itmin.ctr" goto finished; } res = (*fct)(&fb, b, ps); if ((0 == res) || (!(dk4ma_is_finite(fb)))) { back = DK4_ITER_RESULT_E_FCT; #line 375 "dk4itmin.ctr" goto finished; } /* Copy or initialize context */ if (NULL != ctx) { #line 381 "dk4itmin.ctr" DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t)); } else { #line 384 "dk4itmin.ctr" dk4iter_ctx_init(&mctx); } /* Choose algorithm depending on context At this time only golden section available */ back = dk4itmin_golden_section(d, pp, fct, ps, a, fa, b, fb, &mctx); finished: #line 394 "dk4itmin.ctr" return back; } /* vim: set ai sw=4 ts=4 : */