text:
/**
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_interp.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "inout_aper.h"
#include "aper_conf.h"
#include "spc_spc.h"
#include "spce_PET.h"
#include "fringe_conf.h"
#include "ipixcorr_utils.h"

#define AXE_IMAGE_PATH "AXE_IMAGE_PATH"
#define AXE_OUTPUT_PATH "AXE_OUTPUT_PATH"
#define AXE_CONFIG_PATH "AXE_CONFIG_PATH"

int
main (int argc, char *argv[])
{
  char *opt;
  char aper_file[MAXCHAR];
  char aper_file_path[MAXCHAR];

  char conf_file[MAXCHAR];
  char conf_file_path[MAXCHAR];

  char grism_file[MAXCHAR];
  char grism_file_path[MAXCHAR];

  char IPC_file[MAXCHAR];
  char IPC_file_path[MAXCHAR];

  char PET_file[MAXCHAR];
  char PET_file_path[MAXCHAR];

  char IPC_func[MAXCHAR];
  char IPC_func_path[MAXCHAR];
  char ID[60];

  aperture_conf *conf;

  object **oblist;
  observation *obs = NULL;

  interpolator *ipcorr;

  int index;
  //int i;

  ap_pixel *PET;

  fitsfile  *PET_ptr;
  fitsfile  *IPC_ptr;
  FITScards *cards;

  int f_status=0;
  int bckmode = 0;

  int aperID=0, beamID=0;
  //int objindex;

  int spec_OAF=0;
  int point_like=0;
  int origname=0;

  beam act_beam;

  double max_ext;

  if ((argc < 3) || (opt = get_online_option ("help", argc, argv)))
    {
      fprintf (stdout,
         "aXe_FRINGECORR Version %s:\n",RELEASE);
      exit (1);
    }

  fprintf (stdout, "aXe_PETIPC: Starting...\n");

  // reset the parameter index
  index = 0;

  // read in the grism image name
  strcpy(grism_file, argv[++index]);
  build_path (AXE_IMAGE_PATH, grism_file, grism_file_path);

  // read in the configuration file name
  strcpy(conf_file, argv[++index]);
  build_path (AXE_CONFIG_PATH, conf_file, conf_file_path);

  // load the configuration file
  conf = get_aperture_descriptor(conf_file_path);

  // determine the science extension number
  get_extension_numbers(grism_file_path, conf,conf->optkey1,conf->optval1);

  // Determine if we are using the special bck mode
  // In this mode, file names are handled diferently
  if ((opt = get_online_option("bck", argc, argv)))
    bckmode = 1;

  // determine whether the corrected PET should maintain the
  // filename of the original PET
  if ((opt = get_online_option("origname", argc, argv)))
    origname = 1;

  // get or set up the input OAF name
  if ((opt = get_online_option("in_OAF", argc, argv)))
    {
      // copy the parameter value if given
      strcpy(aper_file, opt);

      // mark that there is a special OAF file
      spec_OAF=1;
    }
  else
    {
      // copy the grism name and construct the name
      // if not given as parameter
      strcpy(aper_file, grism_file);
      replace_file_extension (grism_file, aper_file, ".fits",
			      ".OAF", conf->science_numext);
    }
  build_path(AXE_OUTPUT_PATH, aper_file, aper_file_path);

  // get or set up the output PET name
  if ((opt = get_online_option("out_PET", argc, argv)))
    {
      // copy the parameter value if given
      strcpy(IPC_file, opt);
    }
  else
    {
      // copy the grism name and construct the name
      // if not given as parameter
      strcpy(IPC_file, grism_file);

      if (bckmode)
	replace_file_extension(grism_file, IPC_file, ".fits",
			       "_ipc.BCK.PET.fits", conf->science_numext);
      else
	replace_file_extension(grism_file, IPC_file, ".fits",
			       "_ipc.PET.fits", conf->science_numext);
    }
  build_path(AXE_OUTPUT_PATH, IPC_file, IPC_file_path);

  // get or set up the output PET name
  if ((opt = get_online_option("in_PET", argc, argv)))
    {
      // copy the parameter value if given
      strcpy(PET_file, opt);
    }
  else
    {
      // copy the grism name and construct the name
      // if not given as parameter
      strcpy(PET_file, grism_file);
      if (bckmode)
	replace_file_extension(grism_file, PET_file, ".fits",
			       ".BCK.PET.fits", conf->science_numext);
      else
	replace_file_extension(grism_file, PET_file, ".fits",
			       ".PET.fits", conf->science_numext);

    }
  build_path(AXE_OUTPUT_PATH, PET_file, PET_file_path);

  // Build the intrapixel correction file name
  sprintf(IPC_func,"%s",conf->IPIXfunc);

  // overwrite it with the parameter if given
  if ((opt = get_online_option ("IPIXFUNCTION", argc, argv)))
    strcpy(IPC_func,opt);

  if (!strcmp(IPC_file,"None"))
    // check whether the correction function is given;
    // complain if not given!
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "aXe_PETIPC: Correction function must be given!!\n");
  else
    // compose the full pathname
    build_path (AXE_CONFIG_PATH, IPC_func, IPC_func_path);

  // get the maximum extension number
  if ((opt = get_online_option ("max_ext", argc, argv)))
    max_ext = atof(opt);
  else
    max_ext = 0.0;

  // check whether either max_ext is defined
  // OR a special OAF file is given
  if (!max_ext && !spec_OAF)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "aXe_PETIPC: Either maximum extension\n"
		 "OR a special OAF file must be given!\n");

  fprintf (stdout, "aXe_PETIPC: Grism file name:                     %s\n",
	   grism_file_path);
  fprintf (stdout, "aXe_PETIPC: Configuration file name:             %s\n",
	   conf_file_path);
  fprintf (stdout, "aXe_PETIPC: Correction function:                 %s\n",
	   IPC_func_path);
  fprintf (stdout, "aXe_PETIPC: Input OAF file name:                 %s\n",
	   aper_file_path);
  fprintf (stdout, "aXe_PETIPC: Input PET file name:                 %s\n",
	   PET_file_path);
  if (origname)
    fprintf (stdout, "aXe_PETIPC: Output PET file name:                %s\n",
	     PET_file_path);
  else
    fprintf (stdout, "aXe_PETIPC: Output PET file name:                %s\n",
	     IPC_file_path);
  fprintf (stdout, "aXe_PETIPC: Maximal extension for point objects: %f\n",
	   max_ext);

  // load the aperture list
  fprintf (stdout, "aXe_PETIPC: Loading object aperture list...");
  fflush(stdout);
  oblist = file_to_object_list_seq (aper_file_path, obs);
  fprintf (stdout,"%d objects loaded.\n",object_list_size(oblist));
  fflush(stdout);

  // Open the PET file for reading
  fits_open_file (&PET_ptr, PET_file_path, READONLY, &f_status);
  if (f_status)
    {
      ffrprt (stdout, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_PETIPC: Could not open file: %s\n",
		   PET_file_path);
    }

  // open the output PET file
  IPC_ptr = create_PET_opened(IPC_file_path,1);

  // copy the keywords from the input
  // to the output PET
  cards = get_FITS_cards_opened(PET_ptr);
  put_FITS_cards_opened(IPC_ptr, cards);
  free_FITScards(cards);

  // create an interpolator reading in the fits file
  ipcorr = create_interp_ftable(IPC_func_path,2,"X","FX",IPCORR_INTERP_TYPE);

  while ((aperID!=-1) || (beamID!=-1))
    {
      // Get the PET for this object
      PET = get_ALL_from_next_in_PET(PET_ptr, &aperID, &beamID);
      if ((aperID==-1) && (beamID==-1))
	break;

      // PET is empty, skip it
      if (PET==NULL)
	{
	  fprintf(stdout, "Empty PET BEAM %d%c...\n", aperID, BEAM(beamID));
	  continue; }

      // select the actual beam from the beam list
      act_beam = find_beam_in_object_list(oblist, aperID, beamID);

      // determine whether the object is pointlike
      point_like = is_pointlike(act_beam, spec_OAF, max_ext);

      if (point_like)
	{
	  fprintf(stdout, "aXe_PETIPC: Correcting BEAM %d%c...\n", aperID, BEAM(beamID));
	  intpix_corr_pet(act_beam, conf_file_path, ipcorr, PET);
	}
      else
	{
	  fprintf(stdout, "aXe_PETIPC: Skipping BEAM %d%c...\n", aperID, BEAM(beamID));
	}

      sprintf(ID, "%d%c", aperID, BEAM (beamID));
      add_ALL_to_PET(PET, ID, IPC_ptr, 0);

      /* Copy header from OPET extension into this SPC extension */
      cards = get_FITS_cards_opened(PET_ptr);
      put_FITS_cards_opened(IPC_ptr,cards);
      free_FITScards(cards);

      if (PET!=NULL)
	free(PET);
    }

  // free the object list
  if (oblist!=NULL)
    free_oblist (oblist);

  // free the interpolator
  free_interp(ipcorr);

  // free the configuration structure
  free_aperture_conf(conf);

  fits_close_file (PET_ptr, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_PETIPC: " "Error closing PET: %s\n",
		   PET_file_path);
    }

  fits_close_file (IPC_ptr, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_PETIPC: " "Error closing PET: %s\n",
		   IPC_file_path);
    }

  // check whether the old
  // PET name should be retained
  if (origname)
    {
      // delete the original
      // PET file
      if (remove(PET_file_path) != 0)
	fprintf (stdout, "aXe_PETIPC: Problems deleting %s!\n", PET_file_path);

      // rename the corrcted PET file
      if (rename(IPC_file_path, PET_file_path) !=0)
	fprintf (stdout, "aXe_PETIPC: Problems renaming %s to %s!\n", IPC_file_path, PET_file_path);

    }

  // say goodbye
  fprintf (stdout, "aXe_PETIPC: Done...\n");
  exit (0);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/*
Copyright 2017 Jiawei Chiu

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
#pragma once

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <memory>
#include <random>
#include <sstream>
#include <string>
#include <vector>

#include <gflags/gflags.h>
#include <glog/logging.h>

#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <curand.h>
#include <cusolverDn.h>
#include <cusparse.h>

// When BLAS functions are missing, we fallback on Eigen.
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/SparseCore>

#include <cublas_v2.h>
// #include <magma_lapack.h>
#include <magma_v2.h>

#include <lapacke.h>

namespace gi {

using std::abs;
using std::cout;
using std::ifstream;
using std::istream;
using std::istringstream;
using std::max;
using std::min;
using std::ostream;
using std::ofstream;
using std::ostringstream;
using std::sqrt;
using std::string;
using std::unique_ptr;
using std::vector;

void MainInit(int argc, char **argv);

class Timer {
public:
  Timer();
  void reset();
  double elapsed() const;

private:
  typedef std::chrono::high_resolution_clock clock_;
  typedef std::chrono::duration<double, std::ratio<1> > second_;
  std::chrono::time_point<clock_> beg_;
};

struct EngineOptions {
  int omp_num_threads = 1;
  int device = 0;      // For GPU device.
  int num_streams = 4; // Number of CUDA streams.
  int rng_seed = 56150941;
};

#define CUDA_CALL(x)                                                           \
  { CHECK_EQ(x, cudaSuccess) << "CUDA error: " << x; }

#define CURAND_CALL(x)                                                         \
  { CHECK_EQ(x, CURAND_STATUS_SUCCESS) << "CURAND error: " << x; }

#define CUBLAS_CALL(x)                                                         \
  { CHECK_EQ(x, CUBLAS_STATUS_SUCCESS) << "CUBLAS error: " << x; }

#define CUSOLVER_CALL(x)                                                       \
  { CHECK_EQ(x, CUSOLVER_STATUS_SUCCESS) << "CUSOLVER error: " << x; }

#define CUSPARSE_CALL(x)                                                       \
  { CHECK_EQ(x, CUSPARSE_STATUS_SUCCESS) << "CUSPARSE error: " << x; }

class Engine {
public:
  Engine(const EngineOptions &opt);
  ~Engine();

  static Engine *instance() { return instance_; }
  static cudaStream_t stream(int i) { return instance()->stream_[i]; }
  static curandGenerator_t curand() { return instance()->curand_; }
  static cublasHandle_t cublas() { return instance()->cublas_; }
  static cusparseHandle_t cusparse() { return instance()->cusparse_; }
  static cusparseMatDescr_t cusparse_desc() {
    return instance()->cusparse_desc_;
  }
  static cusolverDnHandle_t cusolver_dn() { return instance()->cusolver_dn_; }
  static std::mt19937 &rng() { return instance()->rng_; }
  static float *device_s_one() { return instance()->device_s_one_; }
  static float *device_s_zero() { return instance()->device_s_zero_; }

private:
  vector<cudaStream_t> stream_;
  curandGenerator_t curand_;
  cublasHandle_t cublas_;
  cusparseHandle_t cusparse_;
  cusparseMatDescr_t cusparse_desc_;
  cusolverDnHandle_t cusolver_dn_;
  std::mt19937 rng_;
  float *device_s_one_;
  float *device_s_zero_;

  static Engine *instance_;
};

// Some convenience functions for GPU kernels.

// Returns number of blocks and number of threads given n things to work on.
void BlocksThreads(int max_blocks, int max_threads, int n, int *num_blocks,
                   int *num_threads);
}

ext: h

max_stars_repo_name: tinkerstash/gpuimpute

max_issues_repo_name: tinkerstash/gpuimpute

max_forks_repo_name: tinkerstash/gpuimpute


text:
/* ode-initval/test.c
 *
 * Copyright (C) 2009, 2010 Tuomo Keskitalo
 *
 * 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.
 */

