text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk.c
 * \brief Source file with common variables and functions to optimize
 *   Runge-Kutta methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#if HAVE_MPI
#include <mpi.h>
#endif
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_2_2.h"
#include "rk_3_2.h"
#include "rk_3_3.h"
#include "rk_4_2.h"
#include "rk_4_3.h"
#include "rk_4_4.h"
#include "rk_5_2.h"
#include "rk_5_3.h"
#include "rk_5_4.h"
#include "rk_6_2.h"
#include "rk_6_3.h"
#include "rk_6_4.h"

#define DEBUG_RK 0              ///< macro to debug.

/**
 * Function to print the t-b Runge-Kutta coefficients.
 */
void
rk_print_tb (Optimize * tb,     ///< Optimize struct.
             char *label,       ///< label.
             FILE * file)       ///< file.
{
  long double *x;
  unsigned int i, j, k;
  x = tb->coefficient;
  fprintf (file, "%s: t1=%.19Le\n", label, x[0]);
  for (i = 2, k = 0; i <= tb->nsteps; ++i)
    {
      fprintf (file, "%s: t%u=%.19Le\n", label, i, x[++k]);
      for (j = 0; j < i; ++j)
        fprintf (file, "%s: b%u%u=%.19Le\n", label, i, j, x[++k]);
    }
}

/**
 * Function to print the e Runge-Kutta coefficients.
 */
void
rk_print_e (Optimize * tb,      ///< Optimize struct.
            char *label,        ///< label.
            FILE * file)        ///< file.
{
  long double *x;
  unsigned int i, k, nsteps;
  x = tb->coefficient;
  nsteps = tb->nsteps;
  k = (nsteps + 2) * (nsteps + 1) / 2 - 2;
  for (i = 0; i < nsteps - 1; ++i)
    fprintf (file, "%s: e%u%u=%.19Le\n", label, nsteps, i, x[k++]);
}

/**
 * Function to print in a maxima file the Runge-Kutta coefficients.
 */
static void
rk_print (RK * rk,              ///< RK struct.
          FILE * file)          ///< file.
{
  Optimize *tb, *ac;
  long double *x, *y;
  unsigned int i, j, k, l, nsteps;
  tb = rk->tb;
  ac = rk->ac;
  x = tb->coefficient;
  y = ac->coefficient;
  fprintf (file, "t1:%.19Le;\n", x[0]);
  nsteps = tb->nsteps;
  for (i = 2, k = l = 0; i <= nsteps; ++i)
    {
      fprintf (file, "t%u:%.19Le;\n", i, x[++k]);
      for (j = 0; j < i; ++j)
        fprintf (file, "b%u%u:%.19Le;\n", i, j, x[++k]);
      if (!rk->strong)
        continue;
      for (j = 0; j < i; ++j)
        fprintf (file, "a%u%u:%.19Le;\n", i, j, y[l++]);
      for (j = 0; j < i; ++j)
        fprintf (file, "c%u%u:%.19Le;\n", i, j, y[l++]);
    }
  if (rk->pair)
    for (i = 0; i < nsteps - 1; ++i)
      fprintf (file, "e%u%u:%.19Le;\n", nsteps, i, x[++k]);
}

/**
 * Function to print a maxima format file to check the accuracy order of the
 * Runge-Kutta simple stable methods.
 */
