summaryrefslogtreecommitdiff
path: root/support/dktools/dk4iter.c
diff options
context:
space:
mode:
Diffstat (limited to 'support/dktools/dk4iter.c')
-rw-r--r--support/dktools/dk4iter.c1226
1 files changed, 0 insertions, 1226 deletions
diff --git a/support/dktools/dk4iter.c b/support/dktools/dk4iter.c
deleted file mode 100644
index d958ac2b1e..0000000000
--- a/support/dktools/dk4iter.c
+++ /dev/null
@@ -1,1226 +0,0 @@
-/*
-Copyright (C) 2018-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: dk4iter.ctr
-*/
-
-/** @file dk4iter.c The dk4iter module.
-*/
-
-
-#line 401 "dk4iter.ctr"
-
-
-#include "dk4conf.h"
-
-#include <stdio.h>
-
-#if DK4_HAVE_ASSERT_H
-#ifndef ASSERT_H_INCLUDED
-#include <assert.h>
-#define ASSERT_H_INCLUDED 1
-#endif
-#endif
-
-#if DK4_HAVE_STDLIB_H
-#ifndef STDLIB_H_INCLUDED
-#include <stdlib.h>
-#define STDLIB_H_INCLUDED 1
-#endif
-#endif
-
-#if DK4_HAVE_LIMITS_H
-#ifndef LIMITS_H_INCLUDED
-#include <limits.h>
-#define LIMITS_H_INCLUDED 1
-#endif
-#endif
-
-#if DK4_HAVE_STDINT_H
-#ifndef STDINT_H_INCLUDED
-#include <stdint.h>
-#define STDINT_H_INCLUDED 1
-#endif
-#endif
-
-
-#include "dk4mem.h"
-#include "dk4iter.h"
-#include "dk4math.h"
-
-
-
-
-
-#line 443 "dk4iter.ctr"
-
-
-
-/** Initial function value condition.
-*/
-enum {
- CONDITION_ILLEGAL = 0 ,
- /** fa < 0, fb >= 0.
- */
- CONDITION_FA_LESS_ZERO ,
-
- /** fa <= 0, fb > 0.
- */
- CONDITION_FA_LEQ_ZERO ,
-
- /** fa > 0, fb <= 0.
- */
- CONDITION_FA_GREATER_ZERO ,
-
- /** fa >= 0, fb < 0.
- */
- CONDITION_FA_GEQ_ZERO
-};
-
-
-
-/** Note which interval border was changed in previous step.
-*/
-enum {
- /** No interval border change yet.
- */
- CHANGED_NONE = 0,
-
- /** Previous step changed interval border a.
- */
- CHANGED_A = -1,
-
- /** Previous step changed interval border b.
- */
- CHANGED_B = 1
-};
-
-
-
-/** Flags which restrictions apply to x values.
-*/
-enum {
- /** Minimum specified.
- */
- MINMAX_MINIMUM = 0x0001 ,
-
- /** Maximum specified.
- */
- MINMAX_MAXIMUM = 0x0002 ,
-
- /** Specified minimum is exclusive.
- */
- MINMAX_MIN_EXCL = 0x0004 ,
-
- /** Specified maximum is exclusive.
- */
- MINMAX_MAX_EXCL = 0x0008
-};
-
-
-
-void
-dk4iter_ctx_init(dk4_iter_ctx_t *ctx)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- DK4_MEMRES(ctx,sizeof(dk4_iter_ctx_t));
- ctx->xmin = 0.0;
- ctx->xmax = 0.0;
- ctx->eps_x = 1.0e-8;
- ctx->eps_y = 1.0e-8;
- ctx->maxpass = (unsigned long)(DK4_ITER_PASSES_REGULAR);
- ctx->exact = 0;
- ctx->algo = DK4_ITER_ALG_RF_ANDERSON_BJOERCK;
- ctx->minmax = 0;
- }
-}
-
-
-
-dk4_iter_ctx_t *
-dk4iter_ctx_open(void)
-{
- dk4_iter_ctx_t *back = NULL;
- back = dk4mem_new(dk4_iter_ctx_t,1,NULL);
- if (NULL != back) {
- dk4iter_ctx_init(back);
- }
- return back;
-}
-
-
-
-void
-dk4iter_ctx_close(dk4_iter_ctx_t *ctx)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- dk4mem_free(ctx);
- }
-}
-
-
-
-void
-dk4iter_ctx_set_eps_x(dk4_iter_ctx_t *ctx, double eps)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->eps_x = eps;
- }
-}
-
-
-
-void
-dk4iter_ctx_set_eps_y(dk4_iter_ctx_t *ctx, double eps)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->eps_y = eps;
- }
-}
-
-
-
-void
-dk4iter_ctx_set_maxpass(dk4_iter_ctx_t *ctx, unsigned long passes)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->maxpass = passes;
- }
-}
-
-
-
-void
-dk4iter_ctx_set_exact(dk4_iter_ctx_t *ctx, int flag)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->exact = flag;
- }
-}
-
-
-
-void
-dk4iter_ctx_set_algorithm(dk4_iter_ctx_t *ctx, int algorithm)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->algo = algorithm;
- }
-}
-
-
-
-void
-dk4iter_ctx_set_min(dk4_iter_ctx_t *ctx, double xmin)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->xmin = xmin;
- ctx->minmax |= MINMAX_MINIMUM;
- }
-}
-
-
-
-void
-dk4iter_ctx_set_max(dk4_iter_ctx_t *ctx, double xmax)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->xmax = xmax;
- ctx->minmax |= MINMAX_MAXIMUM;
- }
-}
-
-
-
-void
-dk4iter_ctx_set_exclusive_min(dk4_iter_ctx_t *ctx, double xmin)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->xmin = xmin;
- ctx->minmax |= (MINMAX_MINIMUM | MINMAX_MIN_EXCL);
- }
-}
-
-
-
-void
-dk4iter_ctx_set_exclusive_max(dk4_iter_ctx_t *ctx, double xmax)
-{
-#if DK4_USE_ASSERT
- assert(NULL != ctx);
-#endif
- if (NULL != ctx) {
- ctx->xmax = xmax;
- ctx->minmax |= (MINMAX_MAXIMUM | MINMAX_MAX_EXCL);
- }
-}
-
-
-
-static
-int
-find_condition(double fa, double fb)
-{
- int back = CONDITION_ILLEGAL;
-
- if ((0.0 > fa) && (0.0 <= fb)) {
- back = CONDITION_FA_LESS_ZERO;
- }
- else {
- if ((0.0 < fa) && (0.0 >= fb)) {
- back = CONDITION_FA_GREATER_ZERO;
- }
- else {
- if ((0.0 >= fa) && (0.0 < fb)) {
- back = CONDITION_FA_LEQ_ZERO;
- }
- else {
- if ((0.0 <= fa) && (0.0 > fb)) {
- back = CONDITION_FA_GEQ_ZERO;
- }
- }
- }
- }
- return back;
-}
-
-
-
-/** Calculate gamma value for PEGASUS variant.
- @param fg Function value at border last set.
- @param fx Function value at recent x position.
- @return Gamma value.
-*/
-static
-double
-dk4iter_gamma_pegasus(double fg, double fx)
-{
- double back;
- back = fg / (fg + fx);
- if (0 == dk4ma_is_finite(back)) {
- back = 0.5;
- }
- else {
- if (0.0 >= back) {
- back = 0.5;
- }
-#if 0
- /* 2018-07-06
- Can not happen, fg and fx have same sign.
- */
- else {
- if (1.0 < back) {
- back = 1.0;
- }
- }
-#endif
- }
- return back;
-}
-
-
-
-/** Calculate gamma value for ANDERSON-BJOERCK variant.
- @param fg Function value at border last set.
- @param fx Function value at recent x position.
- @return Gamma value.
-*/
-static
-double
-dk4iter_gamma_anderson_bjoerck(double fg, double fx)
-{
- double back;
- back = 1.0 - fx / fg;
- if (0 == dk4ma_is_finite(back)) {
- back = 0.5;
- }
- else {
- if (0.0 >= back) {
- back = 0.5;
- }
-#if 0
- /* 2018-07-06
- Can not happen as fx and fy have same sign.
- */
- else {
- if (1.0 < back) {
- back = 1.0;
- }
- }
-#endif
- }
- return back;
-}
-
-
-
-int
-dk4iter_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; /* Copy of context */
- double fa = 0.0; /* Function value for border a */
- double fb = 0.0; /* Function value for border b */
- double x = 0.0; /* Text x value */
- double fx = 0.0; /* Function value for x */
- double afx = 0.0; /* Absolute function value for x */
- double xo = 0.0; /* Previous step x value */
- double gamma = 0.5; /* Correction factor */
- unsigned long passno = 0UL; /* Number of current pass */
- int res = 0; /* Operation result */
- int cond = 0; /* Condition of interval borders */
- int cc = 0; /* 1=continue, 0=finished, -1=abort */
- int pc = CHANGED_NONE; /* Previous border change */
- int nc = CHANGED_NONE; /* Next change */
- int back = DK4_ITER_RESULT_E_ARGS;
-
-
-#line 801 "dk4iter.ctr"
- /* Check function call arguments
- */
-#if DK4_USE_ASSERT
- assert(NULL != d);
- assert(NULL != fct);
-#endif
- if ((NULL == d) || (NULL == fct)) {
-
-#line 808 "dk4iter.ctr"
- goto finished;
- }
- if (a == b) {
-
-#line 811 "dk4iter.ctr"
- goto finished;
- }
-
- /* Copy or initialize context
- */
- if (NULL != ctx) {
-
-#line 817 "dk4iter.ctr"
- DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t));
- }
- else {
-
-#line 820 "dk4iter.ctr"
- dk4iter_ctx_init(&mctx);
- }
-
- /* Check algorithm specification
- */
- if (DK4_ITER_ALG_BISECTION != mctx.algo) {
- if (DK4_ITER_ALG_RF_PRIMITIVE != mctx.algo) {
- if (DK4_ITER_ALG_RF_ILLINOIS != mctx.algo) {
- if (DK4_ITER_ALG_RF_PEGASUS != mctx.algo) {
- if (DK4_ITER_ALG_RF_ANDERSON_BJOERCK != mctx.algo) {
- goto finished;
-
-#line 831 "dk4iter.ctr"
- }
- }
- }
- }
- }
-
- /* Calculate borders at beginning, check
- */
- res = (*fct)(&fa, a, ps);
- if ((0 == res) || (!(dk4ma_is_finite(fa)))) {
- back = DK4_ITER_RESULT_E_FCT;
-
-#line 842 "dk4iter.ctr"
- goto finished;
- }
- res = (*fct)(&fb, b, ps);
- if ((0 == res) || (!(dk4ma_is_finite(fb)))) {
- back = DK4_ITER_RESULT_E_FCT;
-
-#line 847 "dk4iter.ctr"
- goto finished;
- }
- cond = find_condition(fa, fb);
- if (CONDITION_ILLEGAL == cond) {
-
-#line 851 "dk4iter.ctr"
- goto finished;
- }
-
- /* Prepare iteration loop
- */
- passno = 0UL;
- x = xo = 0.0;
- cc = 1;
- pc = CHANGED_NONE;
-
- /* Run iteration loop
- */
- do {
- /* Keep previous x position
- */
- xo = x;
- /* Increase pass number, but avoid wrapping
- */
- if (ULONG_MAX > passno) { passno++; }
-
-#line 870 "dk4iter.ctr"
- /* Calculate x position
- */
-
-
-#line 873 "dk4iter.ctr"
- switch (mctx.algo) {
- case DK4_ITER_ALG_RF_PRIMITIVE :
- case DK4_ITER_ALG_RF_ILLINOIS :
- case DK4_ITER_ALG_RF_PEGASUS :
- case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : { /* Regula falsi */
- x = (a * fb - b * fa) / (fb - fa);
- } break;
- default : { /* Bisection */
- x = (a + b) / 2.0;
- } break;
- }
- if (0 == dk4ma_is_finite(x)) {
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
- }
- else {
-
-#line 889 "dk4iter.ctr"
- /* Calculate function for x, if x is usable
- */
- res = (*fct)(&fx, x, ps);
- if (0 == res) {
-
-#line 893 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_FCT;
- cc = -1;
- }
- if (0 == dk4ma_is_finite(fx)) {
-
-#line 897 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
- }
- afx = fabs(fx);
- if (0 == dk4ma_is_finite(afx)) {
-
-#line 902 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
- }
- /* Check whether we are done, if we could continue
- */
- if (1 == cc) {
-
-#line 908 "dk4iter.ctr"
-
-
-#line 909 "dk4iter.ctr"
-
-
-#line 910 "dk4iter.ctr"
-
-
-#line 911 "dk4iter.ctr"
- if (isgreater(mctx.eps_y,0.0)) {
- /* Must check y */
- if (afx < mctx.eps_y) {
-
-#line 914 "dk4iter.ctr"
- if ((x == xo) && (1UL < passno)) {
- cc = 0;
-
-#line 916 "dk4iter.ctr"
- }
- else {
- if (0 == mctx.exact) {
- if(isgreater(mctx.eps_x,0.0)) {
- if (isless(fabs(x - xo),mctx.eps_x)) {
- if (1UL < passno) {
- cc = 0;
-
-#line 923 "dk4iter.ctr"
- }
- }
- }
- else {
- cc = 0;
-
-#line 928 "dk4iter.ctr"
- }
- }
- }
- }
- }
- else {
-
-#line 934 "dk4iter.ctr"
- /* No y check */
- if ((x == xo) && (1UL < passno)) {
- cc = 0;
-
-#line 937 "dk4iter.ctr"
- }
- else {
- if (0 == mctx.exact) {
- if(isgreater(mctx.eps_x,0.0)) {
- if (isless(fabs(x - xo),mctx.eps_x)) {
- if (1UL < passno) {
- cc = 0;
-
-#line 944 "dk4iter.ctr"
- }
- }
- }
- else {
- /* Neither eps_x nor eps_y,
- so we wait for x=xo.
- */
- }
- }
- }
- }
- }
- /* Continue for usable y values only
- */
- if (1 == cc) {
- /* Find direction for next change
- */
- nc = CHANGED_NONE;
- switch (cond) {
- case CONDITION_FA_LESS_ZERO : {
- if (0.0 > fx) {
- nc = CHANGED_A;
- }
- else {
- nc = CHANGED_B;
- }
- } break;
- case CONDITION_FA_LEQ_ZERO : {
- if (0.0 >= fx) {
- nc = CHANGED_A;
- }
- else {
- nc = CHANGED_B;
- }
- } break;
- case CONDITION_FA_GREATER_ZERO : {
- if (0.0 < fx) {
- nc = CHANGED_A;
- }
- else {
- nc = CHANGED_B;
- }
- } break;
- case CONDITION_FA_GEQ_ZERO : {
- if (0.0 <= fx) {
- nc = CHANGED_A;
- }
- else {
- nc = CHANGED_B;
- }
- } break;
- }
- /* Apply interval border change
- */
- switch (nc) {
- case CHANGED_A: {
-
-#line 1000 "dk4iter.ctr"
- if (CHANGED_A == pc) {
- switch (mctx.algo) {
- case DK4_ITER_ALG_RF_ILLINOIS : {
- fb = 0.5 * fb;
- if (0 == dk4ma_is_finite(fb)) {
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
-
-#line 1007 "dk4iter.ctr"
- }
- } break;
- case DK4_ITER_ALG_RF_PEGASUS : {
- gamma = dk4iter_gamma_pegasus(fa, fx);
- fb *= gamma;
- if (0 == dk4ma_is_finite(fb)) {
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
-
-#line 1015 "dk4iter.ctr"
- }
- } break;
- case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : {
- gamma = dk4iter_gamma_anderson_bjoerck(
- fa, fx
- );
- fb *= gamma;
- if (0 == dk4ma_is_finite(fb)) {
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
-
-#line 1025 "dk4iter.ctr"
- }
- } break;
- }
- }
- a = x;
- fa = fx;
- pc = CHANGED_A;
- } break;
- case CHANGED_B: {
-
-#line 1034 "dk4iter.ctr"
- if (CHANGED_B == pc) {
- switch (mctx.algo) {
- case DK4_ITER_ALG_RF_ILLINOIS : {
- fa = 0.5 * fa;
- if (0 == dk4ma_is_finite(fa)) {
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
-
-#line 1041 "dk4iter.ctr"
- }
- } break;
- case DK4_ITER_ALG_RF_PEGASUS : {
- gamma = dk4iter_gamma_pegasus(
- fb, fx
- );
- fa *= gamma;
- if (0 == dk4ma_is_finite(fa)) {
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
-
-#line 1051 "dk4iter.ctr"
- }
- } break;
- case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : {
- gamma = dk4iter_gamma_anderson_bjoerck(
- fb, fx
- );
- fa *= gamma;
- if (0 == dk4ma_is_finite(fa)) {
- back = DK4_ITER_RESULT_E_INFINITE;
- cc = -1;
-
-#line 1061 "dk4iter.ctr"
- }
- } break;
- }
- }
- b = x;
- fb = fx;
- pc = CHANGED_B;
- } break;
- default: {
- /* ERROR: Must not happen */
- back = DK4_ITER_RESULT_E_CONV;
-
-#line 1072 "dk4iter.ctr"
- cc = -1;
- } break;
- }
- }
- }
- /* Check number of passes
- */
- if ((1 == cc) && (0UL < mctx.maxpass) && (passno >= mctx.maxpass)) {
- back = DK4_ITER_RESULT_E_PASSES;
-
-#line 1081 "dk4iter.ctr"
- cc = -1;
- }
- } while (1 == cc);
-
- /* Success
- */
- if (0 == cc) {
- *d = x;
- if (NULL != pp) { *pp = passno; }
- back = DK4_ITER_RESULT_SUCCESS;
- }
-
- finished:
-
-
-#line 1095 "dk4iter.ctr"
- return back;
-}
-
-
-
-/** Run Newton iteration.
- @param d Destination (address of result variable).
- @param pp Address of variable to store number of passes on success.
- @param fct Iteration function, returns two values into the array
- at address d: the function value and the first derivative
- value.
- @param ps Parameter set, may be NULL if fct does not use it.
- @param x0 Start point.
- @param ctx Iteration context, may be NULL.
- @return DK4_ITER_RESULT_SUCCESS on success, one from
- DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
- DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV, DK4_ITER_RESULT_E_FCT,
- or DK4_ITER_RESULT_E_ARGS on error.
-*/
-static
-int
-dk4iter_newton(
- double *d,
- unsigned long *pp,
- dk4_iter_fct_t *fct,
- void const *ps,
- double x0,
- dk4_iter_ctx_t const *ctx
-)
-{
- double v[2];
- double xn = 0.0;
- unsigned long passno = 0UL;
- int back = DK4_ITER_RESULT_E_ARGS;
- int cc = 1;
- int res = 0;
-
-
-#line 1132 "dk4iter.ctr"
-#if DK4_USE_ASSERT
- assert(NULL != d);
- assert(NULL != fct);
-#endif
- res = (*fct)(v, x0, ps);
- if (0 == res) {
- cc = -1;
-
-#line 1139 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_FCT;
- }
- if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) {
- cc = -1;
-
-#line 1143 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_INFINITE;
- }
- while (1 == cc) {
- if (ULONG_MAX > passno) { passno++; }
-
-#line 1147 "dk4iter.ctr"
- xn = x0 - v[0] / v[1];
- if (dk4ma_is_finite(xn)) {
-
-#line 1149 "dk4iter.ctr"
- if (0 != (MINMAX_MINIMUM & (ctx->minmax))) {
- if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) {
- if (xn <= ctx->xmin) {
- cc = -1;
-
-#line 1153 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- else {
- if (xn < ctx->xmin) {
- cc = -1;
-
-#line 1159 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- }
- if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) {
- if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) {
- if (xn >= ctx->xmax) {
- cc = -1;
-
-#line 1167 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- else {
- if (xn > ctx->xmax) {
- cc = -1;
-
-#line 1173 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- }
- if (1 == cc) {
- res = (*fct)(v, xn, ps);
- if (0 == res) {
- cc = -1;
-
-#line 1181 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_FCT;
- }
- if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) {
- cc = -1;
-
-#line 1185 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_INFINITE;
- }
- if (1 == cc) {
-
-
-#line 1189 "dk4iter.ctr"
-
-
-#line 1190 "dk4iter.ctr"
- if (isgreater(ctx->eps_y,0.0)) {
-
-#line 1191 "dk4iter.ctr"
- if (fabs(v[0]) < ctx->eps_y) {
-
-#line 1192 "dk4iter.ctr"
- if (xn == x0) {
- cc = 0;
-
-#line 1194 "dk4iter.ctr"
- }
- else {
- if (0 == ctx->exact) {
-
-#line 1197 "dk4iter.ctr"
- if(isgreater(ctx->eps_x,0.0)) {
- if (isless(fabs(xn - x0),ctx->eps_x)) {
- cc = 0;
-
-#line 1200 "dk4iter.ctr"
- }
- }
- else {
- cc = 0;
-
-#line 1204 "dk4iter.ctr"
- }
- }
- }
- }
- }
- else {
-
-#line 1210 "dk4iter.ctr"
- if (xn == x0) {
- cc = 0;
-
-#line 1212 "dk4iter.ctr"
- }
- else {
- if (0 == ctx->exact) {
- if (isgreater(ctx->eps_x,0.0)) {
- if (isless(fabs(xn - x0),ctx->eps_x)) {
- cc = 0;
-
-#line 1218 "dk4iter.ctr"
- }
- }
- else {
- /* Neither eps_x nor eps_y,
- so we wait for x=x0.
- */
- }
- }
- }
- }
- }
- x0 = xn;
- }
- }
- else {
- cc = -1;
-
-#line 1234 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_INFINITE;
- }
- if ((1 == cc) && (0UL < ctx->maxpass) && (passno >= ctx->maxpass)) {
- cc = -1;
- back = DK4_ITER_RESULT_E_PASSES;
- }
- }
- if (0 == cc) {
- *d = xn;
- if (NULL != pp) { *pp = passno; }
- back = DK4_ITER_RESULT_SUCCESS;
- }
-
-
-#line 1247 "dk4iter.ctr"
- return back;
-}
-
-
-
-/** Run fix point iteration.
- @param d Destination (address of result variable).
- @param pp Address of variable to store number of passes on success.
- @param fct Iteration function, the phi part of phi(x)=x.
- @param ps Parameter set, may be NULL if fct does not use it.
- @param x0 Start point.
- @param ctx Iteration context, may be NULL.
- @return DK4_ITER_RESULT_SUCCESS on success, one from
- DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
- DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV, DK4_ITER_RESULT_E_FCT,
- or DK4_ITER_RESULT_E_ARGS on error.
-*/
-static
-int
-dk4iter_fix_point(
- double *d,
- unsigned long *pp,
- dk4_iter_fct_t *fct,
- void const *ps,
- double x0,
- dk4_iter_ctx_t const *ctx
-)
-{
-#if 0
- double a; /* Border from previous steps */
- double b; /* Border from previous steps */
-#endif
- double xn = 0.0; /* New x value */
- unsigned long passno = 0UL; /* Iteration step number */
- int res = 0; /* Function evaluation result */
- int cc = 1; /* Flag: Can continue */
- int back = DK4_ITER_RESULT_E_ARGS;
-
-
-#line 1285 "dk4iter.ctr"
-#if DK4_USE_ASSERT
- assert(NULL != d);
- assert(NULL != fct);
-#endif
- while (1 == cc) {
- if (ULONG_MAX > passno) { passno++; }
-
-#line 1291 "dk4iter.ctr"
- /*
- Calculate new x, check result
- */
- res = (*fct)(&xn, x0, ps);
- if (0 == res) {
- cc = -1;
-
-#line 1297 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_FCT;
- }
- if (!(dk4ma_is_finite(xn))) {
- cc = -1;
-
-#line 1301 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_INFINITE;
- }
- if (1 == cc) {
-
-#line 1304 "dk4iter.ctr"
- /*
- Check whether specified interval is exceeded
- */
- if (0 != (MINMAX_MINIMUM & (ctx->minmax))) {
- if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) {
- if (xn <= ctx->xmin) {
- cc = -1;
-
-#line 1311 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- else {
- if (xn < ctx->xmin) {
- cc = -1;
-
-#line 1317 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- }
- if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) {
- if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) {
- if (xn >= ctx->xmax) {
- cc = -1;
-
-#line 1325 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- else {
- if (xn > ctx->xmax) {
- cc = -1;
-
-#line 1331 "dk4iter.ctr"
- back = DK4_ITER_RESULT_E_OOR;
- }
- }
- }
- if (1 == cc) {
- if (1 == cc) {
- /*
- Check whether we are finished
- */
- if (xn == x0) {
- cc = 0;
-
-#line 1342 "dk4iter.ctr"
- }
- else {
- if ((0 == ctx->exact) && (0.0 < ctx->eps_x)) {
- if (isless(fabs(xn-x0),ctx->eps_x)) {
- cc = 0;
-
-#line 1347 "dk4iter.ctr"
- }
-#if TRACE_DEBUG
- else {
-
-#line 1350 "dk4iter.ctr"
- }
-#endif
- }
- else {
-
-#line 1354 "dk4iter.ctr"
- }
- }
- }
- }
- x0 = xn;
- }
- /* Stop if too many passes
- */
- if (1 == cc) {
- if ((0UL < ctx->maxpass) && (passno >= ctx->maxpass)) {
- cc = -1;
- back = DK4_ITER_RESULT_E_PASSES;
- }
- }
- }
- if (0 == cc) {
- *d = xn;
- if (NULL != pp) { *pp = passno; }
- back = DK4_ITER_RESULT_SUCCESS;
- }
-
-
-#line 1375 "dk4iter.ctr"
- return back;
-}
-
-
-
-int
-dk4iter_start_point(
- double *d,
- unsigned long *pp,
- dk4_iter_fct_t *fct,
- void const *ps,
- double x0,
- dk4_iter_ctx_t const *ctx
-)
-{
- dk4_iter_ctx_t mctx;
- int back = DK4_ITER_RESULT_E_ARGS;
-
- /* Check function call arguments
- */
-#if DK4_USE_ASSERT
- assert(NULL != d);
- assert(NULL != fct);
-#endif
- if ((NULL == d) || (NULL == fct)) {
- goto finished;
- }
- /* Copy or initialize context
- */
- if (NULL != ctx) {
- DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t));
- }
- else {
- dk4iter_ctx_init(&mctx);
- mctx.algo = DK4_ITER_ALG_NEWTON;
- }
- /* Check x0 in interval
- */
- if (0 != (MINMAX_MINIMUM & (mctx.minmax))) {
- if (0 != (MINMAX_MIN_EXCL & (mctx.minmax))) {
- if (x0 <= mctx.xmin) {
- goto finished;
- }
- }
- else {
- if (x0 < mctx.xmin) {
- goto finished;
- }
- }
- }
- if (0 != (MINMAX_MAXIMUM & (mctx.minmax))) {
- if (0 != (MINMAX_MAX_EXCL & (mctx.minmax))) {
- if (x0 >= mctx.xmax) {
- goto finished;
- }
- }
- else {
- if (x0 > mctx.xmax) {
- goto finished;
- }
- }
- }
- /* Call function for specified algorithm.
- */
- switch (mctx.algo) {
- case DK4_ITER_ALG_NEWTON : {
- back = dk4iter_newton(d, pp, fct, ps, x0, &mctx);
- } break;
- case DK4_ITER_ALG_FIX_POINT : {
- back = dk4iter_fix_point(d, pp, fct, ps, x0, &mctx);
- } break;
- }
-
- finished:
- return back;
-}
-
-
-
-
-/* vim: set ai sw=4 ts=4 : */
-