From 004ea10198960ce0f7395975b027fd112edf960c Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 2 Jan 2020 14:13:23 +0100 Subject: tmp --- ell_grid_solve.c | 260 ++++++++++++++++++++------------------------------ ell_grid_solve.h | 2 + mg2d.c | 114 +++++++++++++++++----- mg2d.h | 7 ++ step_control.h | 283 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 483 insertions(+), 183 deletions(-) create mode 100644 step_control.h diff --git a/ell_grid_solve.c b/ell_grid_solve.c index ecbf6dc..daa62e6 100644 --- a/ell_grid_solve.c +++ b/ell_grid_solve.c @@ -39,92 +39,9 @@ #include "mg2d_constants.h" #include "ndarray.h" #include "residual_calc.h" +#include "step_control.h" #include "timer.h" -typedef struct StepControl { - /* estimated value for a good step */ - double hint; - /* minimum step size, fail if this value is reached */ - double step_min; - - double cur_converge; - double cur_diverge; - double prev_converge; - double prev_diverge; - - int64_t diverge_count; -} StepControl; - -static void sc_reset(StepControl *sc) -{ - sc->hint = DBL_MAX; - sc->cur_converge = DBL_MAX; - sc->cur_diverge = DBL_MAX; - sc->prev_converge = DBL_MAX; - sc->prev_diverge = DBL_MAX; -} - -static void sc_init(StepControl *sc, double hint, double step_min) -{ - sc->hint = hint; - sc->step_min = step_min; - sc->cur_converge = DBL_MAX; - sc->cur_diverge = DBL_MAX; -} - -static void sc_step_converge(StepControl *sc, double step, double conv_fact) -{ - sc->cur_converge = step; - sc->prev_converge = step; - mg2di_assert(sc->cur_diverge == DBL_MAX || step < sc->cur_diverge); - if (sc->prev_diverge != DBL_MAX && step >= sc->prev_diverge) - sc->prev_diverge = DBL_MAX; -} - -static void sc_step_diverge(StepControl *sc, double step) -{ - sc->cur_diverge = step; - sc->prev_diverge = step; - if (sc->cur_converge != DBL_MAX && step <= sc->cur_converge) - sc->cur_converge = DBL_MAX; - if (sc->prev_converge != DBL_MAX && step <= sc->prev_converge) - sc->prev_converge = DBL_MAX; - - sc->diverge_count++; -} - -static double sc_step_get(StepControl *sc) -{ - double step = 0.25; - double guess_converge = sc->cur_converge == DBL_MAX ? sc->prev_converge : sc->cur_converge; - double guess_diverge = sc->cur_diverge == DBL_MAX ? sc->prev_diverge : sc->cur_diverge; - - if (guess_converge != DBL_MAX && guess_diverge != DBL_MAX) { - if (guess_diverge - guess_converge < 0.1 * guess_diverge) - step = guess_converge; - else - step = 0.5 * (guess_diverge + guess_converge); - } else if (guess_diverge != DBL_MAX) { - step = 0.8 * guess_diverge; - } else if (guess_converge != DBL_MAX) { - step = 1.1 * guess_converge; - } else if (sc->hint != DBL_MAX) - step = sc->hint; - - if (sc->step_min != DBL_MAX && step < sc->step_min) { - fprintf(stderr, "%g %g %g %g %g\n", sc->cur_converge, sc->prev_converge, sc->cur_diverge, sc->prev_diverge, sc->hint); - mg2di_assert(0); - return -1.0; - } - - return step; -} - -static int sc_step_confirm(StepControl *sc, double conv_factor) -{ - return 0; -} - static const double relax_factors[FD_STENCIL_MAX] = { [0] = 1.0 / 5, [1] = 1.0 / 7, @@ -500,18 +417,21 @@ static int solve_relax_step(EGSContext *ctx, int export_res) { EGSInternal *priv = ctx->priv; EGSRelaxContext *r = ctx->relax; - int ic = priv->iteration_cur; + + const int it_cur = priv->iteration_cur % HIST_SIZE; + const int it_next = (priv->iteration_cur + 1) % HIST_SIZE; + const int it_next1 = (priv->iteration_cur + 2) % HIST_SIZE; if (export_res) { - if (priv->u_iteration[(ic + 1) % HIST_SIZE] < ic + 1) { - if (priv->residual_iteration < ic) { - residual_calc(ctx, &priv->residual_max[ic % HIST_SIZE], NULL, priv->u_interior[ic % HIST_SIZE], 0); - priv->residual_iteration = ic; + if (priv->u_iteration[it_next] < priv->iteration_cur + 1) { + if (priv->residual_iteration < priv->iteration_cur) { + residual_calc(ctx, &priv->residual_max[it_cur], NULL, priv->u_interior[it_cur], 0); + priv->residual_iteration = priv->iteration_cur; } - mg2di_assert(priv->residual_iteration == ic); + mg2di_assert(priv->residual_iteration == priv->iteration_cur); - priv->step[ic % HIST_SIZE] = sc_step_get(&priv->sc); - if (priv->step[ic % HIST_SIZE] < 0) { + priv->step[it_cur] = sc_step_get(&priv->sc); + if (priv->step[it_cur] < 0) { mg2di_log(&ctx->logger, MG2D_LOG_ERROR, "Minimum step reached."); return -EDOM; @@ -525,8 +445,8 @@ static int solve_relax_step(EGSContext *ctx, int export_res) priv->reflect_skip = 0; - boundaries_apply(ctx, priv->u_interior[(ic + 1) % HIST_SIZE], 0); - priv->u_iteration[(ic + 1) % HIST_SIZE] = ic + 1; + boundaries_apply(ctx, priv->u_interior[it_next], 0); + priv->u_iteration[it_next] = priv->iteration_cur + 1; priv->steps_since_sync++; if (priv->steps_since_sync >= priv->steps_per_sync) { @@ -534,80 +454,95 @@ static int solve_relax_step(EGSContext *ctx, int export_res) priv->steps_since_sync = 0; } } - mg2di_assert(priv->u_iteration[(ic + 1) % HIST_SIZE] == ic + 1); + mg2di_assert(priv->u_iteration[it_next] == priv->iteration_cur + 1); - residual_calc(ctx, &priv->residual_max[(ic + 1) % HIST_SIZE], NULL, - priv->u_interior[(ic + 1) % HIST_SIZE], 0); + residual_calc(ctx, &priv->residual_max[it_next], NULL, + priv->u_interior[it_next], 0); boundaries_sync(ctx, ctx->residual); - priv->residual_iteration = ic + 1; + priv->residual_iteration = priv->iteration_cur + 1; } else { - double res_old, res_new; - int diverged; - - mg2di_assert(priv->u_iteration[(ic + 1) % HIST_SIZE] == ic + 1); - res_old = priv->residual_max[ic % HIST_SIZE]; - - priv->step[(ic + 1) % HIST_SIZE] = sc_step_get(&priv->sc); - if (priv->step[(ic + 1) % HIST_SIZE] < 0) { - mg2di_log(&ctx->logger, MG2D_LOG_ERROR, - "Minimum step reached."); - return -EDOM; + double step_next; + + mg2di_assert(priv->u_iteration[it_next] == priv->iteration_cur + 1); + + /* Adaptive relaxation stepping - this is slightly tricky, since + * residual computation and computing the next iteration are + * merged into a single step. + * + * At this point, we have already computed the next solution iteration + * -- u[it_next], but don't have the residual norm for it, so we haven't + * verified whether the computation converged. We now generate a trial + * next timestep and use it to compute the next residual+norm and + * next+1 solution. + */ + step_next = r->adaptive_step ? sc_step_get(&priv->sc) : + r->relax_factor * r->relax_multiplier;//priv->step[0]; + priv->r.relax_factor = step_next; + //fprintf(stderr, "%zu step next %g\n", ctx->domain_size[0], step_next); + + residual_calc(ctx, &priv->residual_max[it_next], priv->u_interior[it_next1], + priv->u_interior[it_next], step_next); + + if (r->adaptive_step) { + const double res_old = priv->residual_max[it_cur]; + int step_reject; + + /* Now that we have the next residual norm, we can verify whether + * the u[cur] --(step[cur])--> u[next] step was converging and + * accept or reject it. + * + * If it's rejected, we must generate new step[cur] and redo the + * calculations starting from u[cur]. + */ + do { + double res_new = priv->residual_max[it_next]; + double tmp; + + step_reject = sc_step_confirm(&priv->sc, priv->step[it_cur], res_old, res_new); + if (step_reject < 0) { + mg2di_log(&ctx->logger, MG2D_LOG_ERROR, "Error confirming step: %g", step_next); + return step_reject; + } else if (!step_reject) + break; + + priv->step[it_cur] = sc_step_get(&priv->sc); + residual_calc(ctx, &tmp, priv->u_interior[it_next], + priv->u_interior[it_cur], priv->step[it_cur]); + + step_next = sc_step_get(&priv->sc); + fprintf(stderr, "%zu retry step next %g %g\n", ctx->domain_size[0], step_next, res_old / res_new); + residual_calc(ctx, &priv->residual_max[it_next], priv->u_interior[it_next1], + priv->u_interior[it_next], step_next); + res_new = priv->residual_max[it_next]; + } while (step_reject); + + priv->step[it_next] = step_next; } - residual_calc(ctx, &priv->residual_max[(ic + 1) % HIST_SIZE], priv->u_interior[(ic + 2) % HIST_SIZE], - priv->u_interior[(ic + 1) % HIST_SIZE], priv->step[(ic + 1) % HIST_SIZE]); - - res_new = priv->residual_max[(ic + 1) % HIST_SIZE]; - diverged = res_new > 1e-6 && res_new / res_old - 1.0 >= 1e-2 ; - while (diverged) { - double tmp, step; - - fprintf(stderr, "%zu diverge step %g fact %g\n", ctx->domain_size[0], - priv->step[(ic + 1) % HIST_SIZE], res_new / res_old); - - sc_step_diverge(&priv->sc, priv->step[ic % HIST_SIZE]); - - priv->step[ic % HIST_SIZE] = sc_step_get(&priv->sc); - if (priv->step[(ic) % HIST_SIZE] < 0) { - mg2di_log(&ctx->logger, MG2D_LOG_ERROR, - "Minimum step reached."); - return -EDOM; - } - residual_calc(ctx, &tmp, priv->u_interior[(ic + 1) % HIST_SIZE], - priv->u_interior[ic % HIST_SIZE], priv->step[ic % HIST_SIZE]); - - priv->step[(ic + 1) % HIST_SIZE] = priv->step[ic % HIST_SIZE]; - residual_calc(ctx, &priv->residual_max[(ic + 1) % HIST_SIZE], priv->u_interior[(ic + 2) % HIST_SIZE], - priv->u_interior[(ic + 1) % HIST_SIZE], priv->step[(ic + 1) % HIST_SIZE]); - res_new = priv->residual_max[(ic + 1) % HIST_SIZE]; - diverged = res_new > 1e-6 && res_new >= res_old; - } - - sc_step_converge(&priv->sc, priv->step[ic % HIST_SIZE], res_old / res_new); - boundaries_apply(ctx, priv->u_interior[(ic + 2) % HIST_SIZE], 0); - priv->u_iteration[(ic + 2) % HIST_SIZE] = ic + 2; + boundaries_apply(ctx, priv->u_interior[it_next1], 0); + priv->u_iteration[it_next1] = priv->iteration_cur + 2; priv->steps_since_sync++; if (priv->steps_since_sync >= priv->steps_per_sync) { - boundaries_sync(ctx, priv->u_interior[(ic + 2) % HIST_SIZE]); + boundaries_sync(ctx, priv->u_interior[it_next1]); priv->steps_since_sync = 0; } - { - static int i; - if (!(i & 1023)) { - fprintf(stderr, "%zu converge %g diverge %g guess %g diverges %ld\n", ctx->domain_size[0], priv->sc.cur_converge, - priv->sc.prev_diverge, ctx->relax->relax_factor * ctx->relax->relax_multiplier, priv->sc.diverge_count); - } - i++; - } + //if (r->adaptive_step) { + // static int i; + // if (!(i & 1023)) { + // fprintf(stderr, "%zu converge %g %g diverge %g guess %g diverges %ld\n", ctx->domain_size[0], priv->sc.cur_converge, + // priv->sc.cur_converge / ctx->relax->relax_multiplier, + // priv->sc.prev_diverge, ctx->relax->relax_factor * ctx->relax->relax_multiplier, priv->sc.diverge_count); + // } + // i++; + //} } - ic++; - ctx->u = priv->u_interior [ic % HIST_SIZE]; - ctx->u_exterior = priv->u_exterior [ic % HIST_SIZE]; - ctx->residual_max = priv->residual_max[ic % HIST_SIZE]; - priv->iteration_cur = ic; + priv->iteration_cur++; + ctx->u = priv->u_interior [(priv->iteration_cur) % HIST_SIZE]; + ctx->u_exterior = priv->u_exterior [(priv->iteration_cur) % HIST_SIZE]; + ctx->residual_max = priv->residual_max[(priv->iteration_cur) % HIST_SIZE]; return 0; } @@ -1482,11 +1417,16 @@ finish: } priv->steps_since_sync = 0; - priv->step[0] = sc_step_get(&priv->sc); - if (priv->step[0] < 0) { - mg2di_log(&ctx->logger, MG2D_LOG_ERROR, - "Minimum step reached. 0"); - return -EDOM; + if (r->adaptive_step) { + priv->step[0] = sc_step_get(&priv->sc); + if (priv->step[0] < 0) { + mg2di_log(&ctx->logger, MG2D_LOG_ERROR, + "Minimum step reached.\n"); + return -EDOM; + } + } else { + for (int i = 0; i < ARRAY_ELEMS(priv->step); i++) + priv->step[i] = r->relax_factor * r->relax_multiplier; } residual_calc(ctx, &priv->residual_max[0], priv->u_interior[1], priv->u_interior[0], priv->step[0]); diff --git a/ell_grid_solve.h b/ell_grid_solve.h index 095c591..27c2848 100644 --- a/ell_grid_solve.h +++ b/ell_grid_solve.h @@ -88,6 +88,8 @@ typedef struct EGSRelaxContext { */ double relax_multiplier; + int adaptive_step; + Timer timer_correct; } EGSRelaxContext; diff --git a/mg2d.c b/mg2d.c index 06e7b38..7d01761 100644 --- a/mg2d.c +++ b/mg2d.c @@ -37,6 +37,7 @@ #include "mg2d_boundary_internal.h" #include "mg2d_constants.h" #include "ndarray.h" +#include "step_control.h" #include "timer.h" #include "transfer.h" @@ -79,6 +80,8 @@ typedef struct MG2DLevel { struct MG2DLevel *child; + StepControl sc; + /* timings */ int64_t count_cycles; Timer timer_solve; @@ -228,6 +231,7 @@ static int mg_prolong_u(MG2DContext *ctx, MG2DLevel *l) return 0; } +#define SUBGRID_DIVERGE 1 static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact) { enum EGSType solve_type = (allow_exact && level->dg->nb_components == 1 && @@ -264,7 +268,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact) goto finish; } - for (int i = 0; i < ctx->nb_cycles; i++) { + for (int cycle = 0; cycle < ctx->nb_cycles; cycle++) { double res_prev; /* pre-restrict relaxation */ @@ -277,17 +281,19 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact) } if (level->depth < ctx->priv->dh.nb_levels - 1) { - /* restrict the residual as to the coarser-level rhs */ - ret = mg_restrict_rhs(ctx, level); - if (ret < 0) - return ret; - - /* solve on the coarser level */ - if (level->child) { - ret = mg_solve_subgrid(ctx, level->child, 1); + do { + /* restrict the residual as to the coarser-level rhs */ + ret = mg_restrict_rhs(ctx, level); if (ret < 0) return ret; - } + + /* solve on the coarser level */ + if (level->child) { + ret = mg_solve_subgrid(ctx, level->child, 1); + if (ret < 0) + return ret; + } + } while (ret == SUBGRID_DIVERGE); /* prolongate the coarser-level correction */ ret = mg_prolong_u(ctx, level); @@ -313,10 +319,59 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact) } /* post-correct relaxation */ - for (int j = 0; j < ctx->nb_relax_post; j++) { - ret = mg_solve(ctx, level, solve_type, "post-relax", 0); - if (ret < 0) - return ret; + { + double conv_history[4][2] = { { 0. } }; + int conv_hist_next = 0; + + double min_fact = DBL_MAX; + int min_fact_idx = -1; + + for (int j = 0; j < ctx->nb_relax_post; j++) { + conv_history[conv_hist_next][0] = level->solver->residual_max; + + ret = mg_solve(ctx, level, solve_type, "post-relax", 0); + if (ret < 0) + return ret; + + conv_history[conv_hist_next][1] = level->solver->residual_max; + conv_hist_next = (conv_hist_next + 1) % ARRAY_ELEMS(conv_history); + } + + // find the minimum convergence factor in the history window + for (int i = 0; i < ARRAY_ELEMS(conv_history); i++) { + double *h = conv_history[i], f; + if (h[0] == 0.0) + continue; + f = h[0] / h[1]; + if (f < min_fact) { + min_fact = f; + min_fact_idx = i; + } + } + + // use the minimum convergence factor as input for step control + if (level->solver->residual_max > ctx->tol && min_fact_idx >= 0) { + double *h = conv_history[min_fact_idx]; + double step; + + ret = sc_step_confirm(&level->sc, level->solver->relax->relax_factor, + h[0], h[1]); + if (ret < 0) + return ret; + + step = sc_step_get(&level->sc); + if (fabs(step - level->solver->relax->relax_factor) > 1e-6) { + fprintf(stderr, "%d step %g -> %g\n", level->depth, + level->solver->relax->relax_factor, step); + } + if (ret) { + fprintf(stderr, "%d diverge %g -> %g\n", level->depth, + level->solver->relax->relax_factor, step); + level->solver->relax->relax_factor = step; + return SUBGRID_DIVERGE; + } + level->solver->relax->relax_factor = step; + } } level->count_cycles++; @@ -786,9 +841,9 @@ static int mg_levels_init(MG2DContext *ctx) ret = mg2di_ndarray_copy(cur->solver->u, priv->u); if (ret < 0) return ret; - mg2di_ndarray_free(&priv->u); - ctx->u = cur->solver->u->data; - ctx->u_stride = cur->solver->u->stride[0]; + //mg2di_ndarray_free(&priv->u); + //ctx->u = cur->solver->u->data; + //ctx->u_stride = cur->solver->u->stride[0]; } if (priv->rhs) { @@ -865,10 +920,10 @@ static int mg_levels_init(MG2DContext *ctx) if (ctx->fd_stencil == 1) op_restrict = GRID_TRANSFER_FW_2; else - op_restrict = GRID_TRANSFER_FW_3_GENERIC; + op_restrict = GRID_TRANSFER_FW_1; } else { op_restrict = GRID_TRANSFER_LAGRANGE_1; - op_restrict = GRID_TRANSFER_FW_3_GENERIC; + op_restrict = GRID_TRANSFER_FW_1_GENERIC; } ret = mg_setup_restrict(ctx, cur, op_restrict); @@ -904,9 +959,12 @@ static int mg_levels_init(MG2DContext *ctx) cur = priv->root; while (cur) { - cur->solver->relax->relax_factor = ctx->cfl_factor; + //cur->solver->relax->relax_factor = ctx->cfl_factor; + sc_init(&cur->sc, ctx->cfl_factor, ctx->cfl_factor * 1e-3); + cur->solver->relax->relax_factor = sc_step_get(&cur->sc); cur->solver->relax->relax_multiplier = 1.0 / (diff2_max + cur->solver->step[0] * cur->solver->step[1] * diff0_max / 8.0); + cur->solver->relax->adaptive_step = 0;//ctx->adaptive_step; if (cur->depth > 0) { ret = diff_coeffs_fixup(ctx, cur); @@ -1058,6 +1116,7 @@ static int mg_dh_init(DomainHierarchy *dh, const DomainGeometry *root, dh->stepsizes[depth][0] = step_root[0] * (root->domain_size[0] - 1) / (cur_dg->domain_size[0] - 1); dh->stepsizes[depth][1] = step_root[1] * (root->domain_size[1] - 1) / (cur_dg->domain_size[1] - 1); + fprintf(stderr, "level %d %zux%zu\n", depth, cur_dg->domain_size[0], cur_dg->domain_size[1]); dh->nb_levels++; if (next_size <= 4) @@ -1127,9 +1186,15 @@ int mg2d_solve(MG2DContext *ctx) res_prev = res_orig; for (int i = 0; i < ctx->maxiter; i++) { - ret = mg_solve_subgrid(ctx, root, 1); - if (ret < 0) - goto fail; + do { + ret = mg_solve_subgrid(ctx, root, 1); + if (ret < 0) + goto fail; + if (ret == SUBGRID_DIVERGE) { + mg2di_ndarray_copy(s_root->u, priv->u); + mg2di_egs_init(s_root, 0); + } + } while (ret == SUBGRID_DIVERGE); res_cur = s_root->residual_max; @@ -1159,6 +1224,7 @@ int mg2d_solve(MG2DContext *ctx) fail: mg2di_log(&priv->logger, MG2D_LOG_ERROR, "The solver failed to converge\n"); finish: + ctx->residual_max = s_root->residual_max; mg2di_timer_stop(&priv->timer_solve); return ret; } @@ -1291,6 +1357,8 @@ static MG2DLevel *mg_level_alloc(const DomainGeometry *dg, mg2di_timer_init(&level->timer_reinit); mg2di_timer_init(&level->timer_mpi_sync); + sc_reset(&level->sc); + return level; fail: mg_level_free(&level); diff --git a/mg2d.h b/mg2d.h index f7897f6..0c12f22 100644 --- a/mg2d.h +++ b/mg2d.h @@ -237,6 +237,13 @@ typedef struct MG2DContext { * Set by mg2d_solver_alloc[_mpi](). */ size_t local_size[2]; + + double residual_max; + + /** + * Use adaptive stepping for relaxation. + */ + int adaptive_step; } MG2DContext; /** diff --git a/step_control.h b/step_control.h new file mode 100644 index 0000000..187c80b --- /dev/null +++ b/step_control.h @@ -0,0 +1,283 @@ +/* + * Copyright 2020 Anton Khirnov + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef MG2D_STEP_CONTROL_H +#define MG2D_STEP_CONTROL_H + +#include + +#include "common.h" + +#define SC_HIST_SIZE 16 + +typedef struct StepControl { + /* estimated value for a good step */ + double hint; + /* minimum step size, fail if this value is reached */ + double step_min; + + int steps_since_init; + + double step[SC_HIST_SIZE]; + double fact[SC_HIST_SIZE]; + + double conv_epsilon; + + int hist_len; + int idx_bracket; + int idx_diverge; + + int diverge_count; +} StepControl; + +static void sci_hist_clear(StepControl *sc) +{ + for (int i = 0; i < SC_HIST_SIZE; i++) { + sc->step[i] = DBL_MAX; + sc->fact[i] = DBL_MAX; + } + + sc->hist_len = 0; + sc->idx_bracket = -1; + sc->idx_diverge = -1; +} + +static void sc_reset(StepControl *sc) +{ + sci_hist_clear(sc); + sc->hint = DBL_MAX; + sc->step_min = DBL_MAX; + sc->steps_since_init = 0; +} + +static void sc_init(StepControl *sc, double hint, double step_min) +{ + sc->hint = hint; + sc->step_min = step_min; + sc->steps_since_init = 0; + sc->conv_epsilon = hint * 1e-3; +} + +static double sc_step_get(StepControl *sc) +{ + // no history present, try the hint + if (!sc->hist_len) + return (sc->hint == DBL_MAX) ? 0.25 : sc->hint; + + // don't have the bracket yet + if (sc->idx_bracket < 0) { + const double high_step = sc->step[sc->hist_len - 1]; + const double high_fact = sc->fact[sc->hist_len - 1]; + + // don't know where divergence happens + // -> try a higher step + if (high_fact > 1.0) + return high_step * 1.2; + + // convergence factor should decrease towards smaller timesteps + // so just try lower ones until we get a bracket + return MAX(sc->step_min, sc->step[0] * 0.8); + } + + // TODO periodic random steps to try escaping local minima + + // got a bracket, golden-section-search it until convergence + mg2di_assert(sc->hist_len - sc->idx_bracket >= 3); + mg2di_assert(sc->fact[sc->idx_bracket] < sc->fact[sc->idx_bracket + 1]); + mg2di_assert(sc->fact[sc->idx_bracket + 2] < sc->fact[sc->idx_bracket + 1]); + { + const double dist0 = sc->step[sc->idx_bracket + 1] - sc->step[sc->idx_bracket]; + const double dist1 = sc->step[sc->idx_bracket + 2] - sc->step[sc->idx_bracket + 1]; + + // converged + if (dist0 < sc->conv_epsilon && dist1 < sc->conv_epsilon) + return sc->step[sc->idx_bracket + 1]; + + if (dist0 > dist1) + return sc->step[sc->idx_bracket] + 0.61803 * dist0; + else + return sc->step[sc->idx_bracket + 2] - 0.61803 * dist1; + } +} + +static void sci_find_bracket(StepControl *sc) +{ + for (int i = sc->hist_len - 2; i >= 0; i--) { + if (sc->fact[i] < sc->fact[i + 1] && + sc->fact[i + 2] < sc->fact[i + 1]) { + sc->idx_bracket = i; + return; + } + } +} + +static int sc_step_confirm(StepControl *sc, const double step, + const double norm_old, const double norm_new) +{ + const double conv_fact = norm_old / norm_new; + + // signal failure if we reached minimum step + if (step <= sc->step_min) + return -EINVAL; + + // ignore the first relax step, it tends to be weird + sc->steps_since_init++; + //if (sc->steps_since_init < 2) + // return 0; + + // divergence + if (conv_fact <= 1.0) { + //if (step < sc->hint / 2.) + // return 0; + + // if the step is not larger then the largest known diverging step, + // add it into history + if (!sc->hist_len || + sc->idx_diverge < 0 || + sc->step[sc->idx_diverge] > step) { + int idx = sc->idx_diverge >= 0 ? sc->idx_diverge : sc->hist_len; + + // history size overflow, shouldn't happen + if (idx >= SC_HIST_SIZE) + return -ENOMEM; + + // largest known coverging step doesn't converge anymore, + // clear history + if (idx > 0 && sc->step[idx - 1] <= step) { + sci_hist_clear(sc); + idx = 0; + } + + sc->step[idx] = step; + sc->fact[idx] = conv_fact; + sc->idx_diverge = idx; + sc->hist_len = idx + 1; + + // check if we now have a bracket + if (sc->idx_bracket < 0 && sc->hist_len >= 3 && + sc->fact[idx - 1] > conv_fact && + sc->fact[idx - 1] > sc->fact[idx - 2]) { + sc->idx_bracket = idx - 2; + } + + } + + return 1; + } + + // convergence + + // history empty, add first element + if (!sc->hist_len) { + sc->step[0] = step; + sc->fact[0] = conv_fact; + sc->hist_len = 1; + return 0; + } + + // previously diverging step is now converging + // drop all diverging steps from history + if (sc->idx_diverge >= 0 && + step >= sc->step[sc->idx_diverge]) { + const int idx = sc->idx_diverge; + + sc->idx_diverge = -1; + sc->step[idx] = step; + sc->fact[idx] = conv_fact; + sc->hist_len = idx + 1; + + // check if we just broke bracketing + if (sc->idx_bracket >= 0 && sc->idx_bracket + 2 == idx && + sc->fact[sc->idx_bracket + 1] <= conv_fact) + sc->idx_bracket = -1; + + // could not have created a bracket where there was none before, since + // we increased the factor for topmost step + + return 0; + } + + // no bracket + if (sc->idx_bracket < 0) { + int idx_floor = -1; + + for (int i = 0; i < sc->hist_len; i++) { + if (sc->step[i] < step) + idx_floor = i; + else + break; + } + if (idx_floor + 1 < SC_HIST_SIZE) { + memmove(sc->step + idx_floor + 2, sc->step + idx_floor + 1, + sizeof(*sc->step) * (SC_HIST_SIZE - (idx_floor + 2))); + memmove(sc->fact + idx_floor + 2, sc->fact + idx_floor + 1, + sizeof(*sc->fact) * (SC_HIST_SIZE - (idx_floor + 2))); + + sc->step[idx_floor + 1] = step; + sc->fact[idx_floor + 1] = conv_fact; + sc->hist_len = MIN(SC_HIST_SIZE, sc->hist_len + 1); + if (sc->idx_diverge > idx_floor) + sc->idx_diverge = (sc->idx_diverge + 1 < SC_HIST_SIZE) ? sc->idx_diverge + 1 : -1; + } + sci_find_bracket(sc); + + return 0; + } + + // have bracket + + // got a step outside of the bracket for some reason + if (step <= sc->step[sc->idx_bracket] || + step >= sc->step[sc->idx_bracket + 2]) + return 0; + + if (step < sc->step[sc->idx_bracket + 1]) { + // step inside lower interval + if (conv_fact > sc->fact[sc->idx_bracket + 1]) { + // replace upper bound + sc->step[sc->idx_bracket + 2] = sc->step[sc->idx_bracket + 1]; + sc->fact[sc->idx_bracket + 2] = sc->fact[sc->idx_bracket + 1]; + if (sc->idx_diverge == sc->idx_bracket + 2) + sc->idx_diverge = -1; + + sc->step[sc->idx_bracket + 1] = step; + sc->fact[sc->idx_bracket + 1] = conv_fact; + } else { + // replace lower bound + sc->step[sc->idx_bracket] = step; + sc->fact[sc->idx_bracket] = conv_fact; + } + } else { + // step inside upper interval + if (conv_fact > sc->fact[sc->idx_bracket + 1]) { + // replace lower bound + sc->step[sc->idx_bracket] = sc->step[sc->idx_bracket + 1]; + sc->fact[sc->idx_bracket] = sc->fact[sc->idx_bracket + 1]; + sc->step[sc->idx_bracket + 1] = step; + sc->fact[sc->idx_bracket + 1] = conv_fact; + } else { + sc->step[sc->idx_bracket + 2] = step; + sc->fact[sc->idx_bracket + 2] = conv_fact; + if (sc->idx_diverge == sc->idx_bracket + 2) + sc->idx_diverge = -1; + } + } + + return 0; +} + +#endif // MG2D_STEP_CONTROL_H -- cgit v1.2.3