static void
rk_print_maxima (FILE * file,   ///< file.
                 unsigned int nsteps,   ///< steps number.
                 unsigned int ncoefficients,    ///< coefficients number.
                 unsigned int order,    ///< accuracy order.
                 char label)    ///< coefficient label.
{
  unsigned int i, j, k, l;
  // b_{ij}=1 (1st order)
  for (i = 0; i < ncoefficients; ++i)
    fprintf (file, "%c%u%u+", label, nsteps, i);
  fprintf (file, "-1;\n");
  // b_{ij}t_j=1/2 (2nd order)
  for (i = 1; i < ncoefficients; ++i)
    fprintf (file, "%c%u%u*t%u+", label, nsteps, i, i);
  fprintf (file, "-1/2;\n");
  if (order < 2)
    return;
  // b_{ij}t_j^2=1/3 (3rd order)
  for (i = 1; i < ncoefficients; ++i)
    fprintf (file, "%c%u%u*t%u^2+", label, nsteps, i, i);
  fprintf (file, "-1/3;\n");
  if (order < 3)
    return;
  // b_{ij}b_{jk}t_k=1/6 (3rd order)
  for (i = 2; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*(", label, nsteps, i);
      for (j = 1; j < i; ++j)
        fprintf (file, "b%u%u*t%u+", i, j, j);
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/6;\n");
  // b_{ij}t_j^3=1/4 (4th order)
  for (i = 1; i < ncoefficients; ++i)
    fprintf (file, "%c%u%u*t%u^3+", label, nsteps, i, i);
  fprintf (file, "-1/4;\n");
  if (order < 4)
    return;
  // b_{ij}b_{jk}b_{kl}t_l=1/24 (4th order)
  for (i = 3; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*(", label, nsteps, i);
      for (j = 2; j < i; ++j)
        {
          fprintf (file, "b%u%u*(", i, j);
          for (k = 1; k < j; ++k)
            fprintf (file, "b%u%u*t%u+", j, k, k);
          fprintf (file, "0)+");
        }
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/24;\n");
  // b_{ij}b_{jk}t_k^2=1/12 (4th order)
  for (i = 2; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*(", label, nsteps, i);
      for (j = 1; j < i; ++j)
        fprintf (file, "b%u%u*t%u^2+", i, j, j);
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/12;\n");
  // b_{ij}t_jb_{jk}t_k=1/8 (4th order)
  for (i = 2; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*t%u*(", label, nsteps, i, i);
      for (j = 1; j < i; ++j)
        fprintf (file, "b%u%u*t%u+", i, j, j);
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/8;\n");
  // b_{ij}t_j^4=1/5 (5th order)
  for (i = 1; i < ncoefficients; ++i)
    fprintf (file, "%c%u%u*t%u^4+", label, nsteps, i, i);
  fprintf (file, "-1/5;\n");
  if (order < 5)
    return;
  // b_{ij}t_j^2b_{jk}t_k^2=1/10 (5th order)
  for (i = 2; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*t%u^2*(", label, nsteps, i, i);
      for (j = 1; j < i; ++j)
        fprintf (file, "b%u%u*t%u+", i, j, j);
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/10;\n");
  // b_{ij}b_{jk}t_k^3=1/20 (5th order)
  for (i = 2; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*(", label, nsteps, i);
      for (j = 1; j < i; ++j)
        fprintf (file, "b%u%u*t%u^3+", i, j, j);
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/20;\n");
  // b_{ij}(b_{jk}t_k)^2=1/20 (5th order)
  for (i = 2; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*(", label, nsteps, i);
      for (j = 1; j < i; ++j)
        fprintf (file, "b%u%u*t%u+", i, j, j);
      fprintf (file, "0)^2+");
    }
  fprintf (file, "-1/20;\n");
  // b_{ij}t_jb_{jk}t_k^2=1/15 (5th order)
  for (i = 2; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*t%u*(", label, nsteps, i, i);
      for (j = 1; j < i; ++j)
        fprintf (file, "b%u%u*t%u^2+", i, j, j);
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/8;\n");
  // b_{ij}t_jb_{jk}b_{kl}t_l=1/24 (5th order)
  for (i = 3; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*t%u*(", label, nsteps, i, i);
      for (j = 2; j < i; ++j)
        {
          fprintf (file, "b%u%u*(", i, j);
          for (k = 1; k < j; ++k)
            fprintf (file, "b%u%u*t%u+", j, k, k);
          fprintf (file, "0)+");
        }
      fprintf (file, "0)+");
    }
  fprintf (file, "-7/120;\n");
  // b_{ij}b_{jk}b_{kl}t_l^2=1/60 (5th order)
  for (i = 3; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*(", label, nsteps, i);
      for (j = 2; j < i; ++j)
        {
          fprintf (file, "b%u%u*(", i, j);
          for (k = 1; k < j; ++k)
            fprintf (file, "b%u%u*t%u^2+", j, k, k);
          fprintf (file, "0)+");
        }
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/60;\n");
  // b_{ij}b_{jk}b_{kl}b_{lm}t_m=1/120 (5th order)
  for (i = 4; i < ncoefficients; ++i)
    {
      fprintf (file, "%c%u%u*(", label, nsteps, i);
      for (j = 3; j < i; ++j)
        {
          fprintf (file, "b%u%u*(", i, j);
          for (k = 2; k < j; ++k)
            {
              fprintf (file, "b%u%u*(", j, k);
              for (l = 1; l < k; ++l)
                fprintf (file, "b%u%u*t%u+", k, l, l);
              fprintf (file, "0)+");
            }
          fprintf (file, "0)+");
        }
      fprintf (file, "0)+");
    }
  fprintf (file, "-1/120;\n");
  // b_{ij}t_j^5=1/6 (6th order)
  for (i = 1; i < ncoefficients; ++i)
    fprintf (file, "%c%u%u*t%u^5+", label, nsteps, i, i);
  fprintf (file, "-1/6;\n");
}

/**
 * Function to print in a maxima file the a-c Runge-Kutta equations.
 */
static void
ac_print_maxima (FILE * file,   ///< file.
                 unsigned int nsteps)   ///< steps number.
{
  unsigned int i, j, k;
  for (i = 2; i <= nsteps; ++i)
    {
      // a_{i,j}=1
      for (j = 0; j < i; ++j)
        fprintf (file, "a%u%u+", i, j);
      fprintf (file, "-1;\n");
      // a_{i0}c_{i0}+a_{ij}b_{j0}=b_{i0}
      fprintf (file, "a%u0*c%u0+a%u1*t1+", i, i, i);
      for (j = 2; j < i; ++j)
        fprintf (file, "a%u%u*b%u0+", i, j, j);
      fprintf (file, "-b%u0;\n", i);
      // a_{ij}c_{ij}+a_{ik}b_{kj}=b_{ij}
      for (j = 1; j < i; ++j)
        {
          fprintf (file, "a%u%u*c%u%u+", i, j, i, j);
          for (k = j; ++k < i;)
            fprintf (file, "a%u%u*b%u%u+", i, k, k, j);
          fprintf (file, "-b%u%u;\n", i, j);
        }
    }
}

/**
 * Function to get \f$a_{ij}\f$ and \f$c_{ij}\f$ coefficients of the Runge-Kutta
 * 2nd step.
 */
static int
rk_ac_2 (RK * rk)               ///< RK struct.
{
  long double *tb, *ac, *r;
  register long double ac0;
#if DEBUG_RK
  fprintf (stderr, "rk_ac_2: start\n");
#endif
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  r = rk->ac->random_data;
  c21 (ac) = r[0];
  a21 (ac) = b21 (tb) / c21 (ac);
  a20 (ac) = 1.L - a21 (ac);
  ac0 = b20 (tb) - a21 (ac) * t1 (tb);
  c20 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / a20 (ac);
#if DEBUG_RK
  fprintf (stderr, "rk_ac_2: a20=%Lg c20=%Lg\n", a20 (ac), c20 (ac));
  fprintf (stderr, "rk_ac_2: a21=%Lg c21=%Lg\n", a21 (ac), c21 (ac));
  fprintf (stderr, "rk_ac_2: end\n");
#endif
  if (isnan (c20 (ac)) || isnan (a21 (ac)))
    return 0;
  return 1;
}

/**
 * Function to get \f$a_{ij}\f$ and \f$c_{ij}\f$ coefficients of the Runge-Kutta
 * 3rd step.
 */
static int
rk_ac_3 (RK * rk)               ///< RK struct.
{
  long double *tb, *ac, *r;
  register long double ac0;
#if DEBUG_RK
  fprintf (stderr, "rk_ac_3: start\n");
#endif
  if (!rk_ac_2 (rk))
    {
#if DEBUG_RK
      fprintf (stderr, "rk_ac_3: end\n");
#endif
      return 0;
    }
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  r = rk->ac->random_data;
  c31 (ac) = r[1];
  c32 (ac) = r[2];
  a32 (ac) = b32 (tb) / c32 (ac);
  ac0 = b31 (tb) - a32 (ac) * b21 (tb);
  a31 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c31 (ac);
  a30 (ac) = 1.L - a31 (ac) - a32 (ac);
  ac0 = b30 (tb) - a31 (ac) * t1 (tb) - a32 (ac) * b20 (tb);
  c30 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / a30 (ac);
#if DEBUG_RK
  fprintf (stderr, "rk_ac_3: a30=%Lg c30=%Lg\n", a30 (ac), c30 (ac));
  fprintf (stderr, "rk_ac_3: a31=%Lg c31=%Lg\n", a31 (ac), c31 (ac));
  fprintf (stderr, "rk_ac_3: a32=%Lg c32=%Lg\n", a32 (ac), c32 (ac));
  fprintf (stderr, "rk_ac_3: end\n");
#endif
  if (isnan (c30 (ac)) || isnan (a31 (ac)) || isnan (a32 (ac)))
    return 0;
  return 1;
}

/**
 * Function to get \f$a_{ij}\f$ and \f$c_{ij}\f$ coefficients of the Runge-Kutta
 * 4th step.
 */
static int
rk_ac_4 (RK * rk)               ///< RK struct.
{
  long double *tb, *ac, *r;
  register long double ac0;
#if DEBUG_RK
  fprintf (stderr, "rk_ac_4: start\n");
#endif
  if (!rk_ac_3 (rk))
    {
#if DEBUG_RK
      fprintf (stderr, "rk_ac_4: end\n");
#endif
      return 0;
    }
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  r = rk->ac->random_data;
  c41 (ac) = r[3];
  c42 (ac) = r[4];
  c43 (ac) = r[5];
  a43 (ac) = b43 (tb) / c43 (ac);
  ac0 = b42 (tb) - a43 (ac) * b32 (tb);
  a42 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c42 (ac);
  ac0 = b41 (tb) - a42 (ac) * b21 (tb) - a43 (ac) * b31 (tb);
  a41 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c41 (ac);
  a40 (ac) = 1.L - a41 (ac) - a42 (ac) - a43 (ac);
  ac0 = b40 (tb) - a41 (ac) * t1 (tb) - a42 (ac) * b20 (tb)
    - a43 (ac) * b30 (tb);
  c40 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / a40 (ac);
#if DEBUG_RK
  fprintf (stderr, "rk_ac_4: a40=%Lg c40=%Lg\n", a40 (ac), c40 (ac));
  fprintf (stderr, "rk_ac_4: a41=%Lg c41=%Lg\n", a41 (ac), c41 (ac));
  fprintf (stderr, "rk_ac_4: a42=%Lg c42=%Lg\n", a42 (ac), c42 (ac));
  fprintf (stderr, "rk_ac_4: a43=%Lg c43=%Lg\n", a43 (ac), c43 (ac));
  fprintf (stderr, "rk_ac_4: end\n");
#endif
  if (isnan (c40 (ac)) || isnan (a41 (ac)) || isnan (a42 (ac))
      || isnan (a43 (ac)))
    return 0;
  return 1;
}

/**
 * Function to get \f$a_{ij}\f$ and \f$c_{ij}\f$ coefficients of the Runge-Kutta
 * 5th step.
 */
static int
rk_ac_5 (RK * rk)               ///< RK struct.
{
  long double *tb, *ac, *r;
  register long double ac0;
#if DEBUG_RK
  fprintf (stderr, "rk_ac_5: start\n");
#endif
  if (!rk_ac_4 (rk))
    {
#if DEBUG_RK
      fprintf (stderr, "rk_ac_5: end\n");
#endif
      return 0;
    }
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  r = rk->ac->random_data;
  c51 (ac) = r[6];
  c52 (ac) = r[7];
  c53 (ac) = r[8];
  c54 (ac) = r[9];
  a54 (ac) = b54 (tb) / c54 (ac);
  ac0 = b53 (tb) - a54 (ac) * b43 (tb);
  a53 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c53 (ac);
  ac0 = b52 (tb) - a53 (ac) * b32 (tb) - a54 (ac) * b42 (tb);
  a52 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c52 (ac);
  ac0 = b51 (tb) - a52 (ac) * b21 (tb) - a53 (ac) * b31 (tb)
    - a54 (ac) * b41 (tb);
  a51 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c51 (ac);
  a50 (ac) = 1.L - a51 (ac) - a52 (ac) - a53 (ac) - a54 (ac);
  ac0 = b50 (tb) - a51 (ac) * t1 (tb) - a52 (ac) * b20 (tb)
    - a53 (ac) * b30 (tb) - a54 (ac) * b40 (tb);
  c50 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / a50 (ac);
#if DEBUG_RK
  fprintf (stderr, "rk_ac_5: a50=%Lg c50=%Lg\n", a50 (ac), c50 (ac));
  fprintf (stderr, "rk_ac_5: a51=%Lg c51=%Lg\n", a51 (ac), c51 (ac));
  fprintf (stderr, "rk_ac_5: a52=%Lg c52=%Lg\n", a52 (ac), c52 (ac));
  fprintf (stderr, "rk_ac_5: a53=%Lg c53=%Lg\n", a53 (ac), c53 (ac));
  fprintf (stderr, "rk_ac_5: a54=%Lg c54=%Lg\n", a54 (ac), c54 (ac));
  fprintf (stderr, "rk_ac_5: end\n");
#endif
  if (isnan (c50 (ac)) || isnan (a51 (ac)) || isnan (a52 (ac))
      || isnan (a53 (ac)) || isnan (a54 (ac)))
    return 0;
  return 1;
}

/**
 * Function to get \f$a_{ij}\f$ and \f$c_{ij}\f$ coefficients of the Runge-Kutta
 * 6th step.
 */
static int
rk_ac_6 (RK * rk)               ///< RK struct.
{
  long double *tb, *ac, *r;
  register long double ac0;
#if DEBUG_RK
  fprintf (stderr, "rk_ac_6: start\n");
#endif
  if (!rk_ac_5 (rk))
    {
#if DEBUG_RK
      fprintf (stderr, "rk_ac_6: end\n");
#endif
      return 0;
    }
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  r = rk->ac->random_data;
  c61 (ac) = r[10];
  c62 (ac) = r[11];
  c63 (ac) = r[12];
  c64 (ac) = r[13];
  c65 (ac) = r[14];
  a65 (ac) = b65 (tb) / c65 (ac);
  ac0 = b64 (tb) - a65 (ac) * b54 (tb);
  a64 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c64 (ac);
  ac0 = b63 (tb) - a64 (ac) * b43 (tb) - a65 (ac) * b53 (tb);
  a63 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c63 (ac);
  ac0 = b62 (tb) - a63 (ac) * b32 (tb) - a64 (ac) * b42 (tb)
    - a65 (ac) * b52 (tb);
  a62 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c62 (ac);
  ac0 = b61 (tb) - a62 (ac) * b21 (tb) - a63 (ac) * b31 (tb)
    - a64 (ac) * b41 (tb) - a65 (ac) * b51 (tb);
  a61 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / c61 (ac);
  a60 (ac) = 1.L - a61 (ac) - a62 (ac) - a63 (ac) - a64 (ac) - a65 (ac);
  ac0 = b60 (tb) - a61 (ac) * t1 (tb) - a62 (ac) * b20 (tb)
    - a63 (ac) * b30 (tb) - a64 (ac) * b40 (tb) - a65 (ac) * b50 (tb);
  c60 (ac) = fabsl (ac0) < LDBL_EPSILON ? 0.L : ac0 / a60 (ac);
#if DEBUG_RK
  fprintf (stderr, "rk_ac_6: a60=%Lg c60=%Lg\n", a60 (ac), c60 (ac));
  fprintf (stderr, "rk_ac_6: a61=%Lg c61=%Lg\n", a61 (ac), c61 (ac));
  fprintf (stderr, "rk_ac_6: a62=%Lg c62=%Lg\n", a62 (ac), c62 (ac));
  fprintf (stderr, "rk_ac_6: a63=%Lg c63=%Lg\n", a63 (ac), c63 (ac));
  fprintf (stderr, "rk_ac_6: a64=%Lg c64=%Lg\n", a64 (ac), c64 (ac));
  fprintf (stderr, "rk_ac_6: a65=%Lg c65=%Lg\n", a65 (ac), c65 (ac));
  fprintf (stderr, "rk_ac_6: end\n");
#endif
  if (isnan (c60 (ac)) || isnan (a61 (ac)) || isnan (a62 (ac))
      || isnan (a63 (ac)) || isnan (a64 (ac)) || isnan (a65 (ac)))
    return 0;
  return 1;
}

/**
 * Function to get the objective function of 2 steps Runge-Kutta methods.
 *
 * \return objective function value.
 */
static long double
rk_objective_ac_2 (RK * rk)     ///< RK struct.
{
  long double *tb, *ac;
  long double k;
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_2: start\n");
#endif
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  k = fminl (0.L, a20 (ac));
  if (a21 (ac) < 0.L)
    k += a21 (ac);
  if (k < 0.L)
    {
      k = 20.L - k;
      goto end;
    }
  k = fminl (0.L, c20 (ac));
  if (c21 (ac) < 0.L)
    k += c21 (ac);
  if (k < 0.L)
    {
      k = 10.L - k;
      goto end;
    }
  k = 1.L / rk_cfl_2 (tb, ac);
end:
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_2: objective=%Lg\n", k);
  fprintf (stderr, "rk_objective ac_2: end\n");
#endif
  return k;
}

/**
 * Function to get the objective function of 3 steps Runge-Kutta methods.
 *
 * \return objective function value.
 */
static long double
rk_objective_ac_3 (RK * rk)     ///< RK struct.
{
  long double *tb, *ac;
  long double k;
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_3: start\n");
#endif
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  k = fminl (0.L, a20 (ac));
  if (a21 (ac) < 0.L)
    k += a21 (ac);
  if (a30 (ac) < 0.L)
    k += a30 (ac);
  if (a31 (ac) < 0.L)
    k += a31 (ac);
  if (a32 (ac) < 0.L)
    k += a32 (ac);
  if (k < 0.L)
    {
      k = 20.L - k;
      goto end;
    }
  k = fminl (0.L, c20 (ac));
  if (c21 (ac) < 0.L)
    k += c21 (ac);
  if (c30 (ac) < 0.L)
    k += c30 (ac);
  if (c31 (ac) < 0.L)
    k += c31 (ac);
  if (c32 (ac) < 0.L)
    k += c32 (ac);
  if (k < 0.L)
    {
      k = 10.L - k;
      goto end;
    }
  k = 1.L / rk_cfl_3 (tb, ac);
end:
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_3: objective=%Lg\n", k);
  fprintf (stderr, "rk_objective ac_3: end\n");
#endif
  return k;
}

/**
 * Function to get the objective function of 4 steps Runge-Kutta methods.
 *
 * \return objective function value.
 */
static long double
rk_objective_ac_4 (RK * rk)     ///< RK struct.
{
  long double *tb, *ac;
  long double k;
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_4: start\n");
#endif
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  k = fminl (0.L, a20 (ac));
  if (a21 (ac) < 0.L)
    k += a21 (ac);
  if (a30 (ac) < 0.L)
    k += a30 (ac);
  if (a31 (ac) < 0.L)
    k += a31 (ac);
  if (a32 (ac) < 0.L)
    k += a32 (ac);
  if (a40 (ac) < 0.L)
    k += a40 (ac);
  if (a41 (ac) < 0.L)
    k += a41 (ac);
  if (a42 (ac) < 0.L)
    k += a42 (ac);
  if (a43 (ac) < 0.L)
    k += a43 (ac);
  if (k < 0.L)
    {
      k = 20.L - k;
      goto end;
    }
  k = fminl (0.L, c20 (ac));
  if (c21 (ac) < 0.L)
    k += c21 (ac);
  if (c30 (ac) < 0.L)
    k += c30 (ac);
  if (c31 (ac) < 0.L)
    k += c31 (ac);
  if (c32 (ac) < 0.L)
    k += c32 (ac);
  if (c40 (ac) < 0.L)
    k += c40 (ac);
  if (c41 (ac) < 0.L)
    k += c41 (ac);
  if (c42 (ac) < 0.L)
    k += c42 (ac);
  if (c43 (ac) < 0.L)
    k += c43 (ac);
  if (k < 0.L)
    {
      k = 10.L - k;
      goto end;
    }
  k = 1.L / rk_cfl_4 (tb, ac);
end:
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_4: objective=%Lg\n", k);
  fprintf (stderr, "rk_objective ac_4: end\n");
#endif
  return k;
}

/**
 * Function to get the objective function of 5 steps Runge-Kutta methods.
 *
 * \return objective function value.
 */
static long double
rk_objective_ac_5 (RK * rk)     ///< RK struct.
{
  long double *tb, *ac;
  long double k;
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_5: start\n");
#endif
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  k = fminl (0.L, a20 (ac));
  if (a21 (ac) < 0.L)
    k += a21 (ac);
  if (a30 (ac) < 0.L)
    k += a30 (ac);
  if (a31 (ac) < 0.L)
    k += a31 (ac);
  if (a32 (ac) < 0.L)
    k += a32 (ac);
  if (a40 (ac) < 0.L)
    k += a40 (ac);
  if (a41 (ac) < 0.L)
    k += a41 (ac);
  if (a42 (ac) < 0.L)
    k += a42 (ac);
  if (a43 (ac) < 0.L)
    k += a43 (ac);
  if (a50 (ac) < 0.L)
    k += a50 (ac);
  if (a51 (ac) < 0.L)
    k += a51 (ac);
  if (a52 (ac) < 0.L)
    k += a52 (ac);
  if (a53 (ac) < 0.L)
    k += a53 (ac);
  if (a54 (ac) < 0.L)
    k += a54 (ac);
  if (k < 0.L)
    {
      k = 20.L - k;
      goto end;
    }
  k = fminl (0.L, c20 (ac));
  if (c21 (ac) < 0.L)
    k += c21 (ac);
  if (c30 (ac) < 0.L)
    k += c30 (ac);
  if (c31 (ac) < 0.L)
    k += c31 (ac);
  if (c32 (ac) < 0.L)
    k += c32 (ac);
  if (c40 (ac) < 0.L)
    k += c40 (ac);
  if (c41 (ac) < 0.L)
    k += c41 (ac);
  if (c42 (ac) < 0.L)
    k += c42 (ac);
  if (c43 (ac) < 0.L)
    k += c43 (ac);
  if (c50 (ac) < 0.L)
    k += c50 (ac);
  if (c51 (ac) < 0.L)
    k += c51 (ac);
  if (c52 (ac) < 0.L)
    k += c52 (ac);
  if (c53 (ac) < 0.L)
    k += c53 (ac);
  if (c54 (ac) < 0.L)
    k += c54 (ac);
  if (k < 0.L)
    {
      k = 10.L - k;
      goto end;
    }
  k = 1.L / rk_cfl_5 (tb, ac);
end:
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_5: objective=%Lg\n", k);
  fprintf (stderr, "rk_objective ac_5: end\n");
#endif
  return k;
}

/**
 * Function to get the objective function of 6 steps Runge-Kutta methods.
 *
 * \return objective function value.
 */
static long double
rk_objective_ac_6 (RK * rk)     ///< RK struct.
{
  long double *tb, *ac;
  long double k;
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_6: start\n");
#endif
  tb = rk->tb->coefficient;
  ac = rk->ac->coefficient;
  k = fminl (0.L, a20 (ac));
  if (a21 (ac) < 0.L)
    k += a21 (ac);
  if (a30 (ac) < 0.L)
    k += a30 (ac);
  if (a31 (ac) < 0.L)
    k += a31 (ac);
  if (a32 (ac) < 0.L)
    k += a32 (ac);
  if (a40 (ac) < 0.L)
    k += a40 (ac);
  if (a41 (ac) < 0.L)
    k += a41 (ac);
  if (a42 (ac) < 0.L)
    k += a42 (ac);
  if (a43 (ac) < 0.L)
    k += a43 (ac);
  if (a50 (ac) < 0.L)
    k += a50 (ac);
  if (a51 (ac) < 0.L)
    k += a51 (ac);
  if (a52 (ac) < 0.L)
    k += a52 (ac);
  if (a53 (ac) < 0.L)
    k += a53 (ac);
  if (a54 (ac) < 0.L)
    k += a54 (ac);
  if (a60 (ac) < 0.L)
    k += a60 (ac);
  if (a61 (ac) < 0.L)
    k += a61 (ac);
  if (a62 (ac) < 0.L)
    k += a62 (ac);
  if (a63 (ac) < 0.L)
    k += a63 (ac);
  if (a64 (ac) < 0.L)
    k += a64 (ac);
  if (a65 (ac) < 0.L)
    k += a65 (ac);
  if (k < 0.L)
    {
      k = 20.L - k;
      goto end;
    }
  k = fminl (0.L, c20 (ac));
  if (c21 (ac) < 0.L)
    k += c21 (ac);
  if (c30 (ac) < 0.L)
    k += c30 (ac);
  if (c31 (ac) < 0.L)
    k += c31 (ac);
  if (c32 (ac) < 0.L)
    k += c32 (ac);
  if (c40 (ac) < 0.L)
    k += c40 (ac);
  if (c41 (ac) < 0.L)
    k += c41 (ac);
  if (c42 (ac) < 0.L)
    k += c42 (ac);
  if (c43 (ac) < 0.L)
    k += c43 (ac);
  if (c50 (ac) < 0.L)
    k += c50 (ac);
  if (c51 (ac) < 0.L)
    k += c51 (ac);
  if (c52 (ac) < 0.L)
    k += c52 (ac);
  if (c53 (ac) < 0.L)
    k += c53 (ac);
  if (c54 (ac) < 0.L)
    k += c54 (ac);
  if (c60 (ac) < 0.L)
    k += c60 (ac);
  if (c61 (ac) < 0.L)
    k += c61 (ac);
  if (c62 (ac) < 0.L)
    k += c62 (ac);
  if (c63 (ac) < 0.L)
    k += c63 (ac);
  if (c64 (ac) < 0.L)
    k += c64 (ac);
  if (c65 (ac) < 0.L)
    k += c65 (ac);
  if (k < 0.L)
    {
      k = 10.L - k;
      goto end;
    }
  k = 1.L / rk_cfl_6 (tb, ac);
end:
#if DEBUG_RK
  fprintf (stderr, "rk_objective ac_6: objective=%Lg\n", k);
  fprintf (stderr, "rk_objective ac_6: end\n");
#endif
  return k;
}

/**
 * Function to init required variables on a RK struct data.
 */
static inline void
rk_init (RK * rk,               ///< RK struct.
         gsl_rng * rng,         ///< GSL pseudo-random number generator struct.
         unsigned int thread)   ///< thread number.
{
  optimize_init (rk->tb, rng, thread);
  if (rk->strong)
    optimize_init (rk->ac0, rng, 0);
}

/**
 * Function to free the memory allocated by a RK struct.
 */
static inline void
rk_delete (RK * rk)             ///< RK struct.
{
  if (rk->strong)
    optimize_delete (rk->ac0);
  optimize_delete (rk->tb);
}

/**
 * Function to perform every optimization step for the a-c Runge-Kutta
 * coefficients.
 */
static inline void
rk_step_ac (RK * rk)            ///< RK struct.
{
  Optimize *tb, *ac;
  long double *is, *vo, *vo2;
  long double o, o2, v, f;
  unsigned long long int ii, nsimulations;
  unsigned int i, j, k, n, nfree;

#if DEBUG_RK
  fprintf (stderr, "rk_step_ac: start\n");
#endif

  // save optimal values
#if DEBUG_RK
  fprintf (stderr, "rk_step_ac: save optimal values\n");
#endif
  tb = rk->tb;
  ac = rk->ac;
  nfree = ac->nfree;
  o2 = INFINITY;
  vo = (long double *) alloca (nfree * sizeof (long double));
  vo2 = (long double *) alloca (nfree * sizeof (long double));
  memcpy (vo, ac->value_optimal, nfree * sizeof (long double));

  // optimzation algorithm sampling
#if DEBUG_RK
  fprintf (stderr, "rk_step_ac: optimization algorithm sampling\n");
  fprintf (stderr, "rk_step_ac: nsimulations=%Lu nclimbings=%u nfree=%u\n",
           ac->nsimulations, ac->nclimbings, ac->nfree);
#endif
  nsimulations = ac->nsimulations;
  for (ii = 0L; ii < nsimulations; ++ii)
    {

      // random freedom degrees
#if DEBUG_RK
      fprintf (stderr, "rk_step_ac: random freedom degrees\n");
#endif
      optimize_generate_freedom (ac, ii);

      // method coefficients
#if DEBUG_RK
      fprintf (stderr, "rk_step_ac: method coefficients\n");
#endif
      if (!ac->method ((Optimize *) rk))
        o = INFINITY;
      else
        o = ac->objective ((Optimize *) rk);
#if DEBUG_RK
      fprintf (stderr, "rk_step_ac: objective=%Lg o2=%Lg\n", o, o2);
#endif
      if (o < o2)
        {
          o2 = o;
          memcpy (vo, ac->random_data, nfree * sizeof (long double));
        }
      if (file_variables)
        {
          g_mutex_lock (mutex);
          print_variables (tb->random_data, tb->nfree, file_variables);
          print_variables (ac->random_data, nfree, file_variables);
          fprintf (file_variables, "%.19Le\n", o);
          g_mutex_unlock (mutex);
        }
    }

  // array of intervals to climb around the optimal
#if DEBUG_RK
  fprintf (stderr,
           "rk_step_ac: array of intervals to climb around the optimal\n");
  fprintf (stderr, "rk_step_ac: nclimbings=%u climbing_factor=%Lg\n",
           ac->nclimbings, ac->climbing_factor);
#endif
  is = (long double *) alloca (nfree * sizeof (long double));
  for (j = 0; j < nfree; ++j)
    is[j] = ac->interval0[j] * ac->climbing_factor;
#if DEBUG_RK
  for (j = 0; j < nfree; ++j)
    fprintf (stderr, "rk_step_ac: i=%u is=%Lg\n", j, is[j]);
#endif

  // hill climbing algorithm bucle
#if DEBUG_RK
  fprintf (stderr, "rk_step_ac: hill climbing algorithm bucle\n");
#endif
  memcpy (vo2, vo, nfree * sizeof (long double));
  memcpy (ac->random_data, vo, nfree * sizeof (long double));
  n = ac->nclimbings;
  for (i = 0; i < n; ++i)
    {
#if DEBUG_RK
      for (j = 0; j < nfree; ++j)
        fprintf (stderr, "rk_step_ac: j=%u is=%Lg\n", j, is[j]);
#endif
      for (j = k = 0; j < nfree; ++j)
        {
          v = vo[j];
          ac->random_data[j] = v + is[j];
#if DEBUG_RK
          fprintf (stderr, "rk_step_ac: j=%u random=%Lg\n", j,
                   ac->random_data[j]);
#endif
          if (!ac->method ((Optimize *) rk))
            o = INFINITY;
          else
            o = ac->objective ((Optimize *) rk);
#if DEBUG_RK
          fprintf (stderr, "rk_step_ac: k=%u objective=%Lg o2=%Lg\n", k, o, o2);
#endif
          if (o < o2)
            {
              k = 1;
              o2 = o;
              memcpy (vo2, ac->random_data, nfree * sizeof (long double));
            }
          if (file_variables)
            {
              g_mutex_lock (mutex);
              print_variables (tb->random_data, tb->nfree, file_variables);
              print_variables (ac->random_data, nfree, file_variables);
              fprintf (file_variables, "%.19Le\n", o);
              g_mutex_unlock (mutex);
            }
          ac->random_data[j] = fmaxl (0.L, v - is[j]);
          if (!ac->method ((Optimize *) rk))
            o = INFINITY;
          else
            o = ac->objective ((Optimize *) rk);
#if DEBUG_RK
          fprintf (stderr, "rk_step_ac: k=%u objective=%Lg o2=%Lg\n", k, o, o2);
#endif
          if (o < o2)
            {
              k = 1;
              o2 = o;
              memcpy (vo2, ac->random_data, nfree * sizeof (long double));
            }
          if (file_variables)
            {
              g_mutex_lock (mutex);
              print_variables (tb->random_data, tb->nfree, file_variables);
              print_variables (ac->random_data, nfree, file_variables);
              fprintf (file_variables, "%.19Le\n", o);
              g_mutex_unlock (mutex);
            }
          ac->random_data[j] = v;
        }

      // update optimal values and increase or reduce intervals if converging or
      // not
      if (!k)
        f = 0.5L;
      else
        {
          f = 1.2L;
          memcpy (vo, vo2, nfree * sizeof (long double));
        }
      for (j = 0; j < nfree; ++j)
        is[j] *= f;
    }

  // update optimal values
#if DEBUG_RK
  fprintf (stderr, "rk_step_ac: update optimal values\n");
  fprintf (stderr, "rk_step_ac: optimal=%Lg o2=%Lg\n", *ac->optimal, o2);
#endif
  if (o2 < *ac->optimal)
    {
      *ac->optimal = o2;
      memcpy (ac->value_optimal, vo2, nfree * sizeof (long double));
    }

#if DEBUG_RK
  fprintf (stderr, "rk_step_ac: end\n");
#endif
}

/**
 * Function to do the optimization bucle for the a-c Runge-Kutta coefficients.
 */
void
rk_bucle_ac (RK * rk)           ///< RK struct.
{
  Optimize *tb, *ac, *ac0;
  long double *vo;
  long double optimal;
  unsigned int i, nfree;

#if DEBUG_RK
  fprintf (stderr, "rk_bucle_ac: start\n");
#endif

  tb = rk->tb;
  ac = rk->ac;
  ac0 = rk->ac0;
  nfree = ac0->nfree;
  vo = (long double *) alloca (nfree * sizeof (long double));

  // Init some parameters
#if DEBUG_RK
  fprintf (stderr, "rk_bucle_ac: nfree=%u optimal=%Lg\n", nfree, *tb->optimal);
#endif
  *ac0->optimal = optimal = *tb->optimal;
  for (i = 0; i < nfree; ++i)
    vo[i] = ac0->minimum[i] + 0.5L * ac0->interval[i];
  memcpy (ac, ac0, sizeof (Optimize));
  ac->optimal = &optimal;
  ac->value_optimal = vo;
  optimize_init (ac, ac0->rng, 0);
#if DEBUG_RK
  for (i = 0; i < tb->nfree; ++i)
    fprintf (stderr, "rk_bucle_ac: i=%u random=%Lg\n", i, tb->random_data[i]);
  for (i = 0; i < nfree; ++i)
    fprintf (stderr, "rk_bucle_ac: i=%u minimum=%Lg interval=%Lg type=%u\n",
             i, ac->minimum[i], ac->interval[i], ac->random_type[i]);
#endif

  // Iterate
#if DEBUG_RK
  fprintf (stderr, "rk_bucle_ac: iterate\n");
#endif
  for (i = 0; i < ac->niterations; ++i)
    {

      // Optimization step
      rk_step_ac (rk);

      // Updating coefficient intervals to converge
      optimize_converge (ac);

      // Iterate
#if DEBUG_RK
      fprintf (stderr, "Iteration ac %u\n", i);
#endif
    }

  // Check and save optimal
  if (optimal < *ac0->optimal)
    {
#if DEBUG_RK
      fprintf (stderr, "rk_bucle_ac: optimal=%Lg\n", *ac0->optimal);
#endif
      g_mutex_lock (mutex);
      *ac0->optimal = optimal;
      memcpy (ac0->value_optimal, vo, nfree * sizeof (long double));
      g_mutex_unlock (mutex);
#if DEBUG_RK
      fprintf (stderr, "rk_bucle_ac: optimal=%Lg\n", *ac0->optimal);
      for (i = 0; i < ac0->nfree; ++i)
        fprintf (stderr, "rk_bucle_ac: vo%u=%Lg\n", i, ac0->value_optimal[i]);
#endif
    }

  // Free memory
  optimize_delete (ac);

#if DEBUG_RK
  fprintf (stderr, "rk_bucle_ac: end\n");
#endif
}

/**
 * Function to perform every optimization step for the t-b Runge-Kutta
 * coefficients.
 */
static inline void
rk_step_tb (RK * rk)            ///< RK struct.
{
  Optimize *tb;
  long double *is, *vo;
  long double o, v, f;
  unsigned long long int ii, nrandom;
  unsigned int b, i, j, k, n, nfree;

#if DEBUG_RK
  fprintf (stderr, "rk_step_tb: start\n");
#endif

  // save optimal values
#if DEBUG_RK
  fprintf (stderr, "rk_step_tb: save optimal values\n");
#endif
  tb = rk->tb;
  nfree = tb->nfree;
  vo = (long double *) alloca (nfree * sizeof (long double));
  b = (file_variables && !rk->strong) ? 1 : 0;

  // optimization algorithm sampling
#if DEBUG_RK
  fprintf (stderr, "rk_step_tb: optimization algorithm sampling\n");
  fprintf (stderr, "rk_step_tb: nsimulations=%Lu nclimbings=%u\n",
           tb->nsimulations, tb->nclimbings);
#endif
  ii = tb->nsimulations * (rank * nthreads + tb->thread) / (nnodes * nthreads);
  nrandom = tb->nsimulations * (rank * nthreads + tb->thread + 1)
    / (nnodes * nthreads);
  for (; ii < nrandom; ++ii)
    {

      // random freedom degrees
#if DEBUG_RK
      fprintf (stderr, "rk_step_tb: random freedom degrees\n");
#endif
      optimize_generate_freedom (tb, ii);

      // method coefficients
#if DEBUG_RK
      fprintf (stderr, "rk_step_tb: method coefficients\n");
#endif
      if (!tb->method (tb))
        o = INFINITY;
      else
        o = tb->objective (tb);
      if (o < *tb->optimal)
        {
          g_mutex_lock (mutex);
          *tb->optimal = o;
          memcpy (tb->value_optimal, tb->random_data,
                  nfree * sizeof (long double));
          g_mutex_unlock (mutex);
        }
      if (b)
        {
          g_mutex_lock (mutex);
          print_variables (tb->random_data, nfree, file_variables);
          fprintf (file_variables, "%.19Le\n", o);
          g_mutex_unlock (mutex);
        }
    }

  // array of intervals to climb around the optimal
#if DEBUG_RK
  fprintf (stderr,
           "rk_step_tb: array of intervals to climb around the optimal\n");
#endif
  is = (long double *) alloca (nfree * sizeof (long double));
  for (j = 0; j < nfree; ++j)
    is[j] = tb->interval0[j] * tb->climbing_factor;

  // hill climbing algorithm bucle
#if DEBUG_RK
  fprintf (stderr, "rk_step_tb: hill climbing algorithm bucle\n");
#endif
  memcpy (tb->random_data, tb->value_optimal, nfree * sizeof (long double));
  memcpy (vo, tb->value_optimal, nfree * sizeof (long double));
  n = tb->nclimbings;
  for (i = 0; i < n; ++i)
    {
      for (j = k = 0; j < nfree; ++j)
        {
          v = vo[j];
          tb->random_data[j] = v + is[j];
          if (!tb->method (tb))
            o = INFINITY;
          else
            o = tb->objective (tb);
          if (o < *tb->optimal)
            {
              k = 1;
              g_mutex_lock (mutex);
              *tb->optimal = o;
              memcpy (tb->value_optimal, tb->random_data,
                      nfree * sizeof (long double));
              g_mutex_unlock (mutex);
            }
          if (b)
            {
              g_mutex_lock (mutex);
              print_variables (tb->random_data, nfree, file_variables);
              fprintf (file_variables, "%.19Le\n", o);
              g_mutex_unlock (mutex);
            }
          tb->random_data[j] = fmaxl (0.L, v - is[j]);
          if (!tb->method (tb))
            o = INFINITY;
          else
            o = tb->objective (tb);
          if (o < *tb->optimal)
            {
              k = 1;
              g_mutex_lock (mutex);
              *tb->optimal = o;
              memcpy (tb->value_optimal, tb->random_data,
                      nfree * sizeof (long double));
              g_mutex_unlock (mutex);
            }
          if (b)
            {
              g_mutex_lock (mutex);
              print_variables (tb->random_data, nfree, file_variables);
              fprintf (file_variables, "%.19Le\n", o);
              g_mutex_unlock (mutex);
            }
          tb->random_data[j] = v;
        }

      // increase or reduce intervals if converging or not
      if (!k)
        f = 0.5L;
      else
        {
          f = 1.2L;
          memcpy (vo, tb->value_optimal, nfree * sizeof (long double));
        }
      for (j = 0; j < nfree; ++j)
        is[j] *= f;
    }
#if DEBUG_RK
  fprintf (stderr, "rk_step_tb: end\n");
#endif
}

/**
 * Function to do the optimization bucle.
 */
static inline void
rk_bucle_tb (RK * rk)           ///< RK struct.
{
  GThread *thread[nthreads];
  Optimize *tb, *ac;
#if HAVE_MPI
  long double *vo;
  MPI_Status status;
#endif
  unsigned int i, j, nfree, nfree2, strong;

#if DEBUG_RK
  fprintf (stderr, "rk_bucle_tb: start\n");
#endif

  // Allocate local array of optimal values
#if DEBUG_RK
  fprintf (stderr, "rk_bucle_tb: allocate local array of optimal values\n");
#endif
  tb = rk->tb;
  nfree = tb->nfree;
  strong = rk->strong;
  if (strong)
    {
      ac = rk->ac0;
      nfree2 = ac->nfree;
    }
  else
    nfree2 = 0;
#if HAVE_MPI
  vo = (long double *) alloca ((1 + nfree + nfree2) * sizeof (long double));
#endif

  // Init some parameters
#if DEBUG_RK
  fprintf (stderr, "rk_bucle_tb: init some parameters\n");
  fprintf (stderr, "rk_bucle_tb: nfree=%u\n", nfree);
#endif
  *tb->optimal = INFINITY;
  for (i = 0; i < nfree; ++i)
    tb->value_optimal[i] = tb->minimum[i] + 0.5L * tb->interval[i];
  if (strong)
    for (i = 0; i < nfree2; ++i)
      ac->value_optimal[i] = ac->minimum[i] + 0.5L * ac->interval[i];

  // Iterate
#if DEBUG_RK
  fprintf (stderr, "rk_bucle_tb: iterate\n");
#endif
  for (i = 0; i < tb->niterations; ++i)
    {

      // Optimization step parallelized for every node by GThreads
      if (nthreads > 1)
        {
          for (j = 0; j < nthreads; ++j)
            thread[j]
              = g_thread_new (NULL,
                              (GThreadFunc) (void (*)(void)) rk_step_tb,
                              (void *) (rk + j));
          for (j = 0; j < nthreads; ++j)
            g_thread_join (thread[j]);
        }
      else
        rk_step_tb (rk);

#if HAVE_MPI
      if (rank > 0)
        {

          // Secondary nodes send the optimal coefficients to the master node
          vo[0] = *tb->optimal;
          memcpy (vo + 1, tb->value_optimal, nfree * sizeof (long double));
          if (strong)
            memcpy (vo + 1 + nfree, ac->value_optimal,
                    nfree2 * sizeof (long double));
          MPI_Send (vo, 1 + nfree + nfree2, MPI_LONG_DOUBLE, 0, 1,
                    MPI_COMM_WORLD);

          // Secondary nodes receive the optimal coefficients
          MPI_Recv (vo, 1 + nfree + nfree2, MPI_LONG_DOUBLE, 0, 1,
                    MPI_COMM_WORLD, &status);
          *tb->optimal = *ac->optimal = vo[0];
          memcpy (tb->value_optimal, vo + 1, nfree * sizeof (long double));
          if (strong)
            memcpy (ac->value_optimal, vo + 1 + nfree,
                    nfree2 * sizeof (long double));
        }
      else
        {
          printf ("rank=%d optimal=%.19Le\n", rank, *tb->optimal);

          for (j = 1; j < nnodes; ++j)
            {

              // Master node receives the optimal coefficients obtained by
              // secondary nodes
              MPI_Recv (vo, 1 + nfree + nfree2, MPI_LONG_DOUBLE, j, 1,
                        MPI_COMM_WORLD, &status);

              // Master node selects the optimal coefficients
              if (vo[0] < *tb->optimal)
                {
                  *tb->optimal = *ac->optimal = vo[0];
                  memcpy (tb->value_optimal, vo + 1,
                          nfree * sizeof (long double));
                  if (strong)
                    memcpy (ac->value_optimal, vo + 1 + nfree,
                            nfree2 * sizeof (long double));
                }
            }

          // Master node sends the optimal coefficients to secondary nodes
          vo[0] = *tb->optimal;
          memcpy (vo + 1, tb->value_optimal, nfree * sizeof (long double));
          if (strong)
            memcpy (vo + 1 + nfree, ac->value_optimal,
                    nfree2 * sizeof (long double));
          for (j = 1; j < nnodes; ++i)
            MPI_Send (vo, 1 + nfree + nfree2, MPI_LONG_DOUBLE, j, 1,
                      MPI_COMM_WORLD);
        }

#endif

      // Print the optimal coefficients
#if DEBUG_RK
      optimize_print_random (tb, stderr);
      if (strong)
        optimize_print_random (ac, stderr);
      fprintf (stderr, "optimal=%.19Le\n", *tb->optimal);
#endif

      // Updating coefficient intervals to converge
      optimize_converge (tb);

      // Iterate
      printf ("Iteration %u Optimal %.19Le\n", i, *tb->optimal);
    }

#if DEBUG_RK
  fprintf (stderr, "rk_bucle_tb: end\n");
#endif
}

/**
 * Function to select the Runge-Kutta method.
 *
 * \return 1 on success, 0 on error.
 */
static inline int
rk_select (RK * rk,             ///< RK struct.
           unsigned int nsteps, ///< steps number.
           unsigned int order)  ///< accuracy order.
{
  static int (*tb_method[7][6]) (Optimize *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_2_2, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_3_2, &rk_tb_3_3, NULL, NULL},
    {
    NULL, NULL, &rk_tb_4_2, &rk_tb_4_3, &rk_tb_4_4, NULL},
    {
    NULL, NULL, &rk_tb_5_2, &rk_tb_5_3, &rk_tb_5_4, NULL},
    {
    NULL, NULL, &rk_tb_6_2, &rk_tb_6_3, &rk_tb_6_4, NULL}
  };
  static int (*tb_method_t[7][6]) (Optimize *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_2_2t, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_3_2t, &rk_tb_3_3t, NULL, NULL},
    {
    NULL, NULL, &rk_tb_4_2t, &rk_tb_4_3t, &rk_tb_4_4t, NULL},
    {
    NULL, NULL, &rk_tb_5_2t, &rk_tb_5_3t, &rk_tb_5_4t, NULL},
    {
    NULL, NULL, &rk_tb_6_2t, &rk_tb_6_3t, &rk_tb_6_4t, NULL}
  };
  static int (*tb_method_p[7][6]) (Optimize *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_2_2p, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_3_2p, &rk_tb_3_3p, NULL, NULL},
    {
    NULL, NULL, &rk_tb_4_2p, &rk_tb_4_3p, NULL, NULL},
    {
    NULL, NULL, &rk_tb_5_2p, &rk_tb_5_3p, &rk_tb_5_4p, NULL},
    {
    NULL, NULL, &rk_tb_6_2p, &rk_tb_6_3p, &rk_tb_6_4p, NULL}
  };
  static int (*tb_method_tp[7][6]) (Optimize *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_2_2tp, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_tb_3_2tp, &rk_tb_3_3tp, NULL, NULL},
    {
    NULL, NULL, &rk_tb_4_2tp, &rk_tb_4_3tp, NULL, NULL},
    {
    NULL, NULL, &rk_tb_5_2tp, &rk_tb_5_3tp, &rk_tb_5_4tp, NULL},
    {
    NULL, NULL, &rk_tb_6_2tp, &rk_tb_6_3tp, &rk_tb_6_4tp, NULL}
  };
  static long double (*tb_objective[7][6]) (RK *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_2_2, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_2_2, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_3_2, &rk_objective_tb_3_3, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_4_2, &rk_objective_tb_4_3,
        &rk_objective_tb_4_4, NULL},
    {
    NULL, NULL, &rk_objective_tb_5_2, &rk_objective_tb_5_3,
        &rk_objective_tb_5_4, NULL},
    {
    NULL, NULL, &rk_objective_tb_6_2, &rk_objective_tb_6_3,
        &rk_objective_tb_6_4, NULL}
  };
  static long double (*tb_objective_t[7][6]) (RK *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_2_2t, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_2_2t, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_3_2t, &rk_objective_tb_3_3t, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_4_2t, &rk_objective_tb_4_3t,
        &rk_objective_tb_4_4t, NULL},
    {
    NULL, NULL, &rk_objective_tb_5_2t, &rk_objective_tb_5_3t,
        &rk_objective_tb_5_4t, NULL},
    {
    NULL, NULL, &rk_objective_tb_6_2t, &rk_objective_tb_6_3t,
        &rk_objective_tb_6_4t, NULL}
  };
  static long double (*tb_objective_p[7][6]) (RK *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_2_2, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_3_2, &rk_objective_tb_3_3p, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_4_2, &rk_objective_tb_4_3p, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_5_2, &rk_objective_tb_5_3p,
        &rk_objective_tb_5_4p, NULL},
    {
    NULL, NULL, &rk_objective_tb_6_2, &rk_objective_tb_6_3p,
        &rk_objective_tb_6_4p, NULL}
  };
  static long double (*tb_objective_tp[7][6]) (RK *) =
  {
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, NULL, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_2_2t, NULL, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_3_2t, &rk_objective_tb_3_3tp, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_4_2t, &rk_objective_tb_4_3tp, NULL, NULL},
    {
    NULL, NULL, &rk_objective_tb_5_2t, &rk_objective_tb_5_3tp,
        &rk_objective_tb_5_4tp, NULL},
    {
    NULL, NULL, &rk_objective_tb_6_2t, &rk_objective_tb_6_3tp,
        &rk_objective_tb_6_4tp, NULL}
  };
  static int (*ac_method[7]) (RK *) =
  {
  NULL, NULL, &rk_ac_2, &rk_ac_3, &rk_ac_4, &rk_ac_5, &rk_ac_6};
  static long double (*ac_objective[7]) (RK * rk) =
  {
  NULL, NULL, &rk_objective_ac_2, &rk_objective_ac_3, &rk_objective_ac_4,
      &rk_objective_ac_5, &rk_objective_ac_6};
  const unsigned int nequations[6] = { 0, 1, 2, 4, 8, 16 };
  Optimize *tb, *ac;
#if DEBUG_RK
  fprintf (stderr, "rk_select: start\n");
#endif
  tb = rk->tb;
  ac = rk->ac;
  tb->nsteps = nsteps;
  tb->order = order;
  tb->size = nsteps * (nsteps + 3) / 2 - 1;
  tb->nfree = tb->size - nsteps + 1 - nequations[order];
  if (rk->pair)
    {
      tb->size += nsteps - 1;
      if (rk->time_accuracy)
        {
          tb->method = tb_method_tp[nsteps][order];
          tb->objective = (OptimizeObjective) tb_objective_tp[nsteps][order];
        }
      else
        {
          tb->method = tb_method_p[nsteps][order];
          tb->objective = (OptimizeObjective) tb_objective_p[nsteps][order];
        }
    }
  else
    {
      if (rk->time_accuracy)
        {
          tb->method = tb_method_t[nsteps][order];
          tb->objective = (OptimizeObjective) tb_objective_t[nsteps][order];
        }
      else
        {
          tb->method = tb_method[nsteps][order];
          tb->objective = (OptimizeObjective) tb_objective[nsteps][order];
        }
    }
  if (!tb->method)
    goto exit_on_error;
  if (rk->time_accuracy)
    {
      --tb->nfree;
      if (rk->pair)
        switch (nsteps)
          {
          case 5:
            switch (order)
              {
              case 4:
                --tb->nfree;
              }
            break;
          case 6:
            switch (order)
              {
              case 4:
                --tb->nfree;
              }
            break;
          }
    }
  tb->minimum0
    = (long double *) g_slice_alloc (tb->nfree * sizeof (long double));
  tb->interval0
    = (long double *) g_slice_alloc (tb->nfree * sizeof (long double));
  tb->random_type
    = (unsigned int *) g_slice_alloc (tb->nfree * sizeof (unsigned int));
  if (rk->strong)
    {
      ac = rk->ac0;
      ac->size = nsteps * (nsteps + 1) - 2;
      ac->nfree = nsteps * (nsteps - 1) / 2;
      ac->minimum0
        = (long double *) g_slice_alloc (ac->nfree * sizeof (long double));
      ac->interval0
        = (long double *) g_slice_alloc (ac->nfree * sizeof (long double));
      ac->random_type
        = (unsigned int *) g_slice_alloc (ac->nfree * sizeof (unsigned int));
      ac->method = (OptimizeMethod) ac_method[nsteps];
      if (!ac->method)
        goto exit_on_error;
      ac->objective = (OptimizeObjective) ac_objective[nsteps];
    }

#if DEBUG_RK
  fprintf (stderr, "rk_select: end\n");
#endif
  return 1;

exit_on_error:
  error_message = g_strdup (_("Unknown method"));
#if DEBUG_RK
  fprintf (stderr, "rk_select: end\n");
#endif
  return 0;
}

/**
 * Function to read the Runge-Kutta method data on a XML node.
 *
 * \return 1 on success, 0 on error.
 */
int
rk_run (xmlNode * node,         ///< XML node.
        gsl_rng ** rng)         ///< array of gsl_rng structs.
{
  RK rk[nthreads];
  char filename[64];
  Optimize *tb, *ac;
  gchar *buffer;
  xmlChar *prop;
  FILE *file;
  long double *value_optimal, *value_optimal2;
  long double optimal, optimal2;
  int code;
  unsigned int i, j, nsteps, order, nfree, nfree2;

#if DEBUG_RK
  fprintf (stderr, "rk_run: start\n");
#endif

  tb = rk->tb;
  nsteps = xml_node_get_uint (node, XML_STEPS, &code);
  if (code)
    {
      error_message = g_strdup (_("Bad steps number"));
      goto exit_on_error;
    }
  order = xml_node_get_uint (node, XML_ORDER, &code);
  if (code)
    {
      error_message = g_strdup (_("Bad order"));
      goto exit_on_error;
    }
  prop = xmlGetProp (node, XML_STRONG);
  if (!prop || !xmlStrcmp (prop, XML_NO))
    rk->strong = 0;
  else if (!xmlStrcmp (prop, XML_YES))
    rk->strong = 1;
  else
    {
      error_message = g_strdup (_("Bad strong stability"));
      goto exit_on_error;
    }
  xmlFree (prop);
  prop = xmlGetProp (node, XML_PAIR);
  if (!prop || !xmlStrcmp (prop, XML_NO))
    rk->pair = 0;
  else if (!xmlStrcmp (prop, XML_YES))
    rk->pair = 1;
  else
    {
      error_message = g_strdup (_("Bad pair"));
      goto exit_on_error;
    }
  xmlFree (prop);
  prop = xmlGetProp (node, XML_TIME_ACCURACY);
  if (!prop || !xmlStrcmp (prop, XML_NO))
    rk->time_accuracy = 0;
  else if (!xmlStrcmp (prop, XML_YES))
    rk->time_accuracy = 1;
  else
    {
      error_message = g_strdup (_("Bad time accuracy"));
      goto exit_on_error;
    }
  xmlFree (prop);
  if (!rk_select (rk, nsteps, order))
    goto exit_on_error;
  if (!optimize_read (tb, node))
    goto exit_on_error;
  nfree = tb->nfree;
  value_optimal = (long double *) g_slice_alloc (nfree * sizeof (long double));
  optimize_create (tb, &optimal, value_optimal);
  node = node->children;
  for (i = 0; i < nfree; ++i, node = node->next)
    if (!read_variable (node, tb->minimum0, tb->interval0, tb->random_type, i))
      goto exit_on_error;
  if (rk->strong)
    {
      ac = rk->ac0;
      if (!node)
        {
          error_message = g_strdup (_("No a-c coefficients data"));
          goto exit_on_error;
        }
      if (xmlStrcmp (node->name, XML_AC))
        {
          error_message = g_strdup (_("Bad a-c coefficients XML node"));
          goto exit_on_error;
        }
      if (!optimize_read (ac, node))
        {
          buffer = error_message;
          error_message
            = g_strconcat (_("a-c coefficients"), ":\n", error_message, NULL);
          g_free (buffer);
          goto exit_on_error;
        }
      nfree2 = ac->nfree;
      value_optimal2
        = (long double *) g_slice_alloc (nfree2 * sizeof (long double));
      optimize_create (ac, &optimal2, value_optimal2);
      for (i = 0; i < nfree2; ++i)
        {
          node = node->next;
          if (!read_variable (node, ac->minimum0, ac->interval0,
                              ac->random_type, i))
            goto exit_on_error;
        }
    }
  for (i = 1; i < nthreads; ++i)
    memcpy (rk + i, rk, sizeof (RK));
  j = rank * nthreads;
  for (i = 0; i < nthreads; ++i)
    rk_init (rk + i, rng[j + i], i);

  // Method bucle
  printf ("Optimize bucle\n");
  rk_bucle_tb (rk);

  // Print the optimal coefficients
  printf ("Print the optimal coefficients\n");
  memcpy (tb->random_data, tb->value_optimal, nfree * sizeof (long double));
  code = tb->method (tb);
  if (rk->strong)
    {
      memcpy (ac->random_data, ac->value_optimal,
              nfree2 * sizeof (long double));
      memcpy (rk->ac, ac, sizeof (Optimize));
      code = ac->method ((Optimize *) rk);
    }
  snprintf (filename, 64, "rk-%u-%u-%u-%u-%u.mc",
            nsteps, order, rk->time_accuracy, rk->pair, rk->strong);
  file = fopen (filename, "w");
  print_maxima_precision (file);
  rk_print (rk, file);
  rk_print_maxima (file, nsteps, nsteps, order, 'b');
  if (rk->pair)
    rk_print_maxima (file, nsteps, nsteps - 1, order - 1, 'e');
  if (rk->strong)
    ac_print_maxima (file, nsteps);
  fclose (file);
  snprintf (filename, 64, "sed -i 's/e+/b+/g' rk-%u-%u-%u-%u-%u.mc",
            nsteps, order, rk->time_accuracy, rk->pair, rk->strong);
  code = system (filename);
  snprintf (filename, 64, "sed -i 's/e-/b-/g' rk-%u-%u-%u-%u-%u.mc",
            nsteps, order, rk->time_accuracy, rk->pair, rk->strong);
  code = system (filename);

  // Free memory
  if (rk->strong)
    {
      g_slice_free1 (nfree2 * sizeof (unsigned int), ac->random_type);
      g_slice_free1 (nfree2 * sizeof (long double), ac->interval0);
      g_slice_free1 (nfree2 * sizeof (long double), ac->minimum0);
      g_slice_free1 (nfree2 * sizeof (long double), value_optimal2);
    }
  for (i = 0; i < nthreads; ++i)
    rk_delete (rk + i);
  g_slice_free1 (nfree * sizeof (unsigned int), tb->random_type);
  g_slice_free1 (nfree * sizeof (long double), tb->interval0);
  g_slice_free1 (nfree * sizeof (long double), tb->minimum0);
  g_slice_free1 (nfree * sizeof (long double), value_optimal);

#if DEBUG_RK
  fprintf (stderr, "rk_run: end\n");
#endif
  return 1;

exit_on_error:
  buffer = error_message;
  error_message = g_strconcat ("Runge-Kutta:\n", buffer, NULL);
  g_free (buffer);
#if DEBUG_RK
  fprintf (stderr, "rk_run: end\n");
#endif
  return 0;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
// gcc -msse2 -O3 -ftree-vectorize test_dcl.c dclhelpers.c os_generic.c -DFLT=double -lpthread -lcblas && valgrind
// ./a.out

#include "dclhelpers.h"
#include "os_generic.h"
#include <assert.h>
#include <cblas.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

#include "sv_matrix.h"

void fill_random(FLT *A, int ld, int m, int n) {
	assert(ld == n);
	for (int y = 0; y < m; y++)
		for (int x = 0; x < n; x++)
			A[y * ld + x] = (rand()) / (double)RAND_MAX;
}

void test_dcldgemm_speed(const char *name, char transA, char transB, int m, int n, int k, DCL_FLOAT alpha,
						 const DCL_FLOAT *A, int Ac, const DCL_FLOAT *B, int Bc, DCL_FLOAT beta) {
	printf("%s speed test:\n", name);
	double times[2];
	FLT emo[2][m][n];

	for (int z = 0; z < 2; z++) {
		double start = OGGetAbsoluteTime();
		for (int i = 0; i < 100000; i++) {
			dclZero(DMS(emo[z]), m, n);

			if (z == 0) {
				dcldgemm(transA, transB, m, n, k, alpha, A, Ac, B, Bc, beta, DMS(emo[z]));
			} else {

				cblas_dgemm(CblasRowMajor, transA == 1 ? CblasTrans : CblasNoTrans,
							transB == 1 ? CblasTrans : CblasNoTrans, m, n, k, alpha, A, Ac, B, Bc, beta, emo[z][0], n);
			}
		}

		times[z] = OGGetAbsoluteTime() - start;
	}

	dclPrint(DMS(emo[0]), m, n);
	dclPrint(DMS(emo[1]), m, n);
	printf("dcl  Elapsed: %f\n", times[0]);
	printf("cblas Elapsed: %f\n", times[1]);
	printf("%fx difference\n", times[0] / times[1]);

	for (int i = 0; i < m; i++) {
		for (int j = 0; j < k; j++) {
			assert(fabs(emo[0][i][j] - emo[1][i][j]) < .00001);
		}
	}
}

void compareToCblas() {
	srand(0);
	int m = 12;
	int n = 20;
	int k = 20;

	FLT em1[m][n];
	FLT em2[n][k];

	fill_random(DMS(em1), m, n);
	fill_random(DMS(em2), n, k);

	dclPrint(DMS(em1), m, n);
	dclPrint(DMS(em2), n, k);

	test_dcldgemm_speed("Simple", 0, 0, m, n, k, 1.0, DMS(em1), DMS(em2), .1);
}

void compareToCblasTrans() {
	srand(0);
	int m = 12;
	int n = 20;
	int k = n;

	FLT em1[m][n];

	fill_random(DMS(em1), m, n);

	dclPrint(DMS(em1), m, n);

	SvMat Em1 = cvMat(m, n, SV_FLT, em1);
	FLT em1tem1[n][n];
	SvMat Em1tEm1 = cvMat(n, n, SV_FLT, em1tem1);
	cvMulTransposed(&Em1, &Em1tEm1, 1, 0, 1);
	print_mat(&Em1tEm1);

	test_dcldgemm_speed("Trans", 0, 0,
						n, // # of rows in OP(A) == em1' -- 20
						n, // # of cols in OP(B) == em1 -- 20
						m, // # of cols in OP(A) == em1' -- 12
						1.0,
						DMS(em1), // Note that LD stays the same
						DMS(em1), 0);
}

int main() {
	FLT A[2][4] = {{0, 1, 2, 3}, {4, 5, 6, 7}};
	FLT B[4][2];
	dclPrint(A[0], 4, 2, 4);
	dclTransp(B[0], 2, A[0], 4, 2, 4);
	dclPrint(B[0], 2, 4, 2);

	int i, j;
	for (i = 0; i < 8; i++) {
		printf("%f\n", ((float *)(B[0]))[i]);
	}

	FLT M[3][3] = {{.32, 1, 0}, {0, 1, 2}, {1, 0, 1}};
	FLT Mo[3][3];
	dclInv(Mo[0], 3, M[0], 3, 3);
	dclPrint(Mo[0], 3, 3, 3);

	FLT MM[3][3];
	dclMul(MM[0], 3, M[0], 3, Mo[0], 3, 3, 3, 3);

	printf("The following should be an identity matrix\n");
	dclPrint(MM[0], 3, 3, 3);

	{
		FLT A[3][4];
		dclIdentity(DMS(A), 3, 4);
		dclPrint(DMS(A), 3, 4);

		FLT x[4][2] = {
			{7, -7}, {8, -8}, {9, -9}, {10, -10},
		};
		FLT R[4][2];
		dclZero(DMS(R), 4, 2);

		// dclMul(R, 1, A[0], 4, x, 1, 4, 1, 3);
		dcldgemm(0, 0, 3, 4, 2, 1, A[0], 4, x[0], 2, 0, R[0], 2);

		dclPrint(DMS(x), 4, 2);
		dclPrint(DMS(R), 3, 2);

		for (int j = 0; j < 2; j++) {
			for (int i = 0; i < 3; i++) {
				printf("[%d][%d]\n", i, j);
				assert(R[i][j] == x[i][j]);
			}

			assert(fabs(R[3][j]) < .0000001);
		}
	}

	compareToCblas();
	compareToCblasTrans();
}

ext: c

max_stars_repo_name: jdavidberger/libsurvive

max_issues_repo_name: jdavidberger/libsurvive

max_forks_repo_name: jdavidberger/libsurvive


text:
/* ode-initval/rkck.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Runge-Kutta 4(5), Cash-Karp */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

/* Cash-Karp constants */
static const double ah[] = { 1.0 / 5.0, 0.3, 3.0 / 5.0, 1.0, 7.0 / 8.0 };
static const double b21 = 1.0 / 5.0;
static const double b3[] = { 3.0 / 40.0, 9.0 / 40.0 };
static const double b4[] = { 0.3, -0.9, 1.2 };
static const double b5[] = { -11.0 / 54.0, 2.5, -70.0 / 27.0, 35.0 / 27.0 };
static const double b6[] =
  { 1631.0 / 55296.0, 175.0 / 512.0, 575.0 / 13824.0, 44275.0 / 110592.0,
    253.0 / 4096.0 };
static const double c1 = 37.0 / 378.0;
static const double c3 = 250.0 / 621.0;
static const double c4 = 125.0 / 594.0;
static const double c6 = 512.0 / 1771.0;
static const double ec[] = { 0.0,
  /* the first value is the same as c1, above */
  37.0 / 378.0 - 2825.0 / 27648.0,
  0.0,
  /* the first value is the same as c3, above */
  250.0 / 621.0 - 18575.0 / 48384.0,
  /* the first value is the same as c4, above */
  125.0 / 594.0 - 13525.0 / 55296.0,
  -277.00 / 14336.0,
  /* the first value is the same as c6, above */
  512.0 / 1771.0 - 0.25
};

typedef struct
{
  double *k1;
  double *k2;
  double *k3;
  double *k4;
  double *k5;
  double *k6;
  double *y0;
  double *ytmp;
}
rkck_state_t;

static void *
rkck_alloc (size_t dim)
{
  rkck_state_t *state = (rkck_state_t *) malloc (sizeof (rkck_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rkck_state", GSL_ENOMEM);
    }

  state->k1 = (double *) malloc (dim * sizeof (double));

  if (state->k1 == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
    }

  state->k2 = (double *) malloc (dim * sizeof (double));

  if (state->k2 == 0)
    {
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
    }

  state->k3 = (double *) malloc (dim * sizeof (double));

  if (state->k3 == 0)
    {
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k3", GSL_ENOMEM);
    }

  state->k4 = (double *) malloc (dim * sizeof (double));

  if (state->k4 == 0)
    {
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k4", GSL_ENOMEM);
    }

  state->k5 = (double *) malloc (dim * sizeof (double));

  if (state->k5 == 0)
    {
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k5", GSL_ENOMEM);
    }

  state->k6 = (double *) malloc (dim * sizeof (double));

  if (state->k6 == 0)
    {
      free (state->k5);
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k6", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->k6);
      free (state->k5);
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->y0);
      free (state->k6);
      free (state->k5);
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  return state;
}

static int
rkck_apply (void *vstate,
            size_t dim,
            double t,
            double h,
            double y[],
            double yerr[],
            const double dydt_in[],
            double dydt_out[], const gsl_odeiv_system * sys)
{
  rkck_state_t *state = (rkck_state_t *) vstate;

  size_t i;
  int status = 0;

  double *const k1 = state->k1;
  double *const k2 = state->k2;
  double *const k3 = state->k3;
  double *const k4 = state->k4;
  double *const k5 = state->k5;
  double *const k6 = state->k6;
  double *const ytmp = state->ytmp;

  /* k1 step */
  if (dydt_in != NULL)
    {
      DBL_MEMCPY (k1, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
      GSL_STATUS_UPDATE (&status, s);
    }

  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + b21 * h * k1[i];

  /* k2 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[0] * h, ytmp, k2);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b3[0] * k1[i] + b3[1] * k2[i]);

  /* k3 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[1] * h, ytmp, k3);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b4[0] * k1[i] + b4[1] * k2[i] + b4[2] * k3[i]);

  /* k4 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[2] * h, ytmp, k4);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b5[0] * k1[i] + b5[1] * k2[i] + b5[2] * k3[i] +
                  b5[3] * k4[i]);

  /* k5 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[3] * h, ytmp, k5);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b6[0] * k1[i] + b6[1] * k2[i] + b6[2] * k3[i] +
                  b6[3] * k4[i] + b6[4] * k5[i]);

  /* k6 step and final sum */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[4] * h, ytmp, k6);
    GSL_STATUS_UPDATE (&status, s);
  }

  for (i = 0; i < dim; i++)
    {
      const double d_i = c1 * k1[i] + c3 * k3[i] + c4 * k4[i] + c6 * k6[i];
      y[i] += h * d_i;
      if (dydt_out != NULL)
        dydt_out[i] = d_i;
    }

  /* difference between 4th and 5th order */
  for (i = 0; i < dim; i++)
    yerr[i] =
      h * (ec[1] * k1[i] + ec[3] * k3[i] + ec[4] * k4[i] + ec[5] * k5[i] +
           ec[6] * k6[i]);

  return status;
}

static int
rkck_reset (void *vstate, size_t dim)
{
  rkck_state_t *state = (rkck_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k1, dim);
  DBL_ZERO_MEMSET (state->k2, dim);
  DBL_ZERO_MEMSET (state->k3, dim);
  DBL_ZERO_MEMSET (state->k4, dim);
  DBL_ZERO_MEMSET (state->k5, dim);
  DBL_ZERO_MEMSET (state->k6, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);

  return GSL_SUCCESS;
}

static unsigned int
rkck_order (void *vstate)
{
  rkck_state_t *state = (rkck_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 5;
}

static void
rkck_free (void *vstate)
{
  rkck_state_t *state = (rkck_state_t *) vstate;

  free (state->ytmp);
  free (state->y0);
  free (state->k6);
  free (state->k5);
  free (state->k4);
  free (state->k3);
  free (state->k2);
  free (state->k1);
  free (state);
}

static const gsl_odeiv_step_type rkck_type = { "rkck",  /* name */
  1,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &rkck_alloc,
  &rkck_apply,
  &rkck_reset,
  &rkck_order,
  &rkck_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rkck = &rkck_type;


text:
/** @file */

#ifndef __CCL_LSST_SPECS_INCLUDED__
#define __CCL_LSST_SPECS_INCLUDED__

#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>

/**
 * We assume no relevant source redshift distributions exceed z= 5
 */
#define Z_MIN_SOURCES 0.
#define Z_MAX_SOURCES 5.0

CCL_BEGIN_DECLS
/**
 * P(z) function.
 * This is a P(z) function (which can be user defined)
 * with a void* field to contain the parameters to that function.
 */
typedef struct {
        double (* your_pz_func)(double, double, void *, int*); /*< Function returns the likelihood of measuring a z_ph
        * (first double) given a z_spec (second double), with a pointer to additonal arguments and a status flag.*/
        void *  your_pz_params; /*< Additional parameters to be passed into your_pz_func */
} pz_info;

/**
 * dNdz function.
 * This is a dNdz function (which can be user defined)
 * with a void* field to contain the parameters to that function.
 */
typedef struct {
        double (* your_dN_func)(double, void *, int*); /*< Function returns the differential number density of galaxies wrt redshifts,
        * with a pointer to additonal arguments and a status flag.*/
        void *  your_dN_params; /*< Additional parameters to be passed into your_dN_func */
} dNdz_info;

/**
 * dNdz smail parmas structure.
 * This is a convenience parameters structure
 * to hold the three parameters of the Smail et al. analytic dNdz
 */
 typedef struct{
		double alpha;
		double beta;
		double z0;
	} smail_params;

/**
 * Return dNdz in a particular tomographic bin,
    convolved with a photo-z model (defined by the user), and normalized.
 * @param z redshift
 * @param dNdz_type the choice of dN/dz from Chang+
 * @param bin_zmin the minimum redshift of the tomorgraphic bin
 * @param bin_zmax the maximum redshift of the tomographic bin
 * @param photo_info the P(z) info struct
 * @param tomoout the output dN/dz
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * @return void
 */
void ccl_dNdz_tomog(double z, double bin_zmin, double bin_zmax, pz_info * photo_info,  dNdz_info * dN_info, double *tomoout, int *status);

/**
 * This function creates a structure amalgamating the information on an analytic true dNdz, plus some parameters.
 * @param params parameters for the analytic dNdz form
 * @param dNdz_func dNdz function
 * @return a structure with the dNdz and parameters
 */

dNdz_info* ccl_create_dNdz_info(void * params, double(*dNdz_func)(double,void*,int*));

/**
 * This function creates a structure containing the true dNdz for the built-in Smail-type analytic form:
 * dNdz ~ z^alpha exp(- (z/z0)^beta)
 * @param alpha
 * @param z0
 * @param beta
 * @return a structure with the built-in Smail-type dNdz and parameters
 */
 dNdz_info* ccl_create_Smail_dNdz_info(double alpha, double beta, double z0);

/** Free memory holding the structure containing dNdz information.
 * @param dN_info that holds user-defined dNdz and parameters
 * @return void
 */
void ccl_free_dNdz_info(dNdz_info * dN_info);

/**
 * This function creates a structure amalgamating the information on the photo-z model, P(z) plus some parameters.
 * @param params parameters for the P(z) function
 * @param pz_func P(z) function
 * @return a structure with the P(z) and parameters
 */

pz_info* ccl_create_photoz_info(void * params, double(*pz_func)(double, double,void*,int*));

/**
 * This function creates a structure containing the photo-z model for the built-in Gaussian photo-z pdf.
 * @param sigma_z0 The photo-z uncertainty at z=0. The photo-z uncertainty is assumed to scale like (1 + z).
 * @return a structure with the built-in Gaussian P(z) and parameters
 */
pz_info* ccl_create_gaussian_photoz_info(double sigma_z0);

/** Free memory holding the structure containing user-input photoz information.
 * @param my_photoz_info that holds user-defined P(z) and parameters
 * @return void
 */
void ccl_free_photoz_info(pz_info *my_photoz_info);

CCL_END_DECLS

#endif

ext: h

max_stars_repo_name: vrastil/CCL

max_issues_repo_name: vrastil/CCL

max_forks_repo_name: vrastil/CCL


text:
#ifndef NEURAL_MATRIX_H_
#define NEURAL_MATRIX_H_

#include <vector>
#include <functional>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <cassert>

#include <gsl/gsl_blas.h>

namespace neural {
    class Matrix {
        double *_data;
        int _rows;
        int _cols;

    public:
        Matrix(int rows, int cols);
        Matrix(int rows, int cols, double val);
        Matrix(int rows, int cols, double data[]);
        Matrix(int rows, int cols, std::vector<double> data);
        Matrix(std::vector<std::vector<double>> data);

        /**
         * Copy and move semantics
         */
        Matrix(const Matrix &rhs);
        Matrix(Matrix &&rhs);
        Matrix &operator=(const Matrix &rhs);
        Matrix &operator=(Matrix &&rhs);

        ~Matrix();

        /**
         * Addition operator
         */
        const Matrix operator+(const Matrix &rhs) const;

        /**
         * Subtraction operator
         */
        const Matrix operator-(const Matrix &rhs) const;

        /**
         * Matrix multiplication
         */
        const Matrix operator*(const Matrix &rhs) const;

        /**
         * Scale operator
         */
        const Matrix operator*(const double factor) const;

        /**
         * Division operator
         */
        const Matrix operator/(const double factor) const;

        /**
         * Hadamard operator
         */
        const Matrix operator^(const Matrix &rhs) const;

        /**
         * In-place addition operator
         */
        void operator+=(const Matrix &rhs);

        /**
         * In-place subtraction operator
         */
        void operator-=(const Matrix &rhs);

        /**
         * In-place matrix multiplication
         */
        void operator*=(const Matrix &rhs);

        /**
         * In-place scale operator
         */
        void operator*=(const double factor);

        /**
         * In-place division operator
         */
        void operator/=(const double factor);

        /**
         * In-place hadamard operator
         */
        void operator^=(const Matrix &rhs);

        /**
         * Call a function on each item of the matrix
         */
        void map(std::function<double(double&)> function);

        /**
         * Get the transpose of the matrix
         */
        Matrix transpose() const;

        /**
         * Calculate the Euclidean norm of the matrix
         */
        double norm();

        /**
         * Make all values zero
         */
        void zero();

        /**
         * Fetch a value from the matrix
         */
        inline double get_at(const int row, const int col) const {
            return _data[row * _cols + col];
        }

        /**
         * Set a value in the matrix
         */
        inline void set_at(const int row, const int col, const double value) {
            _data[row * _cols + col] = value;
        }

        /**
         * Get the number of rows in the matrix;
         */
        inline int get_rows() const {
            return _rows;
        }

        /**
         * Get the number of columns in the matrix
         */
        inline int get_cols() const {
            return _cols;
        }

        /**
         * Get the pointer to the matrix data
         */
        inline double *data() const {
            return _data;
        }

        /**
         * Pretty print the matrix for debugging
         */
        void print() const;
    };
}

#endif


text:
/* matrix/gsl_matrix_complex_long_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__
#define __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_complex_long_double.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  long double * data;
  gsl_block_complex_long_double * block;
  int owner;
} gsl_matrix_complex_long_double ;

typedef struct
{
  gsl_matrix_complex_long_double matrix;
} _gsl_matrix_complex_long_double_view;

typedef _gsl_matrix_complex_long_double_view gsl_matrix_complex_long_double_view;

typedef struct
{
  gsl_matrix_complex_long_double matrix;
} _gsl_matrix_complex_long_double_const_view;

typedef const _gsl_matrix_complex_long_double_const_view gsl_matrix_complex_long_double_const_view;

/* Allocation */

gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc (const size_t n1, const size_t n2);

gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_calloc (const size_t n1, const size_t n2);

gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc_from_block (gsl_block_complex_long_double * b,
                                           const size_t offset,
                                           const size_t n1, const size_t n2, const size_t d2);

gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc_from_matrix (gsl_matrix_complex_long_double * b,
                                            const size_t k1, const size_t k2,
                                            const size_t n1, const size_t n2);

gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_row_from_matrix (gsl_matrix_complex_long_double * m,
                                                const size_t i);

gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_col_from_matrix (gsl_matrix_complex_long_double * m,
                                                const size_t j);

void gsl_matrix_complex_long_double_free (gsl_matrix_complex_long_double * m);

/* Views */

_gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_submatrix (gsl_matrix_complex_long_double * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_row (gsl_matrix_complex_long_double * m, const size_t i);

_gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_column (gsl_matrix_complex_long_double * m, const size_t j);

_gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_diagonal (gsl_matrix_complex_long_double * m);

_gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subdiagonal (gsl_matrix_complex_long_double * m, const size_t k);

_gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_superdiagonal (gsl_matrix_complex_long_double * m, const size_t k);

_gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subrow (gsl_matrix_complex_long_double * m,
                                 const size_t i, const size_t offset,
                                 const size_t n);

_gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subcolumn (gsl_matrix_complex_long_double * m,
                                    const size_t j, const size_t offset,
                                    const size_t n);

_gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_array (long double * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_array_with_tda (long double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_vector (gsl_vector_complex_long_double * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_vector_with_tda (gsl_vector_complex_long_double * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_submatrix (const gsl_matrix_complex_long_double * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_row (const gsl_matrix_complex_long_double * m,
                            const size_t i);

_gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_column (const gsl_matrix_complex_long_double * m,
                               const size_t j);

_gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_diagonal (const gsl_matrix_complex_long_double * m);

_gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subdiagonal (const gsl_matrix_complex_long_double * m,
                                    const size_t k);

_gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_superdiagonal (const gsl_matrix_complex_long_double * m,
                                      const size_t k);

_gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subrow (const gsl_matrix_complex_long_double * m,
                                       const size_t i, const size_t offset,
                                       const size_t n);

_gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subcolumn (const gsl_matrix_complex_long_double * m,
                                          const size_t j, const size_t offset,
                                          const size_t n);

_gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_array (const long double * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_array_with_tda (const long double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_vector (const gsl_vector_complex_long_double * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_vector_with_tda (const gsl_vector_complex_long_double * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

void gsl_matrix_complex_long_double_set_zero (gsl_matrix_complex_long_double * m);
void gsl_matrix_complex_long_double_set_identity (gsl_matrix_complex_long_double * m);
void gsl_matrix_complex_long_double_set_all (gsl_matrix_complex_long_double * m, gsl_complex_long_double x);

int gsl_matrix_complex_long_double_fread (FILE * stream, gsl_matrix_complex_long_double * m) ;
int gsl_matrix_complex_long_double_fwrite (FILE * stream, const gsl_matrix_complex_long_double * m) ;
int gsl_matrix_complex_long_double_fscanf (FILE * stream, gsl_matrix_complex_long_double * m);
int gsl_matrix_complex_long_double_fprintf (FILE * stream, const gsl_matrix_complex_long_double * m, const char * format);

int gsl_matrix_complex_long_double_memcpy(gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);
int gsl_matrix_complex_long_double_swap(gsl_matrix_complex_long_double * m1, gsl_matrix_complex_long_double * m2);

int gsl_matrix_complex_long_double_swap_rows(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
int gsl_matrix_complex_long_double_swap_columns(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
int gsl_matrix_complex_long_double_swap_rowcol(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);

int gsl_matrix_complex_long_double_transpose (gsl_matrix_complex_long_double * m);
int gsl_matrix_complex_long_double_transpose_memcpy (gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);

int gsl_matrix_complex_long_double_equal (const gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);

int gsl_matrix_complex_long_double_isnull (const gsl_matrix_complex_long_double * m);
int gsl_matrix_complex_long_double_ispos (const gsl_matrix_complex_long_double * m);
int gsl_matrix_complex_long_double_isneg (const gsl_matrix_complex_long_double * m);
int gsl_matrix_complex_long_double_isnonneg (const gsl_matrix_complex_long_double * m);

int gsl_matrix_complex_long_double_add (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
int gsl_matrix_complex_long_double_sub (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
int gsl_matrix_complex_long_double_mul_elements (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
int gsl_matrix_complex_long_double_div_elements (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
int gsl_matrix_complex_long_double_scale (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);
int gsl_matrix_complex_long_double_add_constant (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);
int gsl_matrix_complex_long_double_add_diagonal (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_complex_long_double_get_row(gsl_vector_complex_long_double * v, const gsl_matrix_complex_long_double * m, const size_t i);
int gsl_matrix_complex_long_double_get_col(gsl_vector_complex_long_double * v, const gsl_matrix_complex_long_double * m, const size_t j);
int gsl_matrix_complex_long_double_set_row(gsl_matrix_complex_long_double * m, const size_t i, const gsl_vector_complex_long_double * v);
int gsl_matrix_complex_long_double_set_col(gsl_matrix_complex_long_double * m, const size_t j, const gsl_vector_complex_long_double * v);
/***********************************************************************/

/* inline functions if you are using GCC */

INLINE_DECL gsl_complex_long_double gsl_matrix_complex_long_double_get(const gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
INLINE_DECL void gsl_matrix_complex_long_double_set(gsl_matrix_complex_long_double * m, const size_t i, const size_t j, const gsl_complex_long_double x);

INLINE_DECL gsl_complex_long_double * gsl_matrix_complex_long_double_ptr(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
INLINE_DECL const gsl_complex_long_double * gsl_matrix_complex_long_double_const_ptr(const gsl_matrix_complex_long_double * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex_long_double
gsl_matrix_complex_long_double_get(const gsl_matrix_complex_long_double * m,
                     const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      gsl_complex_long_double zero = {{0,0}};

      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, zero) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, zero) ;
        }
    }
#endif
  return *(gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
void
gsl_matrix_complex_long_double_set(gsl_matrix_complex_long_double * m,
                     const size_t i, const size_t j, const gsl_complex_long_double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  *(gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) = x ;
}

INLINE_FUN
gsl_complex_long_double *
gsl_matrix_complex_long_double_ptr(gsl_matrix_complex_long_double * m,
                             const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
const gsl_complex_long_double *
gsl_matrix_complex_long_double_const_ptr(const gsl_matrix_complex_long_double * m,
                                   const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game