/* Some functions and tests based on test.c by G. Jungman. */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_odeiv2.h>
#include "odeiv_util.h"

/* Maximum number of ODE equations */
#define MAXEQ 15

/* Maximum number of ODE solvers */
#define MAXNS 20

/* Track number of function and jacobian evaluations in tests
   with global variables
 */
int nfe;
int nje;

/**********************************************************/
/* ODE test system definitions                            */
/**********************************************************/

/* RHS for f=2. Solution y = 2 * t + t0 */

int
rhs_linear (double t, const double y[], double f[], void *params)
{
  extern int nfe;
  nfe += 1;

  f[0] = 2.0;

  return GSL_SUCCESS;
}

int
jac_linear (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  extern int nje;
  nje += 1;

  dfdy[0] = 0.0;
  dfdt[0] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_lin = {
  rhs_linear,
  jac_linear,
  1,
  0
};

/* RHS for f=y. Equals y=exp(t) with initial value y(0)=1.0 */

int
rhs_exp (double t, const double y[], double f[], void *params)
{
  extern int nfe;
  nfe += 1;

  f[0] = y[0];

  return GSL_SUCCESS;
}

int
jac_exp (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  extern int nje;
  nje += 1;

  dfdy[0] = y[0];
  dfdt[0] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_exp = {
  rhs_exp,
  jac_exp,
  1,
  0
};

int
rhs_sin (double t, const double y[], double f[], void *params)
{
  extern int nfe;
  nfe += 1;

  f[0] = -y[1];
  f[1] = y[0];

  return GSL_SUCCESS;
}

int
jac_sin (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  extern int nje;
  nje += 1;

  dfdy[0] = 0.0;
  dfdy[1] = -1.0;
  dfdy[2] = 1.0;
  dfdy[3] = 0.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_sin = {
  rhs_sin,
  jac_sin,
  2,
  0
};

/*  Sine/cosine with random failures */

static int rhs_xsin_reset = 0;
static int jac_xsin_reset = 0;

int
rhs_xsin (double t, const double y[], double f[], void *params)
{
  static int n = 0, m = 0;
  extern int nfe;
  nfe += 1;

  if (rhs_xsin_reset)
    {
      rhs_xsin_reset = 0;
      n = 0;
      m = 1;
    }
  n++;

  if (n >= m)
    {
      m = n * 1.3;
      return GSL_EFAILED;
    }

  if (n > 40 && n < 65)
    {
      f[0] = GSL_NAN;
      f[1] = GSL_NAN;
      return GSL_EFAILED;
    }

  f[0] = -y[1];
  f[1] = y[0];

  return GSL_SUCCESS;
}

int
jac_xsin (double t, const double y[], double *dfdy, double dfdt[],
          void *params)
{
  static int n = 0;
  extern int nje;
  nje += 1;

  if (jac_xsin_reset)
    {
      jac_xsin_reset = 0;
      n = 0;
    }
  n++;

  if (n > 50 && n < 55)
    {
      dfdy[0] = GSL_NAN;
      dfdy[1] = GSL_NAN;
      dfdy[2] = GSL_NAN;
      dfdy[3] = GSL_NAN;

      dfdt[0] = GSL_NAN;
      dfdt[1] = GSL_NAN;
      return GSL_EFAILED;
    }

  dfdy[0] = 0.0;
  dfdy[1] = -1.0;
  dfdy[2] = 1.0;
  dfdy[3] = 0.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_xsin = {
  rhs_xsin,
  jac_xsin,
  2,
  0
};

/* RHS for classic stiff example
   dy0 / dt =  998 * y0 + 1998 * y1    y0(0) = 1.0
   dy1 / dt = -999 * y0 - 1999 * y1    y1(0) = 0.0

   solution is
   y0 = 2 * exp(-t) - exp(-1000 * t)
   y1 = - exp(-t) + exp(-1000 * t)
*/

int
rhs_stiff (double t, const double y[], double f[], void *params)
{
  extern int nfe;
  nfe += 1;

  f[0] = 998.0 * y[0] + 1998.0 * y[1];
  f[1] = -999.0 * y[0] - 1999.0 * y[1];

  return GSL_SUCCESS;
}

int
jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
           void *params)
{
  extern int nje;
  nje += 1;

  dfdy[0] = 998.0;
  dfdy[1] = 1998.0;
  dfdy[2] = -999.0;
  dfdy[3] = -1999.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_stiff = {
  rhs_stiff,
  jac_stiff,
  2,
  0
};

/* Cosine function */

int
rhs_cos (double t, const double *y, double *dydt, void *params)
{
  dydt[0] = cos (t);
  return GSL_SUCCESS;
}

int
jac_cos (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  dfdy[0] = 0.0;
  dfdt[0] = -sin (t);

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_cos = {
  rhs_cos,
  jac_cos,
  1,
  0
};

/* Broken problem for testing numerical problems in user function that
   leads to decrease of step size in gsl_odeiv2_evolve below machine
   precision.
 */

int
rhs_broken (double t, const double y[], double f[], void *params)
{
  if (t < 10.0)
    {
      f[0] = 1.0;
    }

  else
    {
      f[0] = GSL_NAN;
      return 123;
    }

  return GSL_SUCCESS;
}

int
jac_broken (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  if (t < 10.0)
    {
      dfdy[0] = 0.0;
      dfdt[0] = 0.0;
    }
  else
    {
      dfdy[0] = GSL_NAN;
      dfdt[0] = GSL_NAN;
      return 123;
    }

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_broken = {
  rhs_broken,
  jac_broken,
  1,
  0
};

/* Immediate user break (at t > 1.5) test sine system */

int
rhs_sin_ub (double t, const double y[], double f[], void *params)
{
  extern int nfe;
  nfe += 1;

  f[0] = -y[1];
  f[1] = y[0];

  if (t > 1.5)
    {
      return GSL_EBADFUNC;
    }

  return GSL_SUCCESS;
}

int
jac_sin_ub (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  extern int nje;
  nje += 1;

  dfdy[0] = 0.0;
  dfdy[1] = -1.0;
  dfdy[2] = 1.0;
  dfdy[3] = 0.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  if (t > 1.5)
    {
      return GSL_EBADFUNC;
    }

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_sin_ub = {
  rhs_sin_ub,
  jac_sin_ub,
  2,
  0
};

/* Badly scaled random function */

int
rhs_br (double t, const double *y, double *dydt, void *params)
{
  dydt[0] = (rand () - RAND_MAX / 2) * 2e100;
  return GSL_SUCCESS;
}

int
jac_br (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
  dfdy[0] = (rand () - RAND_MAX / 2) * 2e100;
  dfdt[0] = (rand () - RAND_MAX / 2) * 2e100;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_br = {
  rhs_br,
  jac_br,
  1,
  0
};

/* stepfn and stepfn2 based on testcases from Frank Reininghaus
   <frank78ac@googlemail.com>
*/

/* Derivative change at t=0, small derivative */

int
rhs_stepfn (double t, const double *y, double *dydt, void *params)
{
  if (t >= 1.0)
    dydt[0] = 1;
  else
    dydt[0] = 0;

  return GSL_SUCCESS;
}

int
jac_stepfn (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  dfdy[0] = 0.0;
  dfdt[0] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_stepfn = {
  rhs_stepfn,
  jac_stepfn,
  1,
  0
};

/* Derivative change at t=0, large derivative */

int
rhs_stepfn2 (double t, const double *y, double *dydt, void *params)
{
  if (t >= 0.0)
    dydt[0] = 1e300;
  else
    dydt[0] = 0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_stepfn2 = {
  rhs_stepfn2,
  jac_stepfn,
  1,
  0
};

/* Volterra-Lotka predator-prey model
   f0 = (a - b * y1) * y0     y0(0) = 2.725
   f1 = (-c + d * y0) * y1    y1(0) = 1.0
 */

int
rhs_vl (double t, const double y[], double f[], void *params)
{
  const double a = -1.0;
  const double b = -1.0;
  const double c = -2.0;
  const double d = -1.0;

  extern int nfe;
  nfe += 1;

  f[0] = (a - b * y[1]) * y[0];
  f[1] = (-c + d * y[0]) * y[1];

  return GSL_SUCCESS;
}

int
jac_vl (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
  const double a = -1.0;
  const double b = -1.0;
  const double c = -2.0;
  const double d = -1.0;

  extern int nje;
  nje += 1;

  dfdy[0] = a - b * y[1];
  dfdy[1] = -b * y[0];
  dfdy[2] = d * y[1];
  dfdy[3] = -c + d * y[0];

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_vl = {
  rhs_vl,
  jac_vl,
  2,
  0
};

/* van Der Pol oscillator
   f0 = y1                           y0(0) = 1.0
   f1 = -y0 + mu * y1 * (1 - y0^2)   y1(0) = 0.0
*/

int
rhs_vanderpol (double t, const double y[], double f[], void *params)
{
  const double mu = 10.0;

  extern int nfe;
  nfe += 1;

  f[0] = y[1];
  f[1] = -y[0] + mu * y[1] * (1.0 - y[0] * y[0]);

  return GSL_SUCCESS;
}

int
jac_vanderpol (double t, const double y[], double *dfdy, double dfdt[],
               void *params)
{
  const double mu = 10.0;

  extern int nje;
  nje += 1;

  dfdy[0] = 0.0;
  dfdy[1] = 1.0;
  dfdy[2] = -2.0 * mu * y[0] * y[1] - 1.0;
  dfdy[3] = mu * (1.0 - y[0] * y[0]);

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_vanderpol = {
  rhs_vanderpol,
  jac_vanderpol,
  2,
  0
};

/* Stiff trigonometric example
   f0 = -50 * (y0 - cos(t))    y0(0) = 0.0
 */

int
rhs_stifftrig (double t, const double y[], double f[], void *params)
{
  extern int nfe;
  nfe += 1;

  f[0] = -50 * (y[0] - cos (t));

  return GSL_SUCCESS;
}

int
jac_stifftrig (double t, const double y[], double *dfdy, double dfdt[],
               void *params)
{
  extern int nje;
  nje += 1;

  dfdy[0] = -50;

  dfdt[0] = -50 * sin (t);

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_stifftrig = {
  rhs_stifftrig,
  jac_stifftrig,
  1,
  0
};

/* The Oregonator - chemical Belusov-Zhabotinskii reaction
   y0(0) = 1.0, y1(0) = 2.0, y2(0) = 3.0
*/

int
rhs_oregonator (double t, const double y[], double f[], void *params)
{
  const double c1 = 77.27;
  const double c2 = 8.375e-6;
  const double c3 = 0.161;

  extern int nfe;
  nfe += 1;

  f[0] = c1 * (y[1] + y[0] * (1 - c2 * y[0] - y[1]));
  f[1] = 1 / c1 * (y[2] - y[1] * (1 + y[0]));
  f[2] = c3 * (y[0] - y[2]);

  return GSL_SUCCESS;
}

int
jac_oregonator (double t, const double y[], double *dfdy, double dfdt[],
                void *params)
{
  const double c1 = 77.27;
  const double c2 = 8.375e-6;
  const double c3 = 0.161;

  extern int nje;
  nje += 1;

  dfdy[0] = c1 * (1 - 2 * c2 * y[0] - y[1]);
  dfdy[1] = c1 * (1 - y[0]);
  dfdy[2] = 0.0;

  dfdy[3] = 1 / c1 * (-y[1]);
  dfdy[4] = 1 / c1 * (-1 - y[0]);
  dfdy[5] = 1 / c1;

  dfdy[6] = c3;
  dfdy[7] = 0.0;
  dfdy[8] = -c3;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_oregonator = {
  rhs_oregonator,
  jac_oregonator,
  3,
  0
};

/* E5 - a stiff badly scaled chemical problem by Enright, Hull &
   Lindberg (1975): Comparing numerical methods for stiff systems of
   ODEs. BIT, vol. 15, pp. 10-48.

   f0 = -a * y0 - b * y0 * y2                            y0(0) = 1.76e-3
   f1 = a * y0 - m * c * y1 * y2                         y1(0) = 0.0
   f2 = a * y0 - b * y0 * y2 - m * c * y1 * y2 + c * y3  y2(0) = 0.0
   f3 = b * y0 * y2 - c * y3                             y3(0) = 0.0
 */

int
rhs_e5 (double t, const double y[], double f[], void *params)
{
  const double a = 7.89e-10;
  const double b = 1.1e7;
  const double c = 1.13e3;
  const double m = 1.0e6;

  extern int nfe;
  nfe += 1;

  f[0] = -a * y[0] - b * y[0] * y[2];
  f[1] = a * y[0] - m * c * y[1] * y[2];
  f[3] = b * y[0] * y[2] - c * y[3];
  f[2] = f[1] - f[3];

  return GSL_SUCCESS;
}

int
jac_e5 (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
  const double a = 7.89e-10;
  const double b = 1.1e7;
  const double c = 1.13e3;
  const double m = 1.0e6;

  extern int nje;
  nje += 1;

  dfdy[0] = -a - b * y[2];
  dfdy[1] = 0.0;
  dfdy[2] = -b * y[0];
  dfdy[3] = 0.0;

  dfdy[4] = a;
  dfdy[5] = -m * c * y[2];
  dfdy[6] = -m * c * y[1];
  dfdy[7] = 0.0;

  dfdy[8] = a - b * y[2];
  dfdy[9] = -m * c * y[2];
  dfdy[10] = -b * y[0] - m * c * y[1];
  dfdy[11] = c;

  dfdy[12] = b * y[2];
  dfdy[13] = 0.0;
  dfdy[14] = b * y[0];
  dfdy[15] = -c;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;
  dfdt[3] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_e5 = {
  rhs_e5,
  jac_e5,
  4,
  0
};

/* Chemical reaction system of H.H. Robertson (1966): The solution of
   a set of reaction rate equations. In: J. Walsh, ed.: Numer.
   Anal., an Introduction, Academ. Press, pp. 178-182.

   f0 = -a * y0 + b * y1 * y2             y0(0) = 1.0
   f1 = a * y0 - b * y1 * y2 - c * y1^2   y1(0) = 0.0
   f2 = c * y1^2                          y2(0) = 0.0
 */

int
rhs_robertson (double t, const double y[], double f[], void *params)
{
  const double a = 0.04;
  const double b = 1.0e4;
  const double c = 3.0e7;

  extern int nfe;
  nfe += 1;

  f[0] = -a * y[0] + b * y[1] * y[2];
  f[2] = c * y[1] * y[1];
  f[1] = -f[0] - f[2];

  return GSL_SUCCESS;
}

int
jac_robertson (double t, const double y[], double *dfdy, double dfdt[],
               void *params)
{
  const double a = 0.04;
  const double b = 1.0e4;
  const double c = 3.0e7;

  extern int nje;
  nje += 1;

  dfdy[0] = -a;
  dfdy[1] = b * y[2];
  dfdy[2] = b * y[1];

  dfdy[3] = a;
  dfdy[4] = -b * y[2] - 2 * c * y[1];
  dfdy[5] = -b * y[1];

  dfdy[6] = 0.0;
  dfdy[7] = 2 * c * y[1];
  dfdy[8] = 0.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_robertson = {
  rhs_robertson,
  jac_robertson,
  3,
  0
};

/* A two-dimensional oscillating Brusselator system.

   f0 = a + y0^2 * y1 - (b + 1) * y0      y0(0) = 1.5
   f1 = b * y0 - y0^2 * y1                y1(0) = 3.0
 */

int
rhs_brusselator (double t, const double y[], double f[], void *params)
{
  const double a = 1.0;
  const double b = 3.0;

  extern int nfe;
  nfe += 1;

  f[0] = a + y[0] * y[0] * y[1] - (b + 1.0) * y[0];
  f[1] = b * y[0] - y[0] * y[0] * y[1];

  return GSL_SUCCESS;
}

int
jac_brusselator (double t, const double y[], double *dfdy, double dfdt[],
                 void *params)
{
  const double b = 3.0;

  extern int nje;
  nje += 1;

  dfdy[0] = 2 * y[0] * y[1] - (b + 1.0);
  dfdy[1] = y[0] * y[0];

  dfdy[2] = b - 2 * y[0] * y[1];
  dfdy[3] = -y[0] * y[0];

  dfdt[0] = 0;
  dfdt[1] = 0;

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_brusselator = {
  rhs_brusselator,
  jac_brusselator,
  2,
  0
};

/* Ring Modulator, stiff ODE of dimension 15.

   Reference: Walter M. Lioen, Jacques J.B. de Swart, Test Set for
   Initial Value Problem Solvers, Release 2.1 September 1999,
   http://ftp.cwi.nl/IVPtestset/software.htm
 */

#define NRINGMOD 15

int
rhs_ringmod (double t, const double y[], double f[], void *params)
{
  const double c = 1.6e-8;
  const double cs = 2e-12;
  const double cp = 1e-8;
  const double r = 25e3;
  const double rp = 50e0;
  const double lh = 4.45e0;
  const double ls1 = 2e-3;
  const double ls2 = 5e-4;
  const double ls3 = 5e-4;
  const double rg1 = 36.3;
  const double rg2 = 17.3;
  const double rg3 = 17.3;
  const double ri = 5e1;
  const double rc = 6e2;
  const double gamma = 40.67286402e-9;
  const double delta = 17.7493332;
  const double pi = 3.141592653589793238462643383;

  const double uin1 = 0.5 * sin (2e3 * pi * t);
  const double uin2 = 2 * sin (2e4 * pi * t);
  const double ud1 = +y[2] - y[4] - y[6] - uin2;
  const double ud2 = -y[3] + y[5] - y[6] - uin2;
  const double ud3 = +y[3] + y[4] + y[6] + uin2;
  const double ud4 = -y[2] - y[5] + y[6] + uin2;

  const double qud1 = gamma * (exp (delta * ud1) - 1.0);
  const double qud2 = gamma * (exp (delta * ud2) - 1.0);
  const double qud3 = gamma * (exp (delta * ud3) - 1.0);
  const double qud4 = gamma * (exp (delta * ud4) - 1.0);

  extern int nfe;
  nfe += 1;

  f[0] = (y[7] - 0.5 * y[9] + 0.5 * y[10] + y[13] - y[0] / r) / c;
  f[1] = (y[8] - 0.5 * y[11] + 0.5 * y[12] + y[14] - y[1] / r) / c;
  f[2] = (y[9] - qud1 + qud4) / cs;
  f[3] = (-y[10] + qud2 - qud3) / cs;
  f[4] = (y[11] + qud1 - qud3) / cs;
  f[5] = (-y[12] - qud2 + qud4) / cs;
  f[6] = (-y[6] / rp + qud1 + qud2 - qud3 - qud4) / cp;
  f[7] = -y[0] / lh;
  f[8] = -y[1] / lh;
  f[9] = (0.5 * y[0] - y[2] - rg2 * y[9]) / ls2;
  f[10] = (-0.5 * y[0] + y[3] - rg3 * y[10]) / ls3;
  f[11] = (0.5 * y[1] - y[4] - rg2 * y[11]) / ls2;
  f[12] = (-0.5 * y[1] + y[5] - rg3 * y[12]) / ls3;
  f[13] = (-y[0] + uin1 - (ri + rg1) * y[13]) / ls1;
  f[14] = (-y[1] - (rc + rg1) * y[14]) / ls1;

  return GSL_SUCCESS;
}

int
jac_ringmod (double t, const double y[], double *dfdy, double dfdt[],
             void *params)
{
  const double c = 1.6e-8;
  const double cs = 2e-12;
  const double cp = 1e-8;
  const double r = 25e3;
  const double rp = 50e0;
  const double lh = 4.45e0;
  const double ls1 = 2e-3;
  const double ls2 = 5e-4;
  const double ls3 = 5e-4;
  const double rg1 = 36.3;
  const double rg2 = 17.3;
  const double rg3 = 17.3;
  const double ri = 5e1;
  const double rc = 6e2;
  const double gamma = 40.67286402e-9;
  const double delta = 17.7493332;
  const double pi = 3.141592653589793238462643383;

  const double uin2 = 2 * sin (2e4 * pi * t);
  const double ud1 = +y[2] - y[4] - y[6] - uin2;
  const double ud2 = -y[3] + y[5] - y[6] - uin2;
  const double ud3 = +y[3] + y[4] + y[6] + uin2;
  const double ud4 = -y[2] - y[5] + y[6] + uin2;
  const double qpud1 = gamma * delta * exp (delta * ud1);
  const double qpud2 = gamma * delta * exp (delta * ud2);
  const double qpud3 = gamma * delta * exp (delta * ud3);
  const double qpud4 = gamma * delta * exp (delta * ud4);

  extern int nje;
  size_t i;

  nje += 1;

  for (i = 0; i < NRINGMOD * NRINGMOD; i++)
    {
      dfdy[i] = 0.0;
    }

  dfdy[0 * NRINGMOD + 0] = -1 / (c * r);
  dfdy[0 * NRINGMOD + 7] = 1 / c;
  dfdy[0 * NRINGMOD + 9] = -0.5 / c;
  dfdy[0 * NRINGMOD + 10] = -dfdy[0 * NRINGMOD + 9];
  dfdy[0 * NRINGMOD + 13] = dfdy[0 * NRINGMOD + 7];
  dfdy[1 * NRINGMOD + 1] = dfdy[0 * NRINGMOD + 0];
  dfdy[1 * NRINGMOD + 8] = dfdy[0 * NRINGMOD + 7];
  dfdy[1 * NRINGMOD + 11] = dfdy[0 * NRINGMOD + 9];
  dfdy[1 * NRINGMOD + 12] = dfdy[0 * NRINGMOD + 10];
  dfdy[1 * NRINGMOD + 14] = dfdy[0 * NRINGMOD + 13];
  dfdy[2 * NRINGMOD + 2] = (-qpud1 - qpud4) / cs;
  dfdy[2 * NRINGMOD + 4] = qpud1 / cs;
  dfdy[2 * NRINGMOD + 5] = -qpud4 / cs;
  dfdy[2 * NRINGMOD + 6] = (qpud1 + qpud4) / cs;
  dfdy[2 * NRINGMOD + 9] = 1 / cs;
  dfdy[3 * NRINGMOD + 3] = (-qpud2 - qpud3) / cs;
  dfdy[3 * NRINGMOD + 4] = -qpud3 / cs;
  dfdy[3 * NRINGMOD + 5] = qpud2 / cs;
  dfdy[3 * NRINGMOD + 6] = (-qpud2 - qpud3) / cs;
  dfdy[3 * NRINGMOD + 10] = -1 / cs;
  dfdy[4 * NRINGMOD + 2] = qpud1 / cs;
  dfdy[4 * NRINGMOD + 3] = -qpud3 / cs;
  dfdy[4 * NRINGMOD + 4] = (-qpud1 - qpud3) / cs;
  dfdy[4 * NRINGMOD + 6] = (-qpud1 - qpud3) / cs;
  dfdy[4 * NRINGMOD + 11] = 1 / cs;
  dfdy[5 * NRINGMOD + 2] = -qpud4 / cs;
  dfdy[5 * NRINGMOD + 3] = qpud2 / cs;
  dfdy[5 * NRINGMOD + 5] = (-qpud2 - qpud4) / cs;
  dfdy[5 * NRINGMOD + 6] = (qpud2 + qpud4) / cs;
  dfdy[5 * NRINGMOD + 12] = -1 / cs;
  dfdy[6 * NRINGMOD + 2] = (qpud1 + qpud4) / cp;
  dfdy[6 * NRINGMOD + 3] = (-qpud2 - qpud3) / cp;
  dfdy[6 * NRINGMOD + 4] = (-qpud1 - qpud3) / cp;
  dfdy[6 * NRINGMOD + 5] = (qpud2 + qpud4) / cp;
  dfdy[6 * NRINGMOD + 6] = (-qpud1 - qpud2 - qpud3 - qpud4 - 1 / rp) / cp;
  dfdy[7 * NRINGMOD + 0] = -1 / lh;
  dfdy[8 * NRINGMOD + 1] = dfdy[7 * NRINGMOD + 0];
  dfdy[9 * NRINGMOD + 0] = 0.5 / ls2;
  dfdy[9 * NRINGMOD + 2] = -1 / ls2;
  dfdy[9 * NRINGMOD + 9] = -rg2 / ls2;
  dfdy[10 * NRINGMOD + 0] = -0.5 / ls3;
  dfdy[10 * NRINGMOD + 3] = 1 / ls3;
  dfdy[10 * NRINGMOD + 10] = -rg3 / ls3;
  dfdy[11 * NRINGMOD + 1] = dfdy[9 * NRINGMOD + 0];
  dfdy[11 * NRINGMOD + 4] = dfdy[9 * NRINGMOD + 2];
  dfdy[11 * NRINGMOD + 11] = dfdy[9 * NRINGMOD + 9];
  dfdy[12 * NRINGMOD + 1] = dfdy[10 * NRINGMOD + 0];
  dfdy[12 * NRINGMOD + 5] = dfdy[10 * NRINGMOD + 3];
  dfdy[12 * NRINGMOD + 12] = dfdy[10 * NRINGMOD + 10];
  dfdy[13 * NRINGMOD + 0] = -1 / ls1;
  dfdy[13 * NRINGMOD + 13] = -(ri + rg1) / ls1;
  dfdy[14 * NRINGMOD + 1] = dfdy[13 * NRINGMOD + 0];
  dfdy[14 * NRINGMOD + 14] = -(rc + rg1) / ls1;

  for (i = 0; i < NRINGMOD; i++)
    {
      dfdt[i] = 0.0;
    }

  return GSL_SUCCESS;
}

gsl_odeiv2_system rhs_func_ringmod = {
  rhs_ringmod,
  jac_ringmod,
  NRINGMOD,
  NULL
};

/**********************************************************/
/* Functions for carrying out tests                       */
/**********************************************************/

void
test_odeiv_stepper (const gsl_odeiv2_step_type * T,
                    const gsl_odeiv2_system * sys, const double h,
                    const double t, const char desc[], const double ystart[],
                    const double yfin[], const double relerr)
{
  /* tests stepper T with one fixed length step advance of system sys
     and compares the result with given values yfin
   */

  double y[MAXEQ] = { 0.0 };
  double yerr[MAXEQ] = { 0.0 };
  double scale_abs[MAXEQ];
  size_t ne = sys->dimension;
  size_t i;
  gsl_odeiv2_driver *d;

  for (i = 0; i < MAXEQ; i++)
    {
      scale_abs[i] = 1.0;
    }

  d = gsl_odeiv2_driver_alloc_scaled_new (sys, T, h, relerr, relerr,
                                          1.0, 0.0, scale_abs);

  DBL_MEMCPY (y, ystart, MAXEQ);

  {
    int s = gsl_odeiv2_step_apply (d->s, t, h, y, yerr, 0, 0, sys);

    if (s != GSL_SUCCESS)
      {
        gsl_test (s, "test_odeiv_stepper: %s step_apply returned %d", desc,
                  s);
      }
  }

  for (i = 0; i < ne; i++)
    {
      gsl_test_rel (y[i], yfin[i], relerr,
                    "%s %s step(%d)", gsl_odeiv2_step_name (d->s), desc, i);
    }

  gsl_odeiv2_driver_free (d);
}

void
test_stepper (const gsl_odeiv2_step_type * T)
{
  /* Tests stepper T with a step of selected systems */

  double y[MAXEQ] = { 0.0 };
  double yfin[MAXEQ] = { 0.0 };

  /* Step length */
  double h;

  /* Required tolerance */
  double err_target;

  /* classic stiff */
  h = 1e-7;
  err_target = 1e-4;
  y[0] = 1.0;
  y[1] = 0.0;

  {
    const double e1 = exp (-h);
    const double e2 = exp (-1000.0 * h);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }

  test_odeiv_stepper (T, &rhs_func_stiff, h, 0.0, "classic_stiff",
                      y, yfin, err_target);

  /* linear */
  h = 1e-1;
  err_target = 1e-10;
  y[0] = 0.58;
  yfin[0] = y[0] + 2 * h;
  test_odeiv_stepper (T, &rhs_func_lin, h, 0.0, "linear",
                      y, yfin, err_target);

  /* exponential */
  h = 1e-4;
  err_target = 1e-8;
  y[0] = exp (2.7);
  yfin[0] = exp (2.7 + h);
  test_odeiv_stepper (T, &rhs_func_exp, h, 2.7, "exponential",
                      y, yfin, err_target);
  /* cosine-sine */
  h = 1e-3;
  err_target = 1e-6;
  y[0] = cos (1.2);
  y[1] = sin (1.2);
  yfin[0] = cos (1.2 + h);
  yfin[1] = sin (1.2 + h);
  test_odeiv_stepper (T, &rhs_func_sin, h, 1.2, "cosine-sine",
                      y, yfin, err_target);
}

void
test_evolve_system (const gsl_odeiv2_step_type * T,
                    const gsl_odeiv2_system * sys,
                    double t0, double t1, double hstart,
                    double y[], double yfin[],
                    double err_target, const char *desc)
{
  /* Tests system sys with stepper T. Step length is controlled by
     error estimation from the stepper.
   */

  int steps = 0;
  size_t i;

  double t = t0;
  double h = hstart;

  /* Tolerance factor in testing errors */
  const double factor = 10;

  gsl_odeiv2_driver *d =
    gsl_odeiv2_driver_alloc_standard_new (sys, T, h, err_target, err_target,
                                          1.0, 0.0);

  double *y_orig = (double *) malloc (sys->dimension * sizeof (double));

  int s = 0;

  extern int nfe, nje;
  nfe = 0;
  nje = 0;

  while (t < t1)
    {
      double t_orig = t;
      memcpy (y_orig, y, sys->dimension * sizeof (double));
      s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, sys, &t, t1, &h, y);

#ifdef DEBUG
      printf ("test_evolve_system at: %.5e %.5e %.5e %d\n", t, y[0], y[1],
              gsl_odeiv2_step_order (d->s));
#endif

      if (s != GSL_SUCCESS)
        {
          /* check that t and y are unchanged */
          gsl_test_abs (t, t_orig, 0.0, "%s t must be restored on failure",
                        gsl_odeiv2_step_name (d->s));

          for (i = 0; i < sys->dimension; i++)
            {
              gsl_test_abs (y[i], y_orig[i], 0.0,
                            "%s y must be restored on failure",
                            gsl_odeiv2_step_name (d->s), desc, i);
            }

          if (sys != &rhs_func_xsin)
            {
              /* apart from xsin, other functions should not return errors */
              gsl_test (s, "%s evolve_apply returned %d",
                        gsl_odeiv2_step_name (d->s), s);
              break;
            }
        }

      if (steps > 100000)
        {
          gsl_test (GSL_EFAILED,
                    "%s evolve_apply reached maxiter",
                    gsl_odeiv2_step_name (d->s));
          break;
        }

      steps++;
    }

  gsl_test (s, "%s %s [%g,%g], %d steps (nfe %d, nje %d) completed",
            gsl_odeiv2_step_name (d->s), desc, t0, t1, steps, nfe, nje);

  /* err_target is target error of one step. Test if stepper has made
     larger error than (tolerance factor times) the number of steps
     times the err_target.
   */

  for (i = 0; i < sys->dimension; i++)
    {
      gsl_test_abs (y[i], yfin[i], factor * d->e->count * err_target,
                    "%s %s evolve(%d)", gsl_odeiv2_step_name (d->s), desc, i);
    }

  free (y_orig);
  gsl_odeiv2_driver_free (d);
}

int
sys_driver (const gsl_odeiv2_step_type * T,
            const gsl_odeiv2_system * sys,
            double t0, double t1, double hstart,
            double y[], double epsabs, double epsrel, const char desc[])
{
  /* This function evolves a system sys with stepper T from t0 to t1.
     Step length is varied via error control with possibly different
     absolute and relative error tolerances.
   */

  int s = 0;
  int steps = 0;

  double t = t0;
  double h = hstart;

  gsl_odeiv2_driver *d =
    gsl_odeiv2_driver_alloc_standard_new (sys, T, h, epsabs, epsrel,
                                          1.0, 0.0);

  extern int nfe, nje;
  nfe = 0;
  nje = 0;

  while (t < t1)
    {
      s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, sys, &t, t1, &h, y);

#ifdef DEBUG
      printf ("sys_driver at: %.5e %.5e %.5e %d\n", t, y[0], y[1],
              gsl_odeiv2_step_order (d->s));
#endif

      if (s != GSL_SUCCESS)
        {
          gsl_test (s, "sys_driver: %s evolve_apply returned %d",
                    gsl_odeiv2_step_name (d->s), s);
          break;
        }

      if (steps > 1e7)
        {
          gsl_test (GSL_EMAXITER,
                    "sys_driver: %s evolve_apply reached maxiter at t=%g",
                    gsl_odeiv2_step_name (d->s), t);
          s = GSL_EMAXITER;
          break;
        }

      steps++;
    }

  gsl_test (s, "%s %s [%g,%g], %d steps (nfe %d, nje %d) completed",
            gsl_odeiv2_step_name (d->s), desc, t0, t1, steps, nfe, nje);

  gsl_odeiv2_driver_free (d);

  return s;
}

void
test_evolve_linear (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Test linear evolve */

  double y[1];
  double yfin[1];

  y[0] = 1.0;
  yfin[0] = 9.0;
  test_evolve_system (T, &rhs_func_lin, 0.0, 4.0, h, y, yfin, err,
                      "linear[0,4]");
}

void
test_evolve_exp (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Test exponential evolve */

  double y[1];
  double yfin[1];

  y[0] = 1.0;
  yfin[0] = exp (2.0);
  test_evolve_system (T, &rhs_func_exp, 0.0, 2.0, h, y, yfin, err,
                      "exp[0,2]");
}

void
test_evolve_sin (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Test sinusoidal evolve */

  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  yfin[0] = cos (2.0);
  yfin[1] = sin (2.0);
  test_evolve_system (T, &rhs_func_sin, 0.0, 2.0, h, y, yfin, err,
                      "sine[0,2]");
}

void
test_evolve_xsin (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Test sinusoidal evolve including a failing window */

  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  yfin[0] = cos (2.0);
  yfin[1] = sin (2.0);
  rhs_xsin_reset = 1;
  jac_xsin_reset = 1;
  test_evolve_system (T, &rhs_func_xsin, 0.0, 2.0, h, y, yfin, err,
                      "sine[0,2] w/errors");
}

void
test_evolve_stiff1 (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Test classical stiff problem evolve, t=[0,1] */

  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 1.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err,
                      "stiff[0,1]");
}

void
test_evolve_stiff5 (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Test classical stiff problem evolve, t=[0,5] */

  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 5.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 5.0, h, y, yfin, err,
                      "stiff[0,5]");
}

void
test_evolve_negative_h (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Test evolution in negative direction */

  /* Tolerance factor in testing errors */
  const double factor = 10;

  const gsl_odeiv2_system sys = rhs_func_cos;

  double t = 0;
  double t1 = -4.0;

  double y = 0.0;
  double yfin = sin (t1);
  gsl_odeiv2_driver *d;

  /* Make initial h negative */
  h = -fabs (h);
  d = gsl_odeiv2_driver_alloc_standard_new (&sys, T, h, err, err, 1.0, 0.0);

  while (t > t1)
    {
      int status = gsl_odeiv2_evolve_apply (d->e, d->c, d->s,
                                            &sys, &t, t1, &h, &y);

      if (status != GSL_SUCCESS)
        {
          gsl_test (status, "%s evolve_apply returned %d for negative h",
                    gsl_odeiv2_step_name (d->s), status);
          break;
        }
    }

  gsl_test_abs (y, yfin, factor * d->e->count * err,
                "%s evolution with negative h", gsl_odeiv2_step_name (d->s));

  gsl_odeiv2_driver_free (d);
}

void
test_broken (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Check for gsl_odeiv2_evolve_apply. The user function fails at
     t>=10, which in this test causes step size to decrease below
     machine precision and return with a failure code.
   */

  /* Tolerance factor in testing errors */
  const double factor = 10;

  const gsl_odeiv2_system sys = rhs_func_broken;

  gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys, T, h, err, err);

  double t = 0;
  double t1 = 100.0;

  double y = 0.0;
  const double final_max_h = GSL_DBL_EPSILON;

  int status;
  while (t < t1)
    {
      status =
        gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, t1, &h, &y);

      if (status != GSL_SUCCESS)
        {
          break;
        }
    }

  gsl_test_abs (h + final_max_h, final_max_h, factor * final_max_h,
                "%s test_broken: step size at break point",
                gsl_odeiv2_step_name (d->s));

  gsl_test_abs (t, 10.0, factor * err,
                "%s test_broken: point of break",
                gsl_odeiv2_step_name (d->s));

  /* GSL_FAILURE results from stepper failure, 123 from user function */

  if (status != GSL_FAILURE && status != 123)
    {
      gsl_test (status, "%s test_broken: evolve return value %d",
                gsl_odeiv2_step_name (d->s), status);
    }
  else
    {
      gsl_test (GSL_SUCCESS, "%s test_broken: evolve return value %d",
                gsl_odeiv2_step_name (d->s), status);
    }

  gsl_odeiv2_driver_free (d);
}

void
test_stepsize_fail (const gsl_odeiv2_step_type * T, double h)
{
  /* Check for gsl_odeiv2_evolve_apply. The user function works
     apparently fine (returns GSL_SUCCESS) but is faulty in this case.
     gsl_odeiv2_evolve_apply therefore decreases the step-size below
     machine precision and therefore stops with GSL_FAILURE.
   */

  /* Error tolerance */
  const double epsabs = 1e-16;
  const double epsrel = 1e-6;

  /* Tolerance factor in testing errors */
  const double factor = 10;

  const double final_max_h = GSL_DBL_EPSILON;
  const gsl_odeiv2_system sys = rhs_func_br;

  gsl_odeiv2_driver *d =
    gsl_odeiv2_driver_alloc_y_new (&sys, T, h, epsabs, epsrel);

  double t = 1.0;
  double t1 = 1e5;
  double y = 0.0;
  int status;

  while (t < t1)
    {
      status =
        gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, t1, &h, &y);

      if (status != GSL_SUCCESS)
        {
          break;
        }
    }

  gsl_test_abs (h + final_max_h, final_max_h, factor * final_max_h,
                "%s test_stepsize_fail: step size at break point",
                gsl_odeiv2_step_name (d->s));

  gsl_test_abs (t, 1.0, 1e-6,
                "%s test_stepsize_fail: point of break",
                gsl_odeiv2_step_name (d->s));

  gsl_test_int (status, GSL_FAILURE,
                "%s test_stepsize_fail: evolve return value",
                gsl_odeiv2_step_name (d->s));

  gsl_odeiv2_driver_free (d);
}

void
test_user_break (const gsl_odeiv2_step_type * T, double h)
{
  /* Tests for user signaled immediate break */

  const double tol = 1e-8;

  gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin_ub, T,
                                                        h, tol, tol);

  double y[] = { 1.0, 0.0 };
  double t = 0.0;
  const double t1 = 8.25;

  int s = gsl_odeiv2_driver_apply (d, &t, t1, y);

  gsl_test ((s - GSL_EBADFUNC), "%s test_user_break returned %d",
            gsl_odeiv2_step_name (d->s), s);

  gsl_odeiv2_driver_free (d);
}

void
test_stepfn (const gsl_odeiv2_step_type * T)
{
  /* Test evolve on a small derivative change at t=0 */

  double epsabs = 1e-16;
  double epsrel = 1e-6;

  const gsl_odeiv2_system sys = rhs_func_stepfn;

  double t = 0.0;
  double h = 1e-6;
  double y = 0.0;
  int i = 0;
  int status = 0;

  gsl_odeiv2_driver *d =
    gsl_odeiv2_driver_alloc_y_new (&sys, T, h, epsabs, epsrel);

  while (t < 2 && i < 1000000)
    {
      status =
        gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, 2, &h, &y);
#ifdef DEBUG
      printf ("test_stepfn at: i=%d status=%d t=%g h=%g y=%g\n",
			  i, status, t, h, y);
#endif

      if (status != GSL_SUCCESS)
        break;

      i++;
    }

  gsl_test (status, "evolve step function, return value (stepfn/%s): %d",
            gsl_odeiv2_step_name (d->s), status);

  gsl_test_abs (t, 2, 1e-16, "evolve step function, t (stepfn/%s)",
                gsl_odeiv2_step_name (d->s));

  gsl_test_rel (y, 1, epsrel, "evolve step function, y (stepfn/%s)",
                gsl_odeiv2_step_name (d->s));

  gsl_odeiv2_driver_free (d);
}

void
test_stepfn2 (const gsl_odeiv2_step_type * T)
{
  /* Test evolve on a large derivative change at t=0 */

  double epsabs = 1e-16;
  double epsrel = 1e-6;

  const gsl_odeiv2_system sys = rhs_func_stepfn2;

  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;

  int i = 0;
  int status;
  const int maxiter = 100000;

  gsl_odeiv2_driver *d =
    gsl_odeiv2_driver_alloc_yp_new (&sys, T, h, epsabs, epsrel);

  while (t < 1.0 && i < maxiter)
    {
      status =
        gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, 1.0, &h, &y);
#ifdef DEBUG
      printf ("test_stepfn2 at: i=%d status=%d t=%g h=%g y=%g\n",
			  i, status, t, h, y);
#endif
      if (status != GSL_SUCCESS)
        break;

      i++;
    }

  if (i >= maxiter)
    printf ("FAIL: evolve big step function, (stepfn2/%s) reached maxiter\n",
            gsl_odeiv2_step_name (d->s));

  gsl_test_abs (t, 1.0, 1e-16, "evolve big step function, t (stepfn2/%s)",
                gsl_odeiv2_step_name (d->s));
  gsl_test_rel (y, 1e300, epsrel, "evolve big step function, y (stepfn2/%s)",
                gsl_odeiv2_step_name (d->s));

  gsl_odeiv2_driver_free (d);
}

void
test_nonstiff_problems (void)
{
  /* Compares output of non-stiff (or only moderately stiff) problems
     with several steppers
   */

  const gsl_odeiv2_step_type *steppers[MAXNS];

  /* Required error tolerance for each stepper. */
  double err_target[MAXNS];

  /* initial values for each ode-solver */
  double y[MAXEQ * MAXNS];

  size_t i, k, p;

  /* Number of problems to test */
#define CONST_NONSTIFF_NPROB 4

  /* Problems, their names and number of equations */
  const gsl_odeiv2_system *prob[] =
    { &rhs_func_vl, &rhs_func_vanderpol, &rhs_func_stifftrig,
    &rhs_func_brusselator
  };
  const char *probname[] =
    { "volterra-lotka", "vanderpol", "stifftrig", "brusselator" };
  const size_t sd[] = { 2, 2, 1, 2 };

  /* Integration interval for problems */
  const double start[CONST_NONSTIFF_NPROB] = { 0.0 };
  const double end[] = { 9.0, 100.0, 1.5, 20.0 };

  const double epsabs = 1e-8;
  const double epsrel = 1e-8;
  const double initstepsize = 1e-5;

  /* Steppers */

  steppers[0] = gsl_odeiv2_step_rk4;
  err_target[0] = 1e-6;
  steppers[1] = gsl_odeiv2_step_rk2;
  err_target[1] = 1e-6;
  steppers[2] = gsl_odeiv2_step_rkf45;
  err_target[2] = 1e-6;
  steppers[3] = gsl_odeiv2_step_rkck;
  err_target[3] = 1e-6;
  steppers[4] = gsl_odeiv2_step_rk8pd;
  err_target[4] = 1e-6;
  steppers[5] = gsl_odeiv2_step_rk1imp;
  err_target[5] = 1e-3;
  steppers[6] = gsl_odeiv2_step_rk2imp;
  err_target[6] = 1e-5;
  steppers[7] = gsl_odeiv2_step_rk4imp;
  err_target[7] = 1e-6;
  steppers[8] = gsl_odeiv2_step_bsimp;
  err_target[8] = 1e-6;
  steppers[9] = gsl_odeiv2_step_msadams;
  err_target[9] = 1e-5;
  steppers[10] = gsl_odeiv2_step_msbdf;
  err_target[10] = 1e-5;
  steppers[11] = 0;

  /* Loop over problems */

  for (p = 0; p < CONST_NONSTIFF_NPROB; p++)
    {
      /* Initialize */

      for (i = 0; i < MAXNS * MAXEQ; i++)
        {
          y[i] = 0.0;
        }

      for (i = 0; i < MAXNS; i++)
        {
          switch (p)
            {
            case 0:
              y[i * sd[p]] = 2.725;
              y[i * sd[p] + 1] = 1.0;
              break;
            case 1:
              y[i * sd[p]] = 1.0;
              y[i * sd[p] + 1] = 0.0;
              break;
            case 2:
              y[i * sd[p]] = 0.0;
              break;
            case 3:
              y[i * sd[p]] = 1.5;
              y[i * sd[p] + 1] = 3.0;
              break;
            default:
              gsl_test (GSL_EFAILED,
                        "test_nonstiff_problems: initialization error\n");
              return;
            }
        }

      /* Call each solver for the problem */

      for (i = 0; steppers[i] != 0; i++)
        {
          int s = sys_driver (steppers[i], prob[p], start[p], end[p],
                              initstepsize, &y[sd[p] * i],
                              epsabs, epsrel, probname[p]);

          if (s != GSL_SUCCESS)
            {
              gsl_test (s, "test_nonstiff_problems %s %s",
                        steppers[i]->name, probname[p]);
            }
        }

      /* Compare results */

      for (i = 1; steppers[i] != 0; i++)
        for (k = 0; k < sd[p]; k++)
          {
            const double val1 = y[k];
            const double val2 = y[sd[p] * i + k];
            gsl_test_rel (val1, val2,
                          (GSL_MAX (err_target[0], err_target[i])),
                          "%s/%s %s [%d]",
                          steppers[0]->name, steppers[i]->name,
                          probname[p], k);
          }
    }
}

void
test_stiff_problems (void)
{
  /* Compares output of stiff problems with several steppers */

  const gsl_odeiv2_step_type *steppers[MAXNS];

  /* Required error tolerance for each stepper. */
  double err_target[MAXNS];

  /* initial values for each ode-solver */
  double y[MAXEQ * MAXNS];

  size_t i, k, p;

  /* Number of problems to test */
#define CONST_STIFF_NPROB 3

  /* Problems, their names and number of equations */
  const gsl_odeiv2_system *prob[] =
    { &rhs_func_oregonator, &rhs_func_e5, &rhs_func_robertson };
  const char *probname[] = { "oregonator", "e5", "robertson" };
  const size_t sd[] = { 3, 4, 3 };

  /* Integration interval for problems */
  const double start[CONST_STIFF_NPROB] = { 0.0 };
  const double end[] = { 360.0, 1e4, 1e4 };

  const double epsabs = 1e-40;
  const double epsrel = 1e-7;
  const double initstepsize = 1e-5;

  /* Steppers */

  steppers[0] = gsl_odeiv2_step_bsimp;
  err_target[0] = 1e-6;
  steppers[1] = gsl_odeiv2_step_rk1imp;
  err_target[1] = 5e-3;
  steppers[2] = gsl_odeiv2_step_rk2imp;
  err_target[2] = 5e-5;
  steppers[3] = gsl_odeiv2_step_rk4imp;
  err_target[3] = 5e-5;
  steppers[4] = gsl_odeiv2_step_msbdf;
  err_target[4] = 1e-4;
  steppers[5] = 0;

  /* Loop over problems */

  for (p = 0; p < CONST_STIFF_NPROB; p++)
    {
      /* Initialize */

      for (i = 0; i < MAXNS * MAXEQ; i++)
        {
          y[i] = 0.0;
        }

      for (i = 0; i < MAXNS; i++)
        {
          switch (p)
            {
            case 0:
              y[i * sd[p]] = 1.0;
              y[i * sd[p] + 1] = 2.0;
              y[i * sd[p] + 2] = 3.0;
              break;
            case 1:
              y[i * sd[p]] = 1.76e-3;
              y[i * sd[p] + 1] = 0.0;
              y[i * sd[p] + 2] = 0.0;
              y[i * sd[p] + 3] = 0.0;
              break;
            case 2:
              y[i * sd[p]] = 1.0;
              y[i * sd[p] + 1] = 0.0;
              y[i * sd[p] + 2] = 0.0;
              break;
            default:
              gsl_test (GSL_EFAILED,
                        "test_stiff_problems: initialization error\n");
              return;
            }
        }

      /* Call each solver for the problem */

      for (i = 0; steppers[i] != 0; i++)
        {
          int s = sys_driver (steppers[i], prob[p], start[p], end[p],
                              initstepsize, &y[sd[p] * i],
                              epsabs, epsrel, probname[p]);

          if (s != GSL_SUCCESS)
            {
              gsl_test (s, "test_stiff_problems %s %s",
                        steppers[i]->name, probname[p]);
            }
        }

      /* Compare results */

      for (i = 1; steppers[i] != 0; i++)
        for (k = 0; k < sd[p]; k++)
          {
            const double val1 = y[k];
            const double val2 = y[sd[p] * i + k];
            gsl_test_rel (val1, val2,
                          (GSL_MAX (err_target[0], err_target[i])),
                          "%s/%s %s [%d]",
                          steppers[0]->name, steppers[i]->name,
                          probname[p], k);
          }
    }
}

void
test_extreme_problems (void)
{
  /* Compares output of numerically particularly demanding problems
     with several steppers */

  const gsl_odeiv2_step_type *steppers[MAXNS];

  /* Required error tolerance for each stepper. */
  double err_target[MAXNS];

  /* initial values for each ode-solver */
  double y[MAXEQ * MAXNS];

  size_t i, k, p;

  /* Number of problems to test */
#define CONST_EXTREME_NPROB 3

  /* Problems, their names and number of equations */
  const gsl_odeiv2_system *prob[] =
    { &rhs_func_e5, &rhs_func_robertson, &rhs_func_ringmod };
  const char *probname[] = { "e5_bigt", "robertson_bigt", "ringmod" };
  const size_t sd[] = { 4, 3, NRINGMOD };

  /* Integration interval for problems */

  const double start[CONST_EXTREME_NPROB] = { 0.0 };
  const double end[CONST_EXTREME_NPROB] = { 1e11, 1e11, 1e-5 };

  const double epsabs[CONST_EXTREME_NPROB] =
    { 1e1 * GSL_DBL_MIN, 1e1 * GSL_DBL_MIN, 1e-12 };
  const double epsrel[CONST_EXTREME_NPROB] = { 1e-12, 1e-12, 1e-12 };
  const double initstepsize[CONST_EXTREME_NPROB] = { 1e-5, 1e-5, 1e-10 };

  /* Steppers */

  steppers[0] = gsl_odeiv2_step_bsimp;
  err_target[0] = 1e-3;
  steppers[1] = gsl_odeiv2_step_msbdf;
  err_target[1] = 1e-3;
  steppers[2] = 0;

  /* Loop over problems */

  for (p = 0; p < CONST_EXTREME_NPROB; p++)
    {
      /* Initialize */

      for (i = 0; i < MAXNS * MAXEQ; i++)
        {
          y[i] = 0.0;
        }

      for (i = 0; i < MAXNS; i++)
        {
          switch (p)
            {
            case 0:
              y[i * sd[p]] = 1.76e-3;
              y[i * sd[p] + 1] = 0.0;
              y[i * sd[p] + 2] = 0.0;
              y[i * sd[p] + 3] = 0.0;
              break;
            case 1:
              y[i * sd[p]] = 1.0;
              y[i * sd[p] + 1] = 0.0;
              y[i * sd[p] + 2] = 0.0;
              break;
            case 2:
              {
                size_t j;
                for (j = 0; j < NRINGMOD; j++)
                  {
                    y[i * sd[p] + j] = 0.0;
                  }
              }
              break;
            default:
              gsl_test (GSL_EFAILED,
                        "test_extreme_problems: initialization error\n");
              return;
            }
        }

      /* Call each solver for the problem */

      for (i = 0; steppers[i] != 0; i++)
        {
          int s = sys_driver (steppers[i], prob[p], start[p], end[p],
                              initstepsize[p], &y[sd[p] * i],
                              epsabs[p], epsrel[p], probname[p]);

          if (s != GSL_SUCCESS)
            {
              printf ("start=%.5e, initstepsize=%.5e\n", start[p],
                      initstepsize[p]);
              gsl_test (s, "test_extreme_problems %s %s", steppers[i]->name,
                        probname[p]);
            }
        }

      /* Compare results */

      for (i = 1; steppers[i] != 0; i++)
        for (k = 0; k < sd[p]; k++)
          {
            const double val1 = y[k];
            const double val2 = y[sd[p] * i + k];
            gsl_test_rel (val1, val2,
                          (GSL_MAX (err_target[0], err_target[i])),
                          "%s/%s %s [%d]",
                          steppers[0]->name, steppers[i]->name,
                          probname[p], k);
          }
    }
}

void
test_driver (const gsl_odeiv2_step_type * T)
{
  /* Tests for gsl_odeiv2_driver object */

  int s;
  const double tol = 1e-8;
  const double hmax = 1e-2;

  double y[] = { 1.0, 0.0 };
  double t = 0.0;
  const double tinit = 0.0;
  const double t1 = 8.25;
  const double t2 = 100;
  const double t3 = -2.5;
  const size_t minsteps = ceil (t1 / hmax);
  const size_t maxsteps = 20;
  const double hmin = 1e-10;

  const unsigned long int nfsteps = 100000;
  const double hfixed = 0.000025;

  gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T,
                                                        1e-3, tol, tol);
  gsl_odeiv2_driver_set_hmax (d, hmax);

  s = gsl_odeiv2_driver_apply (d, &t, t1, y);

  if (s != GSL_SUCCESS)
    {
      gsl_test (s, "%s test_driver apply returned %d",
                gsl_odeiv2_step_name (d->s), s);
    }

  /* Test that result is correct */

  gsl_test_rel (y[0], cos (t1), d->n * tol, "%s test_driver y0",
                gsl_odeiv2_step_name (d->s));
  gsl_test_rel (y[1], sin (t1), d->n * tol, "%s test_driver y1",
                gsl_odeiv2_step_name (d->s));

  /* Test that maximum step length is obeyed */

  if (d->n < minsteps)
    {
      gsl_test (1, "%s test_driver steps %d < minsteps %d \n",
                gsl_odeiv2_step_name (d->s), d->n, minsteps);
    }
  else
    {
      gsl_test (0, "%s test_driver max step length test",
                gsl_odeiv2_step_name (d->s));
    }

  /* Test changing integration direction from forward to backward */

  gsl_odeiv2_driver_reset_hstart (d, -1e-3);

  s = gsl_odeiv2_driver_apply (d, &t, tinit, y);

  if (s != GSL_SUCCESS)
    {
      gsl_test (s, "%s test_driver apply returned %d",
                gsl_odeiv2_step_name (d->s), s);
    }

  gsl_test_rel (y[0], cos (tinit), d->n * tol,
                "%s test_driver y0", gsl_odeiv2_step_name (d->s));
  gsl_test_rel (y[1], sin (tinit), d->n * tol,
                "%s test_driver y1", gsl_odeiv2_step_name (d->s));

  gsl_odeiv2_driver_free (d);

  /* Test that maximum number of steps is obeyed */

  d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T, 1e-3, tol, tol);
  gsl_odeiv2_driver_set_hmax (d, hmax);
  gsl_odeiv2_driver_set_nmax (d, maxsteps);

  s = gsl_odeiv2_driver_apply (d, &t, t2, y);

  if (d->n != maxsteps + 1)
    {
      gsl_test (1, "%s test_driver steps %d, expected %d",
                gsl_odeiv2_step_name (d->s), d->n, maxsteps + 1);
    }
  else
    {
      gsl_test (0, "%s test_driver max steps test",
                gsl_odeiv2_step_name (d->s));
    }

  gsl_odeiv2_driver_free (d);

  /* Test that minimum step length is obeyed */

  d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_broken, T, 1e-3, tol, tol);

  gsl_odeiv2_driver_set_hmin (d, hmin);
  y[0] = 0.0;
  t = 0.0;

  s = gsl_odeiv2_driver_apply (d, &t, t2, y);

  if (s != GSL_ENOPROG)
    {
      gsl_test (1, "%s test_driver min step test returned %d",
                gsl_odeiv2_step_name (d->s), s);
    }
  else
    {
      gsl_test (0, "%s test_driver min step test",
                gsl_odeiv2_step_name (d->s));
    }

  gsl_odeiv2_driver_free (d);

  /* Test negative integration direction */

  d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T, -1e-3, tol, tol);

  y[0] = 1.0;
  y[1] = 0.0;
  t = 2.5;

  s = gsl_odeiv2_driver_apply (d, &t, t3, y);

  {
    const double tol = 1e-3;
    const double test = fabs (t - t3);
    const double val = fabs (sin (-5.0) - y[1]);

    if (test > GSL_DBL_EPSILON)
      {
        gsl_test (1,
                  "%s test_driver negative dir diff %e, expected less than %e",
                  gsl_odeiv2_step_name (d->s), test, GSL_DBL_EPSILON);
      }
    else if (val > tol)
      {
        gsl_test (1,
                  "%s test_driver negative dir val %e, expected less than %e",
                  gsl_odeiv2_step_name (d->s), val, tol);
      }
    else
      {
        gsl_test (s, "%s test_driver negative direction test",
                  gsl_odeiv2_step_name (d->s));
      }
  }

  /* Test driver_apply_fixed_step */

  gsl_odeiv2_driver_reset_hstart (d, 1e-3);
  s = gsl_odeiv2_driver_apply_fixed_step (d, &t, hfixed, nfsteps, y);

  if (s != GSL_SUCCESS)
    {
      gsl_test (s, "%s test_driver apply_fixed_step returned %d",
                gsl_odeiv2_step_name (d->s), s);
    }

  {
    const double tol = 1e-3;
    const double val = fabs (sin (-2.5) - y[1]);

    if (fabs (t) > nfsteps * GSL_DBL_EPSILON)
      {
        gsl_test (1,
                  "%s test_driver apply_fixed_step t %e, expected less than %e",
                  gsl_odeiv2_step_name (d->s), fabs (t),
                  nfsteps * GSL_DBL_EPSILON);
      }
    else if (val > tol)
      {
        gsl_test (1,
                  "%s test_driver apply_fixed_step val %e, expected less than %e",
                  gsl_odeiv2_step_name (d->s), val, tol);
      }
    else
      {
        gsl_test (s, "%s test_driver apply_fixed_step test",
                  gsl_odeiv2_step_name (d->s));
      }
  }

  gsl_odeiv2_driver_free (d);
}

void
benchmark_precision (void)
{
  /* Tests steppers with several error tolerances and evaluates their
     performance and precision.
   */

  /* Maximum number of tolerances to be tested */
#define MAXNT 12

  const gsl_odeiv2_step_type *steppers[MAXNS];

  /* Required error tolerance for each stepper. */
  double err_target[MAXNS];

  /* initial values for each ode-solver */
  double y[MAXEQ * MAXNS * MAXNT];

  /* precise results from e.g. analytical solution */
  double yres[MAXEQ];

  size_t i, j, k, p;

  /* Number of problems to test */
#define CONST_BENCHMARK_PRECISION_NPROB 3

  /* Problems, their names and number of equations */
  const gsl_odeiv2_system *prob[] =
    { &rhs_func_sin, &rhs_func_exp, &rhs_func_stiff };
  const char *probname[] =
    { "rhs_func_sin", "rhs_func_exp", "rhs_func_stiff" };
  const size_t sd[] = { 2, 1, 2 };

  /* Integration interval for problems */
  const double start[] = { 0.0, 0.0, 0.0 };
  const double end[] = { 2.4, 8.4, 1.2 };

  const double epsabs[] = { 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14, 0 };
  const double epsrel[] = { 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14, 0 };
  const double initstepsize = 1e-5;

  /* number of function and jacobian evaluations */
  extern int nfe, nje;
  int nnfe[MAXNS * MAXNT] = { 0.0 };
  int nnje[MAXNS * MAXNT] = { 0.0 };

  /* Steppers */

  steppers[0] = gsl_odeiv2_step_rk4;
  err_target[0] = 1e-6;
  steppers[1] = gsl_odeiv2_step_rk2;
  err_target[1] = 1e-6;
  steppers[2] = gsl_odeiv2_step_rkf45;
  err_target[2] = 1e-6;
  steppers[3] = gsl_odeiv2_step_rkck;
  err_target[3] = 1e-6;
  steppers[4] = gsl_odeiv2_step_rk8pd;
  err_target[4] = 1e-6;
  steppers[5] = gsl_odeiv2_step_rk1imp;
  err_target[5] = 1e-3;
  steppers[6] = gsl_odeiv2_step_rk2imp;
  err_target[6] = 1e-5;
  steppers[7] = gsl_odeiv2_step_rk4imp;
  err_target[7] = 1e-6;
  steppers[8] = gsl_odeiv2_step_bsimp;
  err_target[8] = 1e-6;
  steppers[9] = gsl_odeiv2_step_msadams;
  err_target[9] = 1e-5;
  steppers[10] = gsl_odeiv2_step_msbdf;
  err_target[10] = 1e-5;
  steppers[11] = 0;

  /* Loop over problems */

  for (p = 0; p < CONST_BENCHMARK_PRECISION_NPROB; p++)
    {
      /* Initialize */

      for (i = 0; i < MAXNS * MAXEQ * MAXNT; i++)
        {
          y[i] = 0.0;
        }

      for (i = 0; i < MAXNS * MAXNT; i++)
        {
          switch (p)
            {
            case 0:
              y[i * sd[p]] = 1.0;
              y[i * sd[p] + 1] = 0.0;
              yres[0] = cos (2.4);
              yres[1] = sin (2.4);
              break;
            case 1:
              y[i * sd[p]] = 1.0;
              yres[0] = exp (8.4);
              break;
            case 2:
              y[i * sd[p]] = 1.0;
              y[i * sd[p] + 1] = 0.0;
              yres[0] = 2 * exp (-1.2) - exp (-1000 * 1.2);
              yres[1] = -exp (-1.2) + exp (-1000 * 1.2);
              break;
            default:
              gsl_test (GSL_EFAILED,
                        "benchmark_precision: initialization error\n");
              return;
            }
        }

      /* Call each solver for the problem */

      for (i = 0; steppers[i] != 0; i++)
        for (j = 0; epsrel[j] != 0; j++)
          {
            int s = sys_driver (steppers[i], prob[p], start[p], end[p],
                                initstepsize, &y[sd[p] * (j + i * MAXNT)],
                                epsabs[j], epsrel[j], probname[p]);

            if (s != GSL_SUCCESS)
              {
                for (k = 0; k < sd[p]; k++)
                  {
                    y[sd[p] * (j + i * MAXNT) + k] = GSL_NAN;
                  }
              }

            else
              {
                nnfe[j + i * MAXNT] = nfe;
                nnje[j + i * MAXNT] = nje;
              }
          }

      /* Print results */

      printf
        ("benchmark_precision: diff = (y_true - y) / y_true with %s\n          epsrel:   ",
         probname[p]);

      for (i = 0; epsrel[i] != 0.0; i++)
        {
          printf ("%12.0e ", epsrel[i]);
        }

      printf ("\n");

      for (i = 0; steppers[i] != 0; i++)
        {
          for (k = 0; k < sd[p]; k++)
            {
              printf ("%8s diff[%d]:   ", steppers[i]->name, (int) k);

              for (j = 0; epsrel[j] != 0; j++)
                {
                  const double diff =
                    (yres[k] - y[sd[p] * (j + i * MAXNT) + k]) / yres[k];
                  printf ("%12.5e ", diff);
                }

              printf ("\n");
            }
        }

      printf ("\n");

      /* Print number of function/jacobian evaluations */

      printf
        ("benchmark_precision: number of function/jacobian evaluations with %s\n          epsrel:   ",
         probname[p]);

      for (i = 0; epsrel[i] != 0.0; i++)
        {
          printf ("%12.0e          ", epsrel[i]);
        }

      printf ("\n");

      for (i = 0; steppers[i] != 0; i++)
        {
          printf ("%8s nfe/nje: ", steppers[i]->name);

          for (j = 0; epsrel[j] != 0; j++)
            {
              printf ("%9d %9d | ", nnfe[j + i * MAXNT], nnje[j + i * MAXNT]);
            }

          printf ("\n");
        }
      printf ("\n");

    }
}

void
test_evolve_temp (const gsl_odeiv2_step_type * T, double h, double err)
{
  /* Temporary test */

  double y[3];
  double yfin[3];

  y[0] = 1.0;
  y[1] = 0.0;
  y[2] = 0.0;
  yfin[0] = 0;
  yfin[1] = 0;
  yfin[2] = 0;
  test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err, "temp");
}

/**********************************************************/
/* Main function                                          */
/**********************************************************/

int
main (void)
{

  /* Benchmark routine to compare stepper performance */
  /* benchmark_precision(); return 0; */

  /* Test single problem, for debugging purposes */
  /* test_evolve_temp (gsl_odeiv2_step_msadams, 1e-3, 1e-8); return 0; */
  int i;

  struct ptype
  {
    const gsl_odeiv2_step_type *type;
    double h;
  }
  p[MAXNS];

  struct ptype explicit_stepper[MAXNS];

  p[0].type = gsl_odeiv2_step_rk4;
  p[0].h = 1.0e-3;
  p[1].type = gsl_odeiv2_step_rk2;
  p[1].h = 1.0e-3;
  p[2].type = gsl_odeiv2_step_rkf45;
  p[2].h = 1.0e-3;
  p[3].type = gsl_odeiv2_step_rkck;
  p[3].h = 1.0e-3;
  p[4].type = gsl_odeiv2_step_rk8pd;
  p[4].h = 1.0e-3;
  p[5].type = gsl_odeiv2_step_rk1imp;
  p[5].h = 1.0e-3;
  p[6].type = gsl_odeiv2_step_rk2imp;
  p[6].h = 1.0e-3;
  p[7].type = gsl_odeiv2_step_rk4imp;
  p[7].h = 1.0e-3;
  p[8].type = gsl_odeiv2_step_bsimp;
  p[8].h = 1.0e-3;
  p[9].type = gsl_odeiv2_step_msadams;
  p[9].h = 1.0e-3;
  p[10].type = gsl_odeiv2_step_msbdf;
  p[10].h = 1.0e-3;
  p[11].type = 0;

  gsl_ieee_env_setup ();

  /* Basic tests for all steppers */

  for (i = 0; p[i].type != 0; i++)
    {
      test_stepper (p[i].type);
      test_evolve_linear (p[i].type, p[i].h, 1e-10);
      test_evolve_exp (p[i].type, p[i].h, 1e-5);
      test_evolve_sin (p[i].type, p[i].h, 1e-8);
      test_evolve_xsin (p[i].type, p[i].h, 1e-8);
      test_evolve_xsin (p[i].type, 0.1, 1e-8);
      test_evolve_stiff1 (p[i].type, p[i].h, 1e-7);
      test_evolve_stiff5 (p[i].type, p[i].h, 1e-7);
      test_evolve_negative_h (p[i].type, p[i].h, 1e-7);
      test_broken (p[i].type, p[i].h, 1e-8);
      test_stepsize_fail (p[i].type, p[i].h);
      test_user_break (p[i].type, p[i].h);
      test_driver (p[i].type);
    }

  /* Derivative test for explicit steppers */

  explicit_stepper[0].type = gsl_odeiv2_step_rk4;
  explicit_stepper[0].h = 1.0e-3;
  explicit_stepper[1].type = gsl_odeiv2_step_rk2;
  explicit_stepper[1].h = 1.0e-3;
  explicit_stepper[2].type = gsl_odeiv2_step_rkf45;
  explicit_stepper[2].h = 1.0e-3;
  explicit_stepper[3].type = gsl_odeiv2_step_rkck;
  explicit_stepper[3].h = 1.0e-3;
  explicit_stepper[4].type = gsl_odeiv2_step_rk8pd;
  explicit_stepper[4].h = 1.0e-3;
  explicit_stepper[5].type = gsl_odeiv2_step_msadams;
  explicit_stepper[5].h = 1.0e-3;
  explicit_stepper[6].type = 0;

  for (i = 0; explicit_stepper[i].type != 0; i++)
    {
      test_stepfn (explicit_stepper[i].type);
      test_stepfn2 (explicit_stepper[i].type);
    }

  /* Special tests */

  test_nonstiff_problems ();

  test_stiff_problems ();

  test_extreme_problems ();

  exit (gsl_test_summary ());
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*******************************************************************************
 * Copyright (c) 2015-2018 Skymind, Inc.
 *
 * This program and the accompanying materials are made available under the
 * terms of the Apache License, Version 2.0 which is available at
 * https://www.apache.org/licenses/LICENSE-2.0.
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
 * License for the specific language governing permissions and limitations
 * under the License.
 *
 * SPDX-License-Identifier: Apache-2.0
 ******************************************************************************/

//
//  @author raver119@gmail.com
//

#ifndef LIBND4J_BLAS_HELPER_H
#define LIBND4J_BLAS_HELPER_H

#include <system/pointercast.h>
#include <types/float16.h>
#include <cblas.h>
#include <helpers/logger.h>

#ifdef _WIN32
#define CUBLASWINAPI __stdcall
#define CUSOLVERAPI __stdcall
#else
#define CUBLASWINAPI
#define CUSOLVERAPI
#endif

namespace sd {
    typedef enum{
        CUBLAS_STATUS_SUCCESS         =0,
        CUBLAS_STATUS_NOT_INITIALIZED =1,
        CUBLAS_STATUS_ALLOC_FAILED    =3,
        CUBLAS_STATUS_INVALID_VALUE   =7,
        CUBLAS_STATUS_ARCH_MISMATCH   =8,
        CUBLAS_STATUS_MAPPING_ERROR   =11,
        CUBLAS_STATUS_EXECUTION_FAILED=13,
        CUBLAS_STATUS_INTERNAL_ERROR  =14,
        CUBLAS_STATUS_NOT_SUPPORTED   =15,
        CUBLAS_STATUS_LICENSE_ERROR   =16
    } cublasStatus_t;

    typedef enum {
        CUBLAS_OP_N=0,
        CUBLAS_OP_T=1,
        CUBLAS_OP_C=2
    } cublasOperation_t;

    struct cublasContext;
    typedef struct cublasContext *cublasHandle_t;

    typedef enum
    {
	    CUDA_R_16F= 2,  /* real as a half */
	    CUDA_C_16F= 6,  /* complex as a pair of half numbers */
	    CUDA_R_32F= 0,  /* real as a float */
	    CUDA_C_32F= 4,  /* complex as a pair of float numbers */
	    CUDA_R_64F= 1,  /* real as a double */
	    CUDA_C_64F= 5,  /* complex as a pair of double numbers */
	    CUDA_R_8I = 3,  /* real as a signed char */
	    CUDA_C_8I = 7,  /* complex as a pair of signed char numbers */
	    CUDA_R_8U = 8,  /* real as a unsigned char */
	    CUDA_C_8U = 9,  /* complex as a pair of unsigned char numbers */
	    CUDA_R_32I= 10, /* real as a signed int */
	    CUDA_C_32I= 11, /* complex as a pair of signed int numbers */
	    CUDA_R_32U= 12, /* real as a unsigned int */
	    CUDA_C_32U= 13  /* complex as a pair of unsigned int numbers */
    } cublasDataType_t;

    typedef void (*CblasSgemv)(CBLAS_ORDER Layout,
                 CBLAS_TRANSPOSE TransA, int M, int N,
                 float alpha, float *A, int lda,
                 float *X, int incX, float beta,
                 float *Y, int incY);

    typedef void (*CblasDgemv)(CBLAS_ORDER Layout,
                 CBLAS_TRANSPOSE TransA, int M, int N,
                 double alpha, double *A, int lda,
                 double *X, int incX, double beta,
                 double *Y, int incY);

    typedef void (*CblasSgemm)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE TransA,
                     CBLAS_TRANSPOSE TransB, int M, int N,
                     int K, float alpha, float *A,
                     int lda, float *B, int ldb,
                     float beta, float *C, int ldc);

    typedef void (*CblasDgemm)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE TransA,
                     CBLAS_TRANSPOSE TransB, int M, int N,
                     int K, double alpha, double *A,
                     int lda, double *B, int ldb,
                     double beta, double *C, int ldc);

    typedef void (*CblasSgemmBatch)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE *TransA_Array,
                           CBLAS_TRANSPOSE *TransB_Array, int *M_Array, int *N_Array,
                           int *K_Array, float *alpha_Array, float **A_Array,
                           int *lda_Array, float **B_Array, int *ldb_Array,
                           float *beta_Array, float **C_Array, int *ldc_Array,
                           int group_count, int *group_size);

    typedef void (*CblasDgemmBatch)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE *TransA_Array,
                           CBLAS_TRANSPOSE *TransB_Array, int *M_Array, int *N_Array,
                           int *K_Array, double *alpha_Array, double **A_Array,
                           int *lda_Array, double **B_Array, int* ldb_Array,
                           double *beta_Array, double **C_Array, int *ldc_Array,
                           int group_count, int *group_size);

#ifdef LAPACK_ROW_MAJOR
#undef LAPACK_ROW_MAJOR
#endif

#ifdef LAPACK_COL_MAJOR
#undef LAPACK_COL_MAJOR
#endif
    enum LAPACK_LAYOUT { LAPACK_ROW_MAJOR=101, LAPACK_COL_MAJOR=102 };

    typedef int (*LapackeSgesvd)(LAPACK_LAYOUT matrix_layout, char jobu, char jobvt,
                           int m, int n, float* a, int lda,
                           float* s, float* u, int ldu, float* vt,
                           int ldvt, float* superb);

    typedef int (*LapackeDgesvd)(LAPACK_LAYOUT matrix_layout, char jobu, char jobvt,
                           int m, int n, double* a,
                           int lda, double* s, double* u, int ldu,
                           double* vt, int ldvt, double* superb);

    typedef int (*LapackeSgesdd)(LAPACK_LAYOUT matrix_layout, char jobz, int m,
                           int n, float* a, int lda, float* s,
                           float* u, int ldu, float* vt,
                           int ldvt);
    typedef int (*LapackeDgesdd)(LAPACK_LAYOUT matrix_layout, char jobz, int m,
                           int n, double* a, int lda, double* s,
                           double* u, int ldu, double* vt,
                           int ldvt);

    typedef cublasStatus_t (CUBLASWINAPI *CublasSgemv)(cublasHandle_t handle,
                                                      cublasOperation_t trans,
                                                      int m,
                                                      int n,
                                                      float *alpha, /* host or device pointer */
                                                      float *A,
                                                      int lda,
                                                      float *x,
                                                      int incx,
                                                      float *beta,  /* host or device pointer */
                                                      float *y,
                                                      int incy);

    typedef cublasStatus_t (CUBLASWINAPI *CublasDgemv)(cublasHandle_t handle,
                                                      cublasOperation_t trans,
                                                      int m,
                                                      int n,
                                                      double *alpha, /* host or device pointer */
                                                      double *A,
                                                      int lda,
                                                      double *x,
                                                      int incx,
                                                      double *beta, /* host or device pointer */
                                                      double *y,
                                                      int incy);

    typedef cublasStatus_t (CUBLASWINAPI *CublasHgemm)(cublasHandle_t handle,
                                                      cublasOperation_t transa,
                                                      cublasOperation_t transb,
                                                      int m,
                                                      int n,
                                                      int k,
                                                      __half *alpha, /* host or device pointer */
                                                      __half *A,
                                                      int lda,
                                                      __half *B,
                                                      int ldb,
                                                      __half *beta, /* host or device pointer */
                                                      __half *C,
                                                      int ldc);

    typedef cublasStatus_t (CUBLASWINAPI *CublasSgemm)(cublasHandle_t handle,
                                                      cublasOperation_t transa,
                                                      cublasOperation_t transb,
                                                      int m,
                                                      int n,
                                                      int k,
                                                      float *alpha, /* host or device pointer */
                                                      float *A,
                                                      int lda,
                                                      float *B,
                                                      int ldb,
                                                      float *beta, /* host or device pointer */
                                                      float *C,
                                                      int ldc);

    typedef cublasStatus_t (CUBLASWINAPI *CublasDgemm)(cublasHandle_t handle,
                                                      cublasOperation_t transa,
                                                      cublasOperation_t transb,
                                                      int m,
                                                      int n,
                                                      int k,
                                                      double *alpha, /* host or device pointer */
                                                      double *A,
                                                      int lda,
                                                      double *B,
                                                      int ldb,
                                                      double *beta, /* host or device pointer */
                                                      double *C,
                                                      int ldc);

    typedef cublasStatus_t (CUBLASWINAPI *CublasSgemmEx)(cublasHandle_t handle,
                                                      cublasOperation_t transa,
                                                      cublasOperation_t transb,
                                                      int m,
                                                      int n,
                                                      int k,
                                                      float *alpha, /* host or device pointer */
                                                      void *A,
                                                      cublasDataType_t Atype,
                                                      int lda,
                                                      void *B,
                                                      cublasDataType_t Btype,
                                                      int ldb,
                                                      float *beta, /* host or device pointer */
                                                      void *C,
                                                      cublasDataType_t Ctype,
                                                      int ldc);

    typedef cublasStatus_t (CUBLASWINAPI *CublasHgemmBatched)(cublasHandle_t handle,
                                                          cublasOperation_t transa,
                                                          cublasOperation_t transb,
                                                          int m,
                                                          int n,
                                                          int k,
                                                          __half *alpha,  /* host or device pointer */
                                                          __half *Aarray[],
                                                          int lda,
                                                          __half *Barray[],
                                                          int ldb,
                                                          __half *beta,   /* host or device pointer */
                                                          __half *Carray[],
                                                          int ldc,
                                                          int batchCount);

    typedef cublasStatus_t (CUBLASWINAPI *CublasSgemmBatched)(cublasHandle_t handle,
                                                          cublasOperation_t transa,
                                                          cublasOperation_t transb,
                                                          int m,
                                                          int n,
                                                          int k,
                                                          float *alpha,  /* host or device pointer */
                                                          float *Aarray[],
                                                          int lda,
                                                          float *Barray[],
                                                          int ldb,
                                                          float *beta,   /* host or device pointer */
                                                          float *Carray[],
                                                          int ldc,
                                                          int batchCount);

    typedef cublasStatus_t (CUBLASWINAPI *CublasDgemmBatched)(cublasHandle_t handle,
                                                          cublasOperation_t transa,
                                                          cublasOperation_t transb,
                                                          int m,
                                                          int n,
                                                          int k,
                                                          double *alpha,  /* host or device pointer */
                                                          double *Aarray[],
                                                          int lda,
                                                          double *Barray[],
                                                          int ldb,
                                                          double *beta,  /* host or device pointer */
                                                          double *Carray[],
                                                          int ldc,
                                                          int batchCount);

    typedef enum{
        CUSOLVER_STATUS_SUCCESS=0,
        CUSOLVER_STATUS_NOT_INITIALIZED=1,
        CUSOLVER_STATUS_ALLOC_FAILED=2,
        CUSOLVER_STATUS_INVALID_VALUE=3,
        CUSOLVER_STATUS_ARCH_MISMATCH=4,
        CUSOLVER_STATUS_MAPPING_ERROR=5,
        CUSOLVER_STATUS_EXECUTION_FAILED=6,
        CUSOLVER_STATUS_INTERNAL_ERROR=7,
        CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED=8,
        CUSOLVER_STATUS_NOT_SUPPORTED = 9,
        CUSOLVER_STATUS_ZERO_PIVOT=10,
        CUSOLVER_STATUS_INVALID_LICENSE=11
    } cusolverStatus_t;

    typedef enum {
        CUSOLVER_EIG_TYPE_1=1,
        CUSOLVER_EIG_TYPE_2=2,
        CUSOLVER_EIG_TYPE_3=3
    } cusolverEigType_t ;

    typedef enum {
        CUSOLVER_EIG_MODE_NOVECTOR=0,
        CUSOLVER_EIG_MODE_VECTOR=1
    } cusolverEigMode_t ;

    struct cusolverDnContext;
    typedef struct cusolverDnContext *cusolverDnHandle_t;

    typedef cusolverStatus_t (CUSOLVERAPI *CusolverDnSgesvdBufferSize)(
        cusolverDnHandle_t handle,
        int m,
        int n,
        int *lwork);

    typedef cusolverStatus_t (CUSOLVERAPI *CusolverDnDgesvdBufferSize)(
        cusolverDnHandle_t handle,
        int m,
        int n,
        int *lwork);

    typedef cusolverStatus_t (CUSOLVERAPI *CusolverDnSgesvd)(
        cusolverDnHandle_t handle,
        signed char jobu,
        signed char jobvt,
        int m,
        int n,
        float *A,
        int lda,
        float *S,
        float *U,
        int ldu,
        float *VT,
        int ldvt,
        float *work,
        int lwork,
        float *rwork,
        int  *info);

    typedef cusolverStatus_t (CUSOLVERAPI *CusolverDnDgesvd)(
        cusolverDnHandle_t handle,
        signed char jobu,
        signed char jobvt,
        int m,
        int n,
        double *A,
        int lda,
        double *S,
        double *U,
        int ldu,
        double *VT,
        int ldvt,
        double *work,
        int lwork,
        double *rwork,
        int *info);

    enum BlasFunctions {
        GEMV = 0,
        GEMM = 1,
    };

    class BlasHelper {
    private:
		bool _hasHgemv = false;
		bool _hasHgemm = false;
		bool _hasHgemmBatch = false;

        bool _hasSgemv = false;
        bool _hasSgemm = false;
        bool _hasSgemmBatch = false;

        bool _hasDgemv = false;
        bool _hasDgemm = false;
        bool _hasDgemmBatch = false;

        CblasSgemv cblasSgemv;
        CblasDgemv cblasDgemv;
        CblasSgemm cblasSgemm;
        CblasDgemm cblasDgemm;
        CblasSgemmBatch cblasSgemmBatch;
        CblasDgemmBatch cblasDgemmBatch;
        LapackeSgesvd lapackeSgesvd;
        LapackeDgesvd lapackeDgesvd;
        LapackeSgesdd lapackeSgesdd;
        LapackeDgesdd lapackeDgesdd;

        CublasSgemv cublasSgemv;
        CublasDgemv cublasDgemv;
        CublasHgemm cublasHgemm;
        CublasSgemm cublasSgemm;
        CublasDgemm cublasDgemm;
        CublasSgemmEx cublasSgemmEx;
        CublasHgemmBatched cublasHgemmBatched;
        CublasSgemmBatched cublasSgemmBatched;
        CublasDgemmBatched cublasDgemmBatched;
        CusolverDnSgesvdBufferSize cusolverDnSgesvdBufferSize;
        CusolverDnDgesvdBufferSize cusolverDnDgesvdBufferSize;
        CusolverDnSgesvd cusolverDnSgesvd;
        CusolverDnDgesvd cusolverDnDgesvd;

    public:
        static BlasHelper& getInstance();

        void initializeFunctions(Nd4jPointer *functions);
		void initializeDeviceFunctions(Nd4jPointer *functions);

        template <typename T>
        bool hasGEMV();

        template <typename T>
        bool hasGEMM();

        bool hasGEMM(const sd::DataType dtype);
        bool hasGEMV(const sd::DataType dtype);

        template <typename T>
        bool hasBatchedGEMM();

        CblasSgemv sgemv();
        CblasDgemv dgemv();

        CblasSgemm sgemm();
        CblasDgemm dgemm();

        CblasSgemmBatch sgemmBatched();
        CblasDgemmBatch dgemmBatched();

        LapackeSgesvd sgesvd();
        LapackeDgesvd dgesvd();

        LapackeSgesdd sgesdd();
        LapackeDgesdd dgesdd();

        // destructor
        ~BlasHelper() noexcept;
    };
}

#endif


text:
/*
 *  libDefs.h
 *  This file and its associated source file contain any functions and
 *   statements (within the actual program) that need to be modified
 *   when adding a new library. This includes any work done on complex
 *   numbers, FFT planning and execution, and would also contain the
 *   definition of the appropriate MPI type to represent the complex
 *   number, except that I decided to assume they were all just two
 *   MPI_DOUBLEs. (Bit compatability can be verified.)
 *
 * Defining HAS_AUTO in compilation is an instruction to use the parallel
 *  version of the library as well.
 *
 * We use the less eye-friendly version of the C complex number type
 *  to avoid namespace collisions -- complex is equivalent to _Complex.
 *
 *  Created by Ian Kirker on 14/04/2008.
 *
 */

#ifndef HEADER_LIBDEFS
#define HEADER_LIBDEFS

#include <complex.h>
#include <mpi.h>

#ifndef FFT_fftw2
#ifndef FFT_fftw3
#ifndef FFT_mkl
#ifndef FFT_essl
#ifndef FFT_acml
#define FFT_fftw3 /* Default to make test building easier. */
#endif
#endif
#endif
#endif
#endif

/* Include FFT library of choice */
#ifdef FFT_fftw2
	#ifdef FFTW_TYPE_SPECIFIED
		#include <dfftw.h>
		#ifdef HAS_AUTO
			#include <dfftw_mpi.h>
		#endif
	#else
		#include <fftw.h>
		#ifdef HAS_AUTO
			#include <fftw_mpi.h>
		#endif
	#endif

	#define FFT_NAME "fftw2"
	#define FFT_fftw2_LIBKEY
	typedef fftw_complex complexType;
	typedef fftw_plan oneDplanType;
	typedef fftwnd_plan twoDplanType;
	#ifdef HAS_AUTO
		typedef fftwnd_mpi_plan parallelPlanType;
	#endif
#endif

#ifdef FFT_fftw3
	#include <fftw3.h>
	#define FFT_NAME "fftw3"
	#define FFT_fftw3_LIBKEY 1

	typedef fftw_complex complexType;
	typedef fftw_plan planType;
	#ifdef HAS_AUTO
		typedef fftw_plan parallelPlanType;
	#endif
#endif

#ifdef FFT_essl
	#include <essl.h>
	#define FFT_NAME "essl"
	#define FFT_essl_LIBKEY 2

	typedef dcmplx complexType;
	typedef double planType[20000];

	#ifdef HAS_AUTO
		#include <pessl.h>
		#include <Cblacs.h>
		/* This refers to the BLACS handle you get back from the grid init calls. */
		typedef int parallelPlanType;
	#endif
#endif

#ifdef FFT_acml
	#include <acml.h>
	#define FFT_NAME "acml"
	#define FFT_acml_LIBKEY 3
	typedef doublecomplex complexType;

	/* In ACML, the closest you can get to making a plan is to run
	 * a transform and have it save the data in the array it uses as a buffer,
	 * so we alloc some space later but say that the planType is a pointer now.
	 */
	typedef complexType *planType;
#endif

#ifdef FFT_mkl
	#include <mkl_dfti.h>
	#define FFT_NAME "mkl"
	#define FFT_mkl_LIBKEY 4
	typedef _Complex double complexType;
	typedef DFTI_DESCRIPTOR_HANDLE planType;
	#ifdef HAS_AUTO
		typedef DFTI_DESCRIPTOR_DM_HANDLE parallelPlanType;
	#endif
#endif

void prepareFFTs(complexType *data, int decomp, int use2DFFT, int extent, int domainSize[2], MPI_Comm commColumn);
void performFFTset(complexType *data, complexType *buffer, int extent, int domainSize[2]);
void perform2DFFT(complexType *data, complexType *buffer, int extent, int domainSize[2]);
void performAutomatic3DFFT(complexType *data, complexType *buffer, int extent, int domainSize[2]);
void cleanUpFFTs(int decomp);
int libraryHasAutomaticDecomposition();

void printLib();
void complexSwap(complexType *, complexType *);
void complexSet(complexType *z, double real, double imag);
void complexAssign( complexType *z, complexType w );
double complexAbsNorm(complexType z, complexType w);
_Complex double complexNative( complexType z );

#endif


text:
/* randist/multinomial.c
 *
 * Copyright (C) 2002 Gavin E. Crooks <gec@compbio.berkeley.edu>
 *
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>

/* The multinomial distribution has the form

                                      N!           n_1  n_2      n_K
   prob(n_1, n_2, ... n_K) = -------------------- p_1  p_2  ... p_K
                             (n_1! n_2! ... n_K!)

   where n_1, n_2, ... n_K are nonnegative integers, sum_{k=1,K} n_k = N,
   and p = (p_1, p_2, ..., p_K) is a probability distribution.

   Random variates are generated using the conditional binomial method.
   This scales well with N and does not require a setup step.

   Ref:
   C.S. David, The computer generation of multinomial random variates,
   Comp. Stat. Data Anal. 16 (1993) 205-217
*/

void
gsl_ran_multinomial (const gsl_rng * r, const size_t K,
                     const unsigned int N, const double p[], unsigned int n[])
{
  size_t k;
  double norm = 0.0;
  double sum_p = 0.0;

  unsigned int sum_n = 0;

  /* p[k] may contain non-negative weights that do not sum to 1.0.
   * Even a probability distribution will not exactly sum to 1.0
   * due to rounding errors.
   */

  for (k = 0; k < K; k++)
    {
      norm += p[k];
    }

  for (k = 0; k < K; k++)
    {
      if (p[k] > 0.0)
        {
          n[k] = gsl_ran_binomial (r, p[k] / (norm - sum_p), N - sum_n);
        }
      else
        {
          n[k] = 0;
        }

      sum_p += p[k];
      sum_n += n[k];
    }

}

double
gsl_ran_multinomial_pdf (const size_t K,
                         const double p[], const unsigned int n[])
{
  return exp (gsl_ran_multinomial_lnpdf (K, p, n));
}

double
gsl_ran_multinomial_lnpdf (const size_t K,
                           const double p[], const unsigned int n[])
{
  size_t k;
  unsigned int N = 0;
  double log_pdf = 0.0;
  double norm = 0.0;

  for (k = 0; k < K; k++)
    {
      N += n[k];
    }

  for (k = 0; k < K; k++)
    {
      norm += p[k];
    }

  log_pdf = gsl_sf_lnfact (N);

  for (k = 0; k < K; k++)
    {
      log_pdf -= gsl_sf_lnfact (n[k]);
    }

  for (k = 0; k < K; k++)
    {
      log_pdf += log (p[k] / norm) * n[k];
    }

  return log_pdf;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified