change the way that interations are performed at colorize algorithm

[colorize.git] / lib / multigrid.c

diff --git a/lib/multigrid.c b/lib/multigrid.c
index 4c0724b..f3097ee 100644 (file)
--- a/lib/multigrid.c
+++ b/lib/multigrid.c
@@ -6,8 +6,8 @@
 #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
@@ -1385,7 +1385,7 @@ void draw_flow(darray_t *dst, const darray_t *flow, const darray_t *flow_i)
  * 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;
@@ -1402,7 +1402,8 @@ int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **ma
                        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]);
@@ -1430,9 +1431,12 @@ int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **ma
 
                /* 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 (target_residual_change && (prev_residual-curr_residual)/prev_residual < target_residual_change)
                        break;
-               curr_residual = prev_residual;
+               prev_residual = curr_residual;
 #endif
        }