Solve the discretized system exactly on the coarsest level.

1 files changed, 23 insertions, 12 deletions

diff --git a/mg2d.c b/mg2d.c
index 9216af4..902c3cc 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -297,12 +297,9 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
 
     /* handle coarsest grid */
     if (!level->child) {
-
-        for (int j = 0; j < 16; j++) {
-            ret = mg_relax_step(ctx, level, "coarse-step");
-            if (ret < 0)
-                return ret;
-        }
+        ret = mg_relax_step(ctx, level, "coarse-step");
+        if (ret < 0)
+            return ret;
         level->count_cycles++;
 
         goto finish;
@@ -553,7 +550,7 @@ static void mg_level_free(MG2DLevel **plevel)
     *plevel = NULL;
 }
 
-static MG2DLevel *mg_level_alloc(const size_t domain_size)
+static MG2DLevel *mg_level_alloc(enum EGSType type, const size_t domain_size)
 {
     MG2DLevel *level;
     int ret;
@@ -567,7 +564,7 @@ static MG2DLevel *mg_level_alloc(const size_t domain_size)
         goto fail;
     level->prolong_tmp_stride = domain_size + 1;
 
-    level->solver = mg2di_egs_alloc(EGS_SOLVER_RELAXATION,
+    level->solver = mg2di_egs_alloc(type,
                                     (size_t [2]){domain_size, domain_size});
     if (!level->solver)
         goto fail;
@@ -591,9 +588,11 @@ static int log2i(int n)
 static int mg_levels_alloc(MG2DContext *ctx, size_t domain_size)
 {
     MG2DInternal *priv = ctx->priv;
-    MG2DLevel *cur;
+    MG2DLevel *cur, **last;
+    size_t last_size;
+    unsigned int last_depth;
 
-    priv->root = mg_level_alloc(domain_size);
+    priv->root = mg_level_alloc(EGS_SOLVER_RELAXATION, domain_size);
     if (!priv->root)
         return -ENOMEM;
 
@@ -601,14 +600,14 @@ static int mg_levels_alloc(MG2DContext *ctx, size_t domain_size)
     // the children all have sizes 2**n + 1
     // but we skip the closest lower one if it is too close
     if (domain_size <= 2)
-        return 0;
+        goto finish;
     domain_size = (1 << log2i(domain_size - 2)) + 1;
     if ((double)priv->root->solver->domain_size[0] / domain_size  < 1.5)
         domain_size = (domain_size >> 1) + 1;
 
     cur = priv->root;
     for (int i = 0; i < LEVELS_MAX && domain_size > 4; i++) {
-        cur->child = mg_level_alloc(domain_size);
+        cur->child = mg_level_alloc(EGS_SOLVER_RELAXATION, domain_size);
         if (!cur->child)
             return -ENOMEM;
         cur->child->depth = i + 1;
@@ -617,6 +616,18 @@ static int mg_levels_alloc(MG2DContext *ctx, size_t domain_size)
         domain_size = (domain_size >> 1) + 1;
     }
 
+finish:
+    last = &priv->root;
+    while ((*last)->child)
+        last = &((*last)->child);
+    last_size = (*last)->solver->domain_size[0];
+    last_depth = (*last)->depth;
+    mg_level_free(last);
+    *last = mg_level_alloc(EGS_SOLVER_EXACT, last_size);
+    if (!*last)
+        return -ENOMEM;
+    (*last)->depth = last_depth;
+
     return 0;
 }