#include "multigrid.h"
/* if defined, residual is calculated and loop is aborted, if no improvement
- * is reached. this does not seem to work, so it should be undefined */
-//#define RESIDUAL
+ * is reached or the target resudual change is reached. */
+#define RESIDUAL
/* just process, but do not generate neighbor list and do not solve multigrid */
//#define TEST_NO_RENDER
* before z)!
*
*/
-int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **marks, const darray_t *init, unsigned char **nb_list, darray_t **flows, int scalexyz)
+int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **marks, const darray_t *init, unsigned char **nb_list, darray_t **flows, int scalexyz, double target_residual_change)
{
#ifdef RESIDUAL
double curr_residual, prev_residual = 10000000000000.0;
for (i = 0; i < iters; i++)
smooth(level, values[level], marks[level], init, nb_list[level]);
#ifdef RESIDUAL
- curr_residual += residual(level, values[level], marks[level], init, nb_list[level]);
+ if (target_residual_change)
+ curr_residual += residual(level, values[level], marks[level], init, nb_list[level]);
#endif
if (level < scalexyz)
rc = restrict(values[level], values[level+1], flows[level]);
/* end if residual raises again */
#ifdef RESIDUAL
- if (prev_residual < curr_residual)
+//printf("residual %.8f (change=%%%.4f)\n", curr_residual, (prev_residual-curr_residual)/prev_residual*100);
+ if (target_residual_change)
+ printf(" %.1f%%", (prev_residual-curr_residual)/prev_residual*100); fflush(stdout);
+ if (curr_residual == 0.0) {
+ printf(" 0"); fflush(stdout);
break;
- curr_residual = prev_residual;
+ }
+ if (target_residual_change && (prev_residual-curr_residual)/prev_residual < target_residual_change)
+ break;
+ prev_residual = curr_residual;
#endif
}