From 95bd4cc5a7922254b5ba2b227e58ae6c81c69c15 Mon Sep 17 00:00:00 2001
From: Sergey B Kirpichev <skirpichev@gmail.com>
Date: Sat, 24 Jul 2010 16:12:56 +0400
Subject: [PATCH] Add convergence checks for rk*imp

---
 ode-initval/rk2imp.c |   59 +++++++++++++++++++++++++++++++++---------
 ode-initval/rk4imp.c |   69 +++++++++++++++++++++++++++++++++++++++++++++++---
 2 files changed, 111 insertions(+), 17 deletions(-)

diff --git a/ode-initval/rk2imp.c b/ode-initval/rk2imp.c
index 4e94b08..4f282ab 100644
--- a/ode-initval/rk2imp.c
+++ b/ode-initval/rk2imp.c
@@ -40,6 +40,7 @@
 typedef struct
 {
   double *Y1;
+  double *Y1_tmp;
   double *y0;
   double *ytmp;
   double *y_onestep;
@@ -66,10 +67,20 @@ rk2imp_alloc (size_t dim)
       GSL_ERROR_NULL ("failed to allocate space for Y1", GSL_ENOMEM);
     }
 
+  state->Y1_tmp = (double *) malloc (dim * sizeof (double));
+
+  if (state->Y1_tmp == 0)
+    {
+      free (state->Y1);
+      free (state);
+      GSL_ERROR_NULL ("failed to allocate space for Y1_tmp", GSL_ENOMEM);
+    }
+
   state->ytmp = (double *) malloc (dim * sizeof (double));
 
   if (state->ytmp == 0)
     {
+      free (state->Y1_tmp);
       free (state->Y1);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
@@ -79,8 +90,9 @@ rk2imp_alloc (size_t dim)
 
   if (state->y0 == 0)
     {
-      free (state->Y1);
       free (state->ytmp);
+      free (state->Y1_tmp);
+      free (state->Y1);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
     }
@@ -89,9 +101,10 @@ rk2imp_alloc (size_t dim)
 
   if (state->y_onestep == 0)
     {
-      free (state->Y1);
-      free (state->ytmp);
       free (state->y0);
+      free (state->ytmp);
+      free (state->Y1_tmp);
+      free (state->Y1);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
     }
@@ -101,9 +114,10 @@ rk2imp_alloc (size_t dim)
   if (state->y0_orig == 0)
     {
       free (state->y_onestep);
-      free (state->Y1);
-      free (state->ytmp);
       free (state->y0);
+      free (state->ytmp);
+      free (state->Y1_tmp);
+      free (state->Y1);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for y0_orig", GSL_ENOMEM);
     }
@@ -127,21 +141,22 @@ rk2imp_step (double *y, rk2imp_state_t *state,
   */
 
   const double *y0 = state->y0;
-  double *Y1 = state->Y1;
+  double *Y1 = state->Y1, *Y1_tmp = state->Y1_tmp;
   double *ytmp = state->ytmp;
-  int max_iter=3;
-  int nu;
+  const size_t maxiter = 100;
+  size_t iter = 1;
+  const double normtol = GSL_DBL_EPSILON;
+  double norm;
   size_t i;
 
   /* iterative solution of Y1 = y0 + h/2 * f(t + h/2, Y1) 
      Y1 should include initial values at call.
-
-     Note: This method does not check for convergence of the
-     iterative solution! 
   */
 
-  for (nu = 0; nu < max_iter; nu++)
+  do
     {
+      DBL_MEMCPY (Y1_tmp, Y1, dim);
+
       for (i = 0; i < dim; i++)
         {
           ytmp[i] = y0[i] + 0.5 * h * Y1[i];
@@ -155,8 +170,24 @@ rk2imp_step (double *y, rk2imp_state_t *state,
 	    return s;
 	  }    
       }
+
+      norm = 0.0; /* reset norm */
+
+      for (i = 0; i < dim; i++)
+        {
+          norm = GSL_MAX (norm, fabs (Y1[i] - Y1_tmp[i]));
+        }
+
+      iter++;
+
+      if (iter > maxiter)
+        {
+          GSL_ERROR ("reached iteration limit for "
+                     "fixed-point equation solver", GSL_EMAXITER);
+        }
     }
-  
+  while (norm > normtol);
+
   /* assignment */
 
   for (i = 0; i < dim; i++)
@@ -294,6 +325,7 @@ rk2imp_reset (void *vstate, size_t dim)
   rk2imp_state_t *state = (rk2imp_state_t *) vstate;
 
   DBL_ZERO_MEMSET (state->Y1, dim);
+  DBL_ZERO_MEMSET (state->Y1_tmp, dim);
   DBL_ZERO_MEMSET (state->ytmp, dim);
   DBL_ZERO_MEMSET (state->y0, dim);
   DBL_ZERO_MEMSET (state->y_onestep, dim);
@@ -316,6 +348,7 @@ rk2imp_free (void *vstate)
   rk2imp_state_t *state = (rk2imp_state_t *) vstate;
 
   free (state->Y1);
+  free (state->Y1_tmp);
   free (state->ytmp);
   free (state->y0);
   free (state->y_onestep);
diff --git a/ode-initval/rk4imp.c b/ode-initval/rk4imp.c
index e569b7e..f7e0a1b 100644
--- a/ode-initval/rk4imp.c
+++ b/ode-initval/rk4imp.c
@@ -40,7 +40,9 @@
 typedef struct
 {
   double *k1nu;
+  double *k1nu_tmp;
   double *k2nu;
+  double *k2nu_tmp;
   double *ytmp1;
   double *ytmp2;
   double *y0;
@@ -68,20 +70,43 @@ rk4imp_alloc (size_t dim)
       GSL_ERROR_NULL ("failed to allocate space for k1nu", GSL_ENOMEM);
     }
 
+  state->k1nu_tmp = (double *) malloc (dim * sizeof (double));
+
+  if (state->k1nu_tmp == 0)
+    {
+      free (state->k1nu);
+      free (state);
+      GSL_ERROR_NULL ("failed to allocate space for k1nu_tmp", GSL_ENOMEM);
+    }
+
   state->k2nu = (double *) malloc (dim * sizeof (double));
 
   if (state->k2nu == 0)
     {
+      free (state->k1nu_tmp);
       free (state->k1nu);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for k2nu", GSL_ENOMEM);
     }
 
+  state->k2nu_tmp = (double *) malloc (dim * sizeof (double));
+
+  if (state->k2nu_tmp == 0)
+    {
+      free (state->k2nu);
+      free (state->k1nu_tmp);
+      free (state->k1nu);
+      free (state);
+      GSL_ERROR_NULL ("failed to allocate space for k2nu_tmp", GSL_ENOMEM);
+    }
+
   state->ytmp1 = (double *) malloc (dim * sizeof (double));
 
   if (state->ytmp1 == 0)
     {
+      free (state->k2nu_tmp);
       free (state->k2nu);
+      free (state->k1nu_tmp);
       free (state->k1nu);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for ytmp1", GSL_ENOMEM);
@@ -92,7 +117,9 @@ rk4imp_alloc (size_t dim)
   if (state->ytmp2 == 0)
     {
       free (state->ytmp1);
+      free (state->k2nu_tmp);
       free (state->k2nu);
+      free (state->k1nu_tmp);
       free (state->k1nu);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for ytmp2", GSL_ENOMEM);
@@ -104,7 +131,9 @@ rk4imp_alloc (size_t dim)
     {
       free (state->ytmp2);
       free (state->ytmp1);
+      free (state->k2nu_tmp);
       free (state->k2nu);
+      free (state->k1nu_tmp);
       free (state->k1nu);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
@@ -117,7 +146,9 @@ rk4imp_alloc (size_t dim)
       free (state->y0);
       free (state->ytmp2);
       free (state->ytmp1);
+      free (state->k2nu_tmp);
       free (state->k2nu);
+      free (state->k1nu_tmp);
       free (state->k1nu);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for y0_orig", GSL_ENOMEM);
@@ -131,7 +162,9 @@ rk4imp_alloc (size_t dim)
       free (state->y0);
       free (state->ytmp2);
       free (state->ytmp1);
+      free (state->k2nu_tmp);
       free (state->k2nu);
+      free (state->k1nu_tmp);
       free (state->k1nu);
       free (state);
       GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
@@ -161,12 +194,16 @@ rk4imp_step (double *y, rk4imp_state_t *state,
   */
 
   const double ir3 = 1.0 / M_SQRT3;
-  const int iter_steps = 3;
-  int nu;
+  const size_t maxiter = 100;
+  size_t iter = 1;
+  const double normtol = GSL_DBL_EPSILON;
+  double norm;
   size_t i;
 
   double *const k1nu = state->k1nu;
+  double *const k1nu_tmp = state->k1nu_tmp;
   double *const k2nu = state->k2nu;
+  double *const k2nu_tmp = state->k2nu_tmp;
   double *const ytmp1 = state->ytmp1;
   double *const ytmp2 = state->ytmp2;
 
@@ -176,8 +213,11 @@ rk4imp_step (double *y, rk4imp_state_t *state,
      iterative solution! 
   */
 
-  for (nu = 0; nu < iter_steps; nu++)
+  do
     {
+      DBL_MEMCPY (k1nu_tmp, k1nu, dim);
+      DBL_MEMCPY (k2nu_tmp, k2nu, dim);
+
       for (i = 0; i < dim; i++)
         {
           ytmp1[i] =
@@ -192,7 +232,7 @@ rk4imp_step (double *y, rk4imp_state_t *state,
 	if (s != GSL_SUCCESS)
 	  {
 	    return s;
-	  }    
+	  }
       }
       {
         int s =
@@ -203,7 +243,24 @@ rk4imp_step (double *y, rk4imp_state_t *state,
 	    return s;
 	  }    
       }
+
+      norm = 0.0; /* reset norm */
+
+      for (i = 0; i < dim; i++)
+        {
+          norm = GSL_MAX (norm, fabs (k1nu[i] - k1nu_tmp[i]));
+          norm = GSL_MAX (norm, fabs (k2nu[i] - k2nu_tmp[i]));
+        }
+
+      iter++;
+
+      if (iter > maxiter)
+        {
+          GSL_ERROR ("reached iteration limit for "
+                     "fixed-point equation solver", GSL_EMAXITER);
+        }
     }
+  while (norm > normtol);
 
   /* assignment */
   
@@ -347,7 +404,9 @@ rk4imp_reset (void *vstate, size_t dim)
   DBL_ZERO_MEMSET (state->y0_orig, dim);
   DBL_ZERO_MEMSET (state->y0, dim);
   DBL_ZERO_MEMSET (state->k1nu, dim);
+  DBL_ZERO_MEMSET (state->k1nu_tmp, dim);
   DBL_ZERO_MEMSET (state->k2nu, dim);
+  DBL_ZERO_MEMSET (state->k2nu_tmp, dim);
   DBL_ZERO_MEMSET (state->ytmp1, dim);
   DBL_ZERO_MEMSET (state->ytmp2, dim);
 
@@ -371,7 +430,9 @@ rk4imp_free (void *vstate)
   free (state->y0_orig);
   free (state->y0);
   free (state->k1nu);
+  free (state->k1nu_tmp);
   free (state->k2nu);
+  free (state->k2nu_tmp);
   free (state->ytmp1);
   free (state->ytmp2);
   free (state);
-- 
1.5.6.5

