text:
/*
Ballistic: a software to benchmark ballistic models.

AUTHORS: Javier Burguete Tolosa.

Copyright 2018, AUTHORS.

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

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

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

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

/**
 * \file ballistic.c
 * \brief Source file with the main function.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2018.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <libxml/parser.h>
#include <glib.h>
#include "config.h"
#include "utils.h"
#include "equation.h"
#include "method.h"
#include "runge-kutta.h"
#include "multi-steps.h"

#define DEBUG_BALLISTIC 0       ///< macro to debug the ballistic functions.

long double convergence_factor;
///< convergence factor.
unsigned int ntrajectories;
///< number of projectil trajectories to calculate.
unsigned int convergence;
///< number of convergence steps.

/**
 * Function to read the basic input data.
 *
 * \return 1 on success, 0 on error.
 */
static inline int
convergence_read_xml (xmlNode * node)      ///< XML node.
{
	const char *message[] = {
		"Bad trajectories number",
		"Bad convergence steps",
		"Bad convergence factor"
	};
	int e, error_code;

#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_read_xml: start\n");
#endif
  ntrajectories = xml_node_get_uint (node, XML_TRAJECTORIES, &error_code);
  if (error_code || !ntrajectories)
	  {
			e = 0;
      goto fail;
		}
  convergence = xml_node_get_uint (node, XML_CONVERGENCE, &error_code);
  if (error_code || !convergence)
	  {
			e = 1;
      goto fail;
		}
  convergence_factor = xml_node_get_float (node, XML_FACTOR, &error_code);
  if (error_code || convergence_factor <= 0.)
	  {
			e = 2;
      goto fail;
		}
#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_read_xml: success\n");
  fprintf (stderr, "convergence_read_xml: end\n");
#endif
  return 1;

fail:
	error_add (message[e]);
#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_read_xml: error\n");
  fprintf (stderr, "convergence_read_xml: end\n");
#endif
  return 0;
}

/**
 * Function to open a numerical method on a XML node.
 *
 * \return 0 on error, 1 on Runge-Kutta method, 2 on multi-steps method.
 */
static inline int
method_open_xml (MultiSteps * ms,
		             RungeKutta * rk,
								 xmlNode * node)
{
	char *message[] = {
    "No numerical method XML node",
    "Bad Runge-Kutta data",
    "Bad multi-steps data",
    "Unknown numerical method"
	};
	int e, m;

#if DEBUG_BALLISTIC
  fprintf (stderr, "method_open_xml: start\n");
#endif
  if (!node)
    {
      e = 0;
      goto fail;
    }
  if (!xmlStrcmp (node->name, XML_RUNGE_KUTTA))
    {
      if (!runge_kutta_read_xml (rk, node))
        {
          e = 1;
          goto fail;
        }
      runge_kutta_init_variables (rk);
			m = 1;
    }
  else if (!xmlStrcmp (node->name, XML_MULTI_STEPS))
    {
      if (!multi_steps_read_xml (ms, node))
        {
          e = 2;
          goto fail;
        }
      multi_steps_init_variables (ms);
			m = 2;
    }
  else
    {
      e = 3;
      goto fail;
    }
#if DEBUG_BALLISTIC
  fprintf (stderr, "method_open_xml: success\n");
  fprintf (stderr, "method_open_xml: end\n");
#endif
	return m;

fail:
	error_add (message[e]);
#if DEBUG_BALLISTIC
  fprintf (stderr, "method_open_xml: error\n");
  fprintf (stderr, "method_open_xml: end\n");
#endif
	return 0;
}

/**
 * Function to show the error message.
 */
static void
show_error ()
{
  printf ("ERROR!\n%s", error_message);
  g_free (error_message);
	error_message = NULL;
}

/**
 * Function to print the solution values.
 */
static void
print_solution (char *label,    ///< label.
                long double *r0,        ///< position vector.
                long double *r1)        ///< velocity vector.
{
  printf ("%s\n", label);
  printf ("x = %.19Le\n", r0[0]);
  printf ("y = %.19Le\n", r0[1]);
  printf ("z = %.19Le\n", r0[2]);
  printf ("vx = %.19Le\n", r1[0]);
  printf ("vy = %.19Le\n", r1[1]);
  printf ("vz = %.19Le\n", r1[2]);
}

/**
 * Function to print the numerical errors.
 */
static void
print_error (char *label,       ///< label.
             long double *r1,   ///< numerical solution vector.
             long double *r2)   ///< analytical solution vector.
{
  printf ("%s = %.19Le\n", label, distance (r1, r2));
}

/**
 * Function to perform a convergence analysis of a method.
 *
 * \return 0 on success, error code on error.
 */
static inline int
convergence_run (xmlNode * node,   ///< XML node.
                 char *output)  ///< results file name.
{
	const char *message[] = {
		NULL,
		"Bad convergence data",
		"No equation XML node",
		"Unknown numerical method",
		"Bad numerical method data"
	};
  MultiSteps ms[1];
  RungeKutta rk[1];
  Equation eq[1];
  Method *m;
  gsl_rng *rng;
  FILE *file;
  long double sr0[3], sr1[3];
  long double t, l0r0, l2r0, l0r1, l2r1, e;
	int er, me;
  unsigned int i, j;
#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_run: start\n");
#endif
	er = 0;
  if (!convergence_read_xml (node))
	  {
      er = 1;
			goto fail;
		}
#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_run: initing equation\n");
#endif
	node = node->children;
	if (!node)
	  {
			er = 2;
			goto fail;
		}
  if (!equation_read_xml (eq, node, 0))
	  {
			er = 3;
      goto fail;
		}
#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_run: initing method\n");
#endif
	node = node->next;
	me = method_open_xml (ms, rk, node);
	switch (me)
	  {
		case 1:
      m = RUNGE_KUTTA_METHOD (rk);
			break;
		case 2:
      m = MULTI_STEPS_METHOD (ms);
			break;
		default:
			er = 4;
			goto fail;
		}
  rng = gsl_rng_alloc (gsl_rng_taus2);
  file = fopen (output, "w");
  for (j = 0; j < convergence; ++j)
    {
      gsl_rng_set (rng, 0l);
      nevaluations = 0l;
      l0r0 = l2r0 = l0r1 = l2r1 = 0.L;
      for (i = 0; i < ntrajectories; ++i)
        {
#if DEBUG_BALLISTIC
          fprintf (stderr, "convergence_run: initing equation data\n");
#endif
          equation_init (eq, rng);
#if DEBUG_BALLISTIC
          fprintf (stderr, "convergence_run: initing variables\n");
#endif
          equation_solution (eq, r0, r1, 0.);
          equation_acceleration (eq, r0, r1, r2, 0.L);
#if DEBUG_BALLISTIC
          fprintf (stderr, "convergence_run: running\n");
#endif
          if (me == 1)
            t = runge_kutta_run (rk, eq);
          else
            t = multi_steps_run (ms, eq);
#if DEBUG_BALLISTIC
          fprintf (stderr, "convergence_run: solutions\n");
          print_solution ("Numerical solution", r0, r1);
          printf ("Time = %.19Le\n", t);
#endif
          switch (eq->land_type)
            {
            case 0:
              equation_solution (eq, sr0, sr1, eq->tf);
              break;
            default:
              t = equation_solve (eq, sr0, sr1);
            }
#if DEBUG_BALLISTIC
          print_solution ("Analytical solution", sr0, sr1);
          printf ("Time = %.19Le\n", t);
          print_error ("Position error", r0, sr0);
          print_error ("Velocity error", r1, sr1);
#endif
          e = distance (r0, sr0);
          l0r0 = fmaxl (l0r0, e);
          l2r0 += e * e;
          e = distance (r1, sr1);
          l0r1 = fmaxl (l0r1, e);
          l2r1 += e * e;

        }
#if DEBUG_BALLISTIC
      fprintf (stderr, "convergence_run: saving results\n");
#endif
      l2r0 = sqrtl (l2r0 / ntrajectories);
      l2r1 = sqrtl (l2r1 / ntrajectories);
      fprintf (file, "%lu %.19Le %.19Le %.19Le %.19Le %.19Le %.19Le\n",
               nevaluations, l0r0, l2r0, l0r1, l2r1, kt, m->emt);
      switch (eq->size_type)
        {
        case 0:
          dt *= convergence_factor;
          break;
        default:
          kt *= convergence_factor;
        }
      m->emt *= convergence_factor;
      if (me == 2)
        RUNGE_KUTTA_METHOD (MULTI_STEPS_RUNGE_KUTTA (ms))->emt
          *= convergence_factor;
    }
  fclose (file);
  printf ("Time = %.19Le\n", t);
#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_run: deleting method\n");
#endif
  if (me == 1)
    runge_kutta_delete (rk);
  else
    multi_steps_delete (ms);
  gsl_rng_free (rng);
fail:
	if (er)
	  {
			error_add (message[er]);
		  show_error ();
		}
#if DEBUG_BALLISTIC
  fprintf (stderr, "convergence_run: end\n");
#endif
  return 0;
}

/**
 * Function to calculate a ballistic trajectory.
 *
 * \return 0 on success, error code on error.
 */
static inline int
ballistic_run (xmlNode * node)    ///< XML node.
{
  const char *message[] = {
    NULL,
    "Unable to open the XML root element",
    "Bad XML file",
    "No equation XML node",
    "Bad equation data",
  };
  MultiSteps ms[1];
  RungeKutta rk[1];
  Equation eq[1];
  long double sr0[3], sr1[3];
  long double t;
  int e, m;
#if DEBUG_BALLISTIC
  fprintf (stderr, "ballistic_run: start\n");
#endif
  e = 0;
  node = node->children;
  if (!node)
    {
      e = 3;
      goto end;
    }
  if (!equation_read_xml (eq, node, 1))
    {
      e = 4;
      goto end;
    }
  node = node->next;
  if (!node)
    {
      e = 5;
      goto end;
    }
  if (!xmlStrcmp (node->name, XML_RUNGE_KUTTA))
    {
      if (!runge_kutta_read_xml (rk, node))
        {
          e = 6;
          goto end;
        }
      runge_kutta_init_variables (rk);
			m = 1;
    }
  else if (!xmlStrcmp (node->name, XML_MULTI_STEPS))
    {
      if (!multi_steps_read_xml (ms, node))
        {
          e = 7;
          goto end;
        }
      multi_steps_init_variables (ms);
			m = 2;
    }
  else
    {
      e = 8;
      goto end;
    }
  nevaluations = 0l;
#if DEBUG_BALLISTIC
  fprintf (stderr, "ballistic_run: initing variables\n");
#endif
  equation_solution (eq, r0, r1, 0.);
  equation_acceleration (eq, r0, r1, r2, 0.L);
#if DEBUG_BALLISTIC
  fprintf (stderr, "ballistic_run: running\n");
#endif
  switch (m)
	  {
		case 1:
			t = runge_kutta_run (rk, eq);
			break;
		default:
			t = multi_steps_run (ms, eq);
		}
#if DEBUG_BALLISTIC
  fprintf (stderr, "ballistic_run: solutions\n");
#endif
  print_solution ("Numerical solution", r0, r1);
  printf ("Time = %.19Le\n", t);
  switch (eq->land_type)
    {
    case 0:
      equation_solution (eq, sr0, sr1, eq->tf);
      break;
    default:
      t = equation_solve (eq, sr0, sr1);
    }
  print_solution ("Analytical solution", sr0, sr1);
  printf ("Time = %.19Le\n", t);
  print_error ("Position error", r0, sr0);
  print_error ("Velocity error", r1, sr1);
  printf ("Time = %.19Le\n", t);
#if DEBUG_BALLISTIC
  fprintf (stderr, "ballistic_run: deleting method\n");
#endif
  switch (m)
	  {
		case 1:
      runge_kutta_delete (rk);
			break;
		default:
      multi_steps_delete (ms);
		}
end:
#if DEBUG_BALLISTIC
  fprintf (stderr, "ballistic_run: end\n");
#endif
	if (e)
    error_add (message[e]);
  return e;
}

/**
 * Main function
 *
 * \return 0 on success, error code otherwise.
 */
int
main (int argn,                 ///< number of arguments.
      char **argc)              ///< array of argument chars.
{
	const char *message[] = {
		NULL,
    "The syntax is:\n./ballistic input_file output_file\n",
		"Unable to open the input file",
		"Bad XML root element",
	  "Bad ballistic run",
	  "Bad convergence run",
		"Unknown model"
	};
  xmlDoc *doc;
	xmlNode *node;
  int e;
#if DEBUG_BALLISTIC
  fprintf (stderr, "main: start\n");
#endif
	e = 0;
	if (argn != 3)
	  {
			e = 1;
			goto end;
		}
  xmlKeepBlanksDefault (0);
  doc = xmlParseFile (argc[1]);
  if (!doc)
	  {
			e = 2;
			goto end;
		}
  node = xmlDocGetRootElement (doc);
  if (!node)
    {
      e = 3;
      goto end;
    }
  if (!xmlStrcmp (node->name, XML_BALLISTIC))
	  {
      e = ballistic_run (node);
			if (e)
			  {
					e = 4;
					goto end;
				}
		}
	else if (!xmlStrcmp (node->name, XML_CONVERGENCE))
	  {
      e = convergence_run (node, argc[2]);
			if (e)
			  {
					e = 5;
					goto end;
				}
		}
	else
	  e = 6;
  xmlFreeDoc (doc);
end:
  if (e)
    {
			error_add (message[e]);
      show_error ();
    }
  return e;
}

ext: c

max_stars_repo_name: jburguete/ballistic

max_issues_repo_name: jburguete/ballistic

max_forks_repo_name: jburguete/ballistic


text:
/* ode-initval/odeiv2.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

/* Author:  G. Jungman */
/* Modified by Tuomo Keskitalo */

#ifndef __GSL_ODEIV2_H__
#define __GSL_ODEIV2_H__

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_types.h>

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

__BEGIN_DECLS
/* Description of a system of ODEs.
 *
 * y' = f(t,y) = dydt(t, y)
 *
 * The system is specified by giving the right-hand-side
 * of the equation and possibly a jacobian function.
 *
 * Some methods require the jacobian function, which calculates
 * the matrix dfdy and the vector dfdt. The matrix dfdy conforms
 * to the GSL standard, being a continuous range of floating point
 * values, in row-order.
 *
 * As with GSL function objects, user-supplied parameter
 * data is also present.
 */
  typedef struct
{
  int (*function) (double t, const double y[], double dydt[], void *params);
  int (*jacobian) (double t, const double y[], double *dfdy, double dfdt[],
                   void *params);
  size_t dimension;
  void *params;
}
gsl_odeiv2_system;

/* Function evaluation macros */

#define GSL_ODEIV_FN_EVAL(S,t,y,f)  (*((S)->function))(t,y,f,(S)->params)
#define GSL_ODEIV_JA_EVAL(S,t,y,dfdy,dfdt)  (*((S)->jacobian))(t,y,dfdy,dfdt,(S)->params)

/* Type definitions */

typedef struct gsl_odeiv2_step_struct gsl_odeiv2_step;
typedef struct gsl_odeiv2_control_struct gsl_odeiv2_control;
typedef struct gsl_odeiv2_evolve_struct gsl_odeiv2_evolve;
typedef struct gsl_odeiv2_driver_struct gsl_odeiv2_driver;

/* Stepper object
 *
 * Opaque object for stepping an ODE system from t to t+h.
 * In general the object has some state which facilitates
 * iterating the stepping operation.
 */

typedef struct
{
  const char *name;
  int can_use_dydt_in;
  int gives_exact_dydt_out;
  void *(*alloc) (size_t dim);
  int (*apply) (void *state, size_t dim, double t, double h, double y[],
                double yerr[], const double dydt_in[], double dydt_out[],
                const gsl_odeiv2_system * dydt);
  int (*set_driver) (void *state, const gsl_odeiv2_driver * d);
  int (*reset) (void *state, size_t dim);
  unsigned int (*order) (void *state);
  void (*free) (void *state);
}
gsl_odeiv2_step_type;

struct gsl_odeiv2_step_struct
{
  const gsl_odeiv2_step_type *type;
  size_t dimension;
  void *state;
};

/* Available stepper types */

GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk2;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk4;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rkf45;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rkck;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk8pd;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk2imp;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk4imp;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_bsimp;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk1imp;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_msadams;
GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_msbdf;

/* Stepper object methods */

gsl_odeiv2_step *gsl_odeiv2_step_alloc (const gsl_odeiv2_step_type * T,
                                        size_t dim);
int gsl_odeiv2_step_reset (gsl_odeiv2_step * s);
void gsl_odeiv2_step_free (gsl_odeiv2_step * s);
const char *gsl_odeiv2_step_name (const gsl_odeiv2_step * s);
unsigned int gsl_odeiv2_step_order (const gsl_odeiv2_step * s);
int gsl_odeiv2_step_apply (gsl_odeiv2_step * s, double t, double h,
                           double y[], double yerr[], const double dydt_in[],
                           double dydt_out[], const gsl_odeiv2_system * dydt);
int gsl_odeiv2_step_set_driver (gsl_odeiv2_step * s,
                                const gsl_odeiv2_driver * d);

/* Step size control object. */

typedef struct
{
  const char *name;
  void *(*alloc) (void);
  int (*init) (void *state, double eps_abs, double eps_rel, double a_y,
               double a_dydt);
  int (*hadjust) (void *state, size_t dim, unsigned int ord, const double y[],
                  const double yerr[], const double yp[], double *h);
  int (*errlevel) (void *state, const double y, const double dydt,
                   const double h, const size_t ind, double *errlev);
  int (*set_driver) (void *state, const gsl_odeiv2_driver * d);
  void (*free) (void *state);
}
gsl_odeiv2_control_type;

struct gsl_odeiv2_control_struct
{
  const gsl_odeiv2_control_type *type;
  void *state;
};

/* Possible return values for an hadjust() evolution method */

#define GSL_ODEIV_HADJ_INC   1  /* step was increased */
#define GSL_ODEIV_HADJ_NIL   0  /* step unchanged     */
#define GSL_ODEIV_HADJ_DEC (-1) /* step decreased     */

/* General step size control methods.
 *
 * The hadjust() method controls the adjustment of
 * step size given the result of a step and the error.
 * Valid hadjust() methods must return one of the codes below.
 * errlevel function calculates the desired error level D0.
 *
 * The general data can be used by specializations
 * to store state and control their heuristics.
 */

gsl_odeiv2_control *gsl_odeiv2_control_alloc (const gsl_odeiv2_control_type *
                                              T);
int gsl_odeiv2_control_init (gsl_odeiv2_control * c, double eps_abs,
                             double eps_rel, double a_y, double a_dydt);
void gsl_odeiv2_control_free (gsl_odeiv2_control * c);
int gsl_odeiv2_control_hadjust (gsl_odeiv2_control * c, gsl_odeiv2_step * s,
                                const double y[], const double yerr[],
                                const double dydt[], double *h);
const char *gsl_odeiv2_control_name (const gsl_odeiv2_control * c);
int gsl_odeiv2_control_errlevel (gsl_odeiv2_control * c, const double y,
                                 const double dydt, const double h,
                                 const size_t ind, double *errlev);
int gsl_odeiv2_control_set_driver (gsl_odeiv2_control * c,
                                   const gsl_odeiv2_driver * d);

/* Available control object constructors.
 *
 * The standard control object is a four parameter heuristic
 * defined as follows:
 *    D0 = eps_abs + eps_rel * (a_y |y| + a_dydt h |y'|)
 *    D1 = |yerr|
 *    q  = consistency order of method (q=4 for 4(5) embedded RK)
 *    S  = safety factor (0.9 say)
 *
 *                      /  (D0/D1)^(1/(q+1))  D0 >= D1
 *    h_NEW = S h_OLD * |
 *                      \  (D0/D1)^(1/q)      D0 < D1
 *
 * This encompasses all the standard error scaling methods.
 *
 * The y method is the standard method with a_y=1, a_dydt=0.
 * The yp method is the standard method with a_y=0, a_dydt=1.
 */

gsl_odeiv2_control *gsl_odeiv2_control_standard_new (double eps_abs,
                                                     double eps_rel,
                                                     double a_y,
                                                     double a_dydt);
gsl_odeiv2_control *gsl_odeiv2_control_y_new (double eps_abs, double eps_rel);
gsl_odeiv2_control *gsl_odeiv2_control_yp_new (double eps_abs,
                                               double eps_rel);

/* This controller computes errors using different absolute errors for
 * each component
 *
 *    D0 = eps_abs * scale_abs[i] + eps_rel * (a_y |y| + a_dydt h |y'|)
 */

gsl_odeiv2_control *gsl_odeiv2_control_scaled_new (double eps_abs,
                                                   double eps_rel, double a_y,
                                                   double a_dydt,
                                                   const double scale_abs[],
                                                   size_t dim);

/* Evolution object */

struct gsl_odeiv2_evolve_struct
{
  size_t dimension;
  double *y0;
  double *yerr;
  double *dydt_in;
  double *dydt_out;
  double last_step;
  unsigned long int count;
  unsigned long int failed_steps;
  const gsl_odeiv2_driver *driver;
};

/* Evolution object methods */

gsl_odeiv2_evolve *gsl_odeiv2_evolve_alloc (size_t dim);
int gsl_odeiv2_evolve_apply (gsl_odeiv2_evolve * e, gsl_odeiv2_control * con,
                             gsl_odeiv2_step * step,
                             const gsl_odeiv2_system * dydt, double *t,
                             double t1, double *h, double y[]);
int gsl_odeiv2_evolve_apply_fixed_step (gsl_odeiv2_evolve * e,
                                        gsl_odeiv2_control * con,
                                        gsl_odeiv2_step * step,
                                        const gsl_odeiv2_system * dydt,
                                        double *t, const double h0,
                                        double y[]);
int gsl_odeiv2_evolve_reset (gsl_odeiv2_evolve * e);
void gsl_odeiv2_evolve_free (gsl_odeiv2_evolve * e);
int gsl_odeiv2_evolve_set_driver (gsl_odeiv2_evolve * e,
                                  const gsl_odeiv2_driver * d);

/* Driver object
 *
 * This is a high level wrapper for step, control and
 * evolve objects.
 */

struct gsl_odeiv2_driver_struct
{
  const gsl_odeiv2_system *sys; /* ODE system */
  gsl_odeiv2_step *s;           /* stepper object */
  gsl_odeiv2_control *c;        /* control object */
  gsl_odeiv2_evolve *e;         /* evolve object */
  double h;                     /* step size */
  double hmin;                  /* minimum step size allowed */
  double hmax;                  /* maximum step size allowed */
  unsigned long int n;          /* number of steps taken */
  unsigned long int nmax;       /* Maximum number of steps allowed */
};

/* Driver object methods */

gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system *
                                                  sys,
                                                  const gsl_odeiv2_step_type *
                                                  T, const double hstart,
                                                  const double epsabs,
                                                  const double epsrel);
gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system *
                                                   sys,
                                                   const gsl_odeiv2_step_type
                                                   * T, const double hstart,
                                                   const double epsabs,
                                                   const double epsrel);
gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system
                                                       * sys,
                                                       const
                                                       gsl_odeiv2_step_type *
                                                       T, const double hstart,
                                                       const double epsabs,
                                                       const double epsrel,
                                                       const double a_y,
                                                       const double a_dydt,
                                                       const double
                                                       scale_abs[]);
gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_standard_new (const
                                                         gsl_odeiv2_system *
                                                         sys,
                                                         const
                                                         gsl_odeiv2_step_type
                                                         * T,
                                                         const double hstart,
                                                         const double epsabs,
                                                         const double epsrel,
                                                         const double a_y,
                                                         const double a_dydt);
int gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin);
int gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax);
int gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d,
                                const unsigned long int nmax);
int gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double *t,
                             const double t1, double y[]);
int gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double *t,
                                        const double h,
                                        const unsigned long int n,
                                        double y[]);
int gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d);
void gsl_odeiv2_driver_free (gsl_odeiv2_driver * state);

__END_DECLS
#endif /* __GSL_ODEIV2_H__ */

ext: h

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
/*
** compute histogram density function
**
** G.Lohmann, Nov 2012
*/

#include <viaio/Vlib.h>

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>

#define SQR(x) ((x)*(x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

/* kernel density estimation */
double GaussKernel(double x)
{
  return exp(-0.5*x*x)/sqrt(2.0*M_PI);
}

double VKernelDensity(double x,double h,double nh,gsl_histogram *histogram)
{
  size_t i;
  double sum,upper,lower;

  sum = 0;
  for (i=0; i<gsl_histogram_bins(histogram); i++) {
    double nx = gsl_histogram_get (histogram,i);
    gsl_histogram_get_range (histogram,i,&lower,&upper);
    double xi = lower + (upper-lower)*0.5;
    double u = (x-xi)/h;
    sum += nx * GaussKernel(u);
  }
  return sum/nh;
}

void VPrintHistogram(gsl_histogram *histogram,int numperm,VString filename)
{
  size_t i;
  double z=0,mx=0,lower=0,upper=0,fz=0;
  extern double GaussKernel(double x);
  extern double VKernelDensity(double,double,double,gsl_histogram *);

  gsl_histogram_get_range (histogram,0,&lower,&upper);
  double alpha = 1.06;  /* silverman's rule */
  double sig = gsl_histogram_sigma(histogram);
  double kx = gsl_histogram_sum(histogram);
  double h  = alpha * sig * pow(kx,-0.2);
  double nh = h*kx;

  FILE *fph = fopen(filename,"w");
  if (!fph) VError(" err opening hist file");
  for (i=0; i<gsl_histogram_bins(histogram); i++) {
    gsl_histogram_get_range (histogram,i,&lower,&upper);
    z  = lower;
    fz = VKernelDensity(z,h,nh,histogram);
    mx = gsl_histogram_get (histogram,i);
    fprintf(fph,"%lf %lf %.2lf\n",lower,fz,mx);
  }
  fclose(fph);
}

/* update histogram */
void HistoUpdate(float *A,size_t nvox,gsl_histogram *hist)
{
  double tiny = 1.0e-8;
  size_t i,j;
  double hmin = gsl_histogram_min(hist);
  double hmax = gsl_histogram_max(hist);

  for (i=0; i<nvox; i++) {
    for (j=0; j<i; j++) {
      const size_t k=j+i*(i+1)/2;
      double u = (double)A[k];
      if (ABS(u) < tiny) continue;
      if (u < hmin) u = hmin;
      if (u > hmax-tiny) u = hmax-tiny;
      gsl_histogram_increment (hist,u);
    }
  }
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/**
 *
 * @file qwrapper_dgemm.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:56 2014
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dgemm(Quark *quark, Quark_Task_Flags *task_flags,
                      PLASMA_enum transA, int transB,
                      int m, int n, int k, int nb,
                      double alpha, const double *A, int lda,
                                                const double *B, int ldb,
                      double beta, double *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_dgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(double),         &alpha,     VALUE,
        sizeof(double)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(double)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(double),         &beta,      VALUE,
        sizeof(double)*nb*nb,    C,                 INOUT,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dgemm2( Quark *quark, Quark_Task_Flags *task_flags,
                        PLASMA_enum transA, int transB,
                        int m, int n, int k, int nb,
                        double alpha, const double *A, int lda,
                        const double *B, int ldb,
                        double beta, double *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_dgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(double),         &alpha,     VALUE,
        sizeof(double)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(double)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(double),         &beta,      VALUE,
        sizeof(double)*nb*nb,    C,                 INOUT | LOCALITY | GATHERV,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgemm_quark = PCORE_dgemm_quark
#define CORE_dgemm_quark PCORE_dgemm_quark
#endif
void CORE_dgemm_quark(Quark *quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int m;
    int n;
    int k;
    double alpha;
    double *A;
    int lda;
    double *B;
    int ldb;
    double beta;
    double *C;
    int ldc;

    quark_unpack_args_13(quark, transA, transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
    cblas_dgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        m, n, k,
        (alpha), A, lda,
        B, ldb,
        (beta), C, ldc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dgemm_f2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         double alpha, const double *A, int lda,
                                                   const double *B, int ldb,
                         double beta, double *C, int ldc,
                         double *fake1, int szefake1, int flag1,
                         double *fake2, int szefake2, int flag2)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_dgemm_f2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(double),         &alpha,     VALUE,
        sizeof(double)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(double)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(double),         &beta,      VALUE,
        sizeof(double)*nb*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(double)*szefake1, fake1,             flag1,
        sizeof(double)*szefake2, fake2,             flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgemm_f2_quark = PCORE_dgemm_f2_quark
#define CORE_dgemm_f2_quark PCORE_dgemm_f2_quark
#endif
void CORE_dgemm_f2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    double alpha;
    double *A;
    int LDA;
    double *B;
    int LDB;
    double beta;
    double *C;
    int LDC;
    void *fake1, *fake2;

    quark_unpack_args_15(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1, fake2);
    cblas_dgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        B, LDB,
        (beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dgemm_p2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         double alpha, const double *A, int lda,
                         const double **B, int ldb,
                         double beta, double *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_dgemm_p2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(double),         &alpha,     VALUE,
        sizeof(double)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(double*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(double),         &beta,      VALUE,
        sizeof(double)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgemm_p2_quark = PCORE_dgemm_p2_quark
#define CORE_dgemm_p2_quark PCORE_dgemm_p2_quark
#endif
void CORE_dgemm_p2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    double alpha;
    double *A;
    int LDA;
    double **B;
    int LDB;
    double beta;
    double *C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_dgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        *B, LDB,
        (beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dgemm_p3(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           double alpha, const double *A, int lda,
                           const double *B, int ldb,
                           double beta, double **C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_dgemm_p3_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(double),         &alpha,     VALUE,
        sizeof(double)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(double)*ldb*nb,   B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(double),         &beta,      VALUE,
        sizeof(double*),         C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgemm_p3_quark = PCORE_dgemm_p3_quark
#define CORE_dgemm_p3_quark PCORE_dgemm_p3_quark
#endif
void CORE_dgemm_p3_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    double alpha;
    double *A;
    int LDA;
    double *B;
    int LDB;
    double beta;
    double **C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_dgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        B, LDB,
        (beta), *C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dgemm_p2f1(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           double alpha, const double *A, int lda,
                           const double **B, int ldb,
                           double beta, double *C, int ldc,
                           double *fake1, int szefake1, int flag1)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_dgemm_p2f1_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(double),         &alpha,     VALUE,
        sizeof(double)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(double*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(double),         &beta,      VALUE,
        sizeof(double)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(double)*szefake1, fake1,             flag1,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgemm_p2f1_quark = PCORE_dgemm_p2f1_quark
#define CORE_dgemm_p2f1_quark PCORE_dgemm_p2f1_quark
#endif
void CORE_dgemm_p2f1_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    double alpha;
    double *A;
    int LDA;
    double **B;
    int LDB;
    double beta;
    double *C;
    int LDC;
    void *fake1;

    quark_unpack_args_14(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1);
    cblas_dgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        (alpha), A, LDA,
        *B, LDB,
        (beta), C, LDC);
}


text:
/**
 *
 * @file testing_cher2k.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:45:19 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_cmain.h"

static int check_solution(PLASMA_enum uplo, PLASMA_enum trans, int N, int K,
                          PLASMA_Complex32_t alpha, PLASMA_Complex32_t *A, int LDA,
                          PLASMA_Complex32_t *B, int LDB,
                          float beta,  PLASMA_Complex32_t *Cref, PLASMA_Complex32_t *Cplasma, int LDC);

int testing_cher2k(int argc, char **argv)
{
    /* Check for number of arguments*/
    if ( argc != 7 ){
        USAGE("HER2K", "alpha beta M N LDA LDB LDC",
              "   - alpha : alpha coefficient\n"
              "   - beta : beta coefficient\n"
              "   - N : number of columns and rows of matrix C and number of row of matrix A and B\n"
              "   - K : number of columns of matrix A and B\n"
              "   - LDA : leading dimension of matrix A\n"
              "   - LDB : leading dimension of matrix B\n"
              "   - LDC : leading dimension of matrix C\n");
        return -1;
    }

    PLASMA_Complex32_t alpha = (PLASMA_Complex32_t) atol(argv[0]);
    float beta  = (float) atol(argv[1]);
    int N     = atoi(argv[2]);
    int K     = atoi(argv[3]);
    int LDA   = atoi(argv[4]);
    int LDB   = atoi(argv[5]);
    int LDC   = atoi(argv[6]);
    int NKmax = max(N, K);

    float eps;
    int info_solution;
    int u, t;
    size_t LDAxK = LDA*NKmax;
    size_t LDBxK = LDB*NKmax;
    size_t LDCxN = LDC*N;

    PLASMA_Complex32_t *A      = (PLASMA_Complex32_t *)malloc(LDAxK*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *B      = (PLASMA_Complex32_t *)malloc(LDBxK*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *C      = (PLASMA_Complex32_t *)malloc(LDCxN*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *Cinit  = (PLASMA_Complex32_t *)malloc(LDCxN*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *Cfinal = (PLASMA_Complex32_t *)malloc(LDCxN*sizeof(PLASMA_Complex32_t));

    /* Check if unable to allocate memory */
    if ( (!A) || (!B) || (!Cinit) || (!Cfinal) ){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = LAPACKE_slamch_work('e');

    printf("\n");
    printf("------ TESTS FOR PLASMA CHER2K ROUTINE -------  \n");
    printf("            Size of the Matrix C %d by %d\n", N, K);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 10.\n");

    /*----------------------------------------------------------
    *  TESTING CHER2K
    */

    /* Initialize A,B */
    LAPACKE_clarnv_work(IONE, ISEED, LDAxK, A);
    LAPACKE_clarnv_work(IONE, ISEED, LDBxK, B);

    /* Initialize C */
    PLASMA_cplghe( (float)0., N, C, LDC, 51 );

    for (u=0; u<2; u++) {
        for (t=0; t<3; t++) {
            if (trans[t] == PlasmaTrans) continue;

            memcpy(Cinit,  C, LDCxN*sizeof(PLASMA_Complex32_t));
            memcpy(Cfinal, C, LDCxN*sizeof(PLASMA_Complex32_t));

            /* PLASMA CHER2K */
            PLASMA_cher2k(uplo[u], trans[t], N, K, alpha, A, LDA, B, LDB, beta, Cfinal, LDC);

            /* Check the solution */
            info_solution = check_solution(uplo[u], trans[t], N, K,
                                           alpha, A, LDA, B, LDB, beta, Cinit, Cfinal, LDC);

            if (info_solution == 0) {
                printf("***************************************************\n");
                printf(" ---- TESTING CHER2K (%5s, %s) ........... PASSED !\n", uplostr[u], transstr[t]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING CHER2K (%5s, %s) ... FAILED !\n", uplostr[u], transstr[t]);
                printf("************************************************\n");
            }
        }
    }

    free(A); free(B); free(C);
    free(Cinit); free(Cfinal);

    return 0;
}

/*--------------------------------------------------------------
 * Check the solution
 */

static int check_solution(PLASMA_enum uplo, PLASMA_enum trans, int N, int K,
                          PLASMA_Complex32_t alpha, PLASMA_Complex32_t *A, int LDA,
                          PLASMA_Complex32_t *B, int LDB,
                          float beta, PLASMA_Complex32_t *Cref, PLASMA_Complex32_t *Cplasma, int LDC)
{
    int info_solution;
    float Anorm, Bnorm, Cinitnorm, Cplasmanorm, Clapacknorm, Rnorm, result;
    float eps;
    PLASMA_Complex32_t beta_const;

    float *work = (float *)malloc(max(N, K)* sizeof(float));

    beta_const  = -1.0;
    Anorm       = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm),
                                      (trans == PlasmaNoTrans) ? N : K,
                                      (trans == PlasmaNoTrans) ? K : N, A, LDA, work);
    Bnorm       = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm),
                                      (trans == PlasmaNoTrans) ? N : K,
                                      (trans == PlasmaNoTrans) ? K : N, B, LDB, work);
    Cinitnorm   = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Cref,    LDC, work);
    Cplasmanorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Cplasma, LDC, work);

    cblas_cher2k(CblasColMajor, (CBLAS_UPLO)uplo, (CBLAS_TRANSPOSE)trans,
                 N, K, CBLAS_SADDR(alpha), A, LDA, B, LDB, (beta), Cref, LDC);

    Clapacknorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Cref, LDC, work);

    cblas_caxpy(LDC*N, CBLAS_SADDR(beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Cref, LDC, work);

    eps = LAPACKE_slamch_work('e');

    printf("Rnorm %e, Anorm %e, Cinitnorm %e, Cplasmanorm %e, Clapacknorm %e\n",
           Rnorm, Anorm, Cinitnorm, Cplasmanorm, Clapacknorm);

    result = Rnorm / ((Anorm + Bnorm + Cinitnorm) * N * eps);
    printf("============\n");
    printf("Checking the norm of the difference against reference CHER2K \n");
    printf("-- ||Cplasma - Clapack||_oo/((||A||_oo+||C||_oo).N.eps) = %e \n", result);

    if (  isnan(Rnorm) || isinf(Rnorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
         printf("-- The solution is suspicious ! \n");
         info_solution = 1;
    }
    else {
         printf("-- The solution is CORRECT ! \n");
         info_solution= 0 ;
    }

    free(work);

    return info_solution;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"

ccl_f2d_t *ccl_f2d_t_copy(ccl_f2d_t *f2d_o, int *status)
{
  int s2dstatus=0;
  ccl_f2d_t *f2d = malloc(sizeof(ccl_f2d_t));
  if (f2d == NULL)
    *status = CCL_ERROR_MEMORY;

  if(*status==0) {
    f2d->lkmin = f2d_o->lkmin;
    f2d->lkmax = f2d_o->lkmax;
    f2d->amin = f2d_o->amin;
    f2d->amax = f2d_o->amax;
    f2d->is_factorizable = f2d_o->is_factorizable;
    f2d->is_k_constant = f2d_o->is_k_constant;
    f2d->is_a_constant = f2d_o->is_a_constant;
    f2d->extrap_linear_growth = f2d_o->extrap_linear_growth;
    f2d->extrap_order_lok = f2d_o->extrap_order_lok;
    f2d->extrap_order_hik = f2d_o->extrap_order_hik;
    f2d->is_log = f2d_o->is_log;
    f2d->growth = f2d_o->growth;
    f2d->growth_factor_0 = f2d_o->growth_factor_0;
    f2d->growth_exponent = f2d_o->growth_exponent;

    if(f2d_o->fk != NULL) {
      f2d->fk = gsl_spline_alloc(gsl_interp_cspline,
                                 f2d_o->fk->size);
      if(f2d->fk == NULL)
        *status = CCL_ERROR_MEMORY;

      if(*status==0) {
        s2dstatus |= gsl_spline_init(f2d->fk, f2d_o->fk->x,
                                     f2d_o->fk->y, f2d_o->fk->size);
        if(s2dstatus)
          *status = CCL_ERROR_SPLINE;
      }
    }
    else
      f2d->fk = NULL;
  }

  if(*status==0) {
    if(f2d_o->fa != NULL) {
      f2d->fa = gsl_spline_alloc(gsl_interp_cspline,
                                 f2d_o->fa->size);
      if(f2d->fa == NULL)
        *status = CCL_ERROR_MEMORY;

      if(*status==0) {
        s2dstatus |= gsl_spline_init(f2d->fa, f2d_o->fa->x,
                                     f2d_o->fa->y, f2d_o->fa->size);
        if(s2dstatus)
          *status = CCL_ERROR_SPLINE;
      }
    }
    else
      f2d->fa = NULL;
  }

  if(*status==0) {
    if(f2d_o->fka != NULL) {
      f2d->fka = gsl_spline2d_alloc(gsl_interp2d_bicubic,
                                    f2d_o->fka->interp_object.xsize,
                                    f2d_o->fka->interp_object.ysize);
      if(f2d->fka == NULL)
        *status = CCL_ERROR_MEMORY;

      if(*status==0) {
        s2dstatus |= gsl_spline2d_init(f2d->fka, f2d_o->fka->xarr,
                                       f2d_o->fka->yarr, f2d_o->fka->zarr,
                                       f2d_o->fka->interp_object.xsize,
                                       f2d_o->fka->interp_object.ysize);
        if(s2dstatus)
          *status = CCL_ERROR_SPLINE;
      }
    }
    else
      f2d->fka = NULL;
  }

  return f2d;
}

ccl_f2d_t *ccl_f2d_t_new(int na,double *a_arr,
                         int nk,double *lk_arr,
                         double *fka_arr,
                         double *fk_arr,
                         double *fa_arr,
                         int is_factorizable,
                         int extrap_order_lok,
                         int extrap_order_hik,
                         ccl_f2d_extrap_growth_t extrap_linear_growth,
                         int is_fka_log,
                         double (*growth)(double),
                         double growth_factor_0,
                         int growth_exponent,
                         ccl_f2d_interp_t interp_type,
                         int *status) {
  int s2dstatus=0;
  ccl_f2d_t *f2d = malloc(sizeof(ccl_f2d_t));
  if (f2d == NULL)
    *status = CCL_ERROR_MEMORY;

  if (*status == 0) {
    is_factorizable = is_factorizable || (a_arr == NULL) || (lk_arr == NULL) || (fka_arr == NULL);
    f2d->is_factorizable = is_factorizable;
    f2d->is_k_constant = ((lk_arr == NULL) || ((fka_arr == NULL) && (fk_arr == NULL)));
    f2d->is_a_constant = ((a_arr == NULL) || ((fka_arr == NULL) && (fa_arr == NULL)));
    f2d->extrap_order_lok = extrap_order_lok;
    f2d->extrap_order_hik = extrap_order_hik;
    f2d->extrap_linear_growth = extrap_linear_growth;
    f2d->is_log = is_fka_log;
    f2d->growth = growth;
    f2d->growth_factor_0 = growth_factor_0;
    f2d->growth_exponent = growth_exponent;
    f2d->fka = NULL;
    f2d->fk = NULL;
    f2d->fa = NULL;

    if (!(f2d->is_k_constant)) { //If it's not constant
      f2d->lkmin = lk_arr[0];
      f2d->lkmax = lk_arr[nk-1];
    }
    if (!(f2d->is_a_constant)) {
      f2d->amin = a_arr[0];
      f2d->amax = a_arr[na-1];
    }
  }

  if ((extrap_order_lok > 2) || (extrap_order_lok < 0) || (extrap_order_hik > 2) || (extrap_order_hik < 0))
    *status = CCL_ERROR_INCONSISTENT;

  if ((extrap_linear_growth != ccl_f2d_cclgrowth) &&
      (extrap_linear_growth != ccl_f2d_customgrowth) &&
      (extrap_linear_growth != ccl_f2d_constantgrowth) &&
      (extrap_linear_growth != ccl_f2d_no_extrapol))
    *status = CCL_ERROR_INCONSISTENT;

  if(*status == 0) {
    switch(interp_type) {
    case(ccl_f2d_3):
      if (f2d->is_factorizable) {
        // Do not allocate spline if constant
        if(f2d->is_k_constant)
          f2d->fk = NULL;
        else { //Otherwise allocate and check
          f2d->fk = gsl_spline_alloc(gsl_interp_cspline, nk);
          if(f2d->fk == NULL)
            *status = CCL_ERROR_MEMORY;
        }

        // Do not allocate spline if constant
        if (f2d->is_a_constant)
          f2d->fa = NULL;
        else { //Otherwise allocate and check
          f2d->fa = gsl_spline_alloc(gsl_interp_cspline, na);
          if (f2d->fa == NULL)
            *status = CCL_ERROR_MEMORY;
        }
      }
      else {
        // Do not allocate spline if constant
        if ((f2d->is_k_constant) || (f2d->is_a_constant))
          f2d->fka = NULL;
        else { //Otherwise allocate and check
          f2d->fka = gsl_spline2d_alloc(gsl_interp2d_bicubic, nk, na);
          if (f2d->fka == NULL)
            *status = CCL_ERROR_MEMORY;
        }
      }
      break;

    default:
      f2d->fk = NULL;
      f2d->fa = NULL;
      f2d->fka = NULL;
    }
  }

  if (*status == 0) {
    if (f2d->is_factorizable) {
      if (f2d->fk != NULL)
        s2dstatus |= gsl_spline_init(f2d->fk, lk_arr, fk_arr, nk);
      if (f2d->fa != NULL)
        s2dstatus |= gsl_spline_init(f2d->fa, a_arr, fa_arr, na);
    }
    else {
      if (f2d->fka != NULL)
        s2dstatus=gsl_spline2d_init(f2d->fka, lk_arr, a_arr, fka_arr, nk, na);
    }
    if (s2dstatus)
      *status = CCL_ERROR_SPLINE;
  }

  return f2d;
}

double ccl_f2d_t_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, int *status) {
  int is_hiz, is_loz;
  double a_ev = a;
  if (f2d->is_a_constant) {
    is_hiz = 0;
    is_loz = 0;
  }
  else {
    is_hiz = a < f2d->amin;
    is_loz = a > f2d->amax;
    if (is_loz) { // Are we above the interpolation range in a?
      if (f2d->extrap_linear_growth == ccl_f2d_no_extrapol) {
        *status=CCL_ERROR_SPLINE_EV;
        return NAN;
      }
      a_ev = f2d->amax;
    }
    else if (is_hiz) { // Are we below the interpolation range in a?
      if (f2d->extrap_linear_growth == ccl_f2d_no_extrapol) {
        *status=CCL_ERROR_SPLINE_EV;
        return NAN;
      }
      a_ev = f2d->amin;
    }
  }

  int is_hik, is_lok;
  double fka_pre, fka_post;
  double lk_ev = lk;
  if (f2d->is_k_constant) {
    is_hik = 0;
    is_lok = 0;
  }
  else {
    is_hik = lk > f2d->lkmax;
    is_lok = lk < f2d->lkmin;
    if (is_hik) // Are we above the interpolation range in k?
      lk_ev = f2d->lkmax;
    else if (is_lok) // Are we below the interpolation range in k?
      lk_ev = f2d->lkmin;
  }

  // Evaluate spline
  int spstatus=0;
  if (f2d->is_factorizable) {
    double fk, fa;
    if (f2d->fk == NULL) {
      if (f2d->is_log)
        fk = 0;
      else
        fk = 1;
    }
    else
      spstatus |= gsl_spline_eval_e(f2d->fk, lk_ev, NULL, &fk);

    if (f2d->fa == NULL) {
      if (f2d->is_log)
        fa = 0;
      else
        fa = 1;
    }
    else
      spstatus |= gsl_spline_eval_e(f2d->fa, a_ev, NULL, &fa);
      if (f2d->is_log)
        fka_pre = fk+fa;
      else
        fka_pre = fk*fa;
  }
  else {
    if (f2d->fka == NULL) {
      if (f2d->is_log)
        fka_pre = 0;
      else
        fka_pre = 1;
    }
    else
      spstatus = gsl_spline2d_eval_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &fka_pre);
  }

  if (spstatus) {
    *status = CCL_ERROR_SPLINE_EV;
    return NAN;
  }

  // Now extrapolate in k if needed
  if (is_hik) {
    fka_post = fka_pre;
    if (f2d->extrap_order_hik > 0) {
      double pd;
      double dlk = lk-lk_ev;
      if (f2d->is_factorizable)
        spstatus = gsl_spline_eval_deriv_e(f2d->fk, lk_ev, NULL, &pd);
      else
        spstatus = gsl_spline2d_eval_deriv_x_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pd);
      if (spstatus) {
        *status = CCL_ERROR_SPLINE_EV;
        return NAN;
      }
      fka_post += pd*dlk;
      if (f2d->extrap_order_hik > 1) {
        if (f2d->is_factorizable)
          spstatus = gsl_spline_eval_deriv2_e(f2d->fk, lk_ev, NULL, &pd);
        else
          spstatus = gsl_spline2d_eval_deriv_xx_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pd);
        if (spstatus) {
          *status=CCL_ERROR_SPLINE_EV;
          return NAN;
        }
        fka_post += pd*dlk*dlk*0.5;
      }
    }
  }
  else if (is_lok) {
    fka_post = fka_pre;
    if (f2d->extrap_order_lok > 0) {
      double pd;
      double dlk = lk-lk_ev;
      if (f2d->is_factorizable)
        spstatus = gsl_spline_eval_deriv_e(f2d->fk, lk_ev, NULL, &pd);
      else
        spstatus = gsl_spline2d_eval_deriv_x_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pd);
      if (spstatus) {
        *status = CCL_ERROR_SPLINE_EV;
        return NAN;
      }
      fka_post += pd*dlk;

      if (f2d->extrap_order_lok > 1) {
        if (f2d->is_factorizable)
          spstatus = gsl_spline_eval_deriv2_e(f2d->fk, lk_ev, NULL, &pd);
        else
          spstatus = gsl_spline2d_eval_deriv_xx_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pd);
        if (spstatus) {
          *status = CCL_ERROR_SPLINE_EV;
          return NAN;
        }
        fka_post += pd*dlk*dlk*0.5;
      }
    }
  }
  else
    fka_post = fka_pre;

  // Exponentiate if needed
  if (f2d->is_log)
    fka_post = exp(fka_post);

  // Extrapolate in a if needed
  if (is_hiz) {
    double gz;
    if (f2d->extrap_linear_growth == ccl_f2d_cclgrowth) { // Use CCL's growth function
      ccl_cosmology *csm = (ccl_cosmology *)cosmo;
      if (!csm->computed_growth) {
        *status = CCL_ERROR_GROWTH_INIT;
        ccl_cosmology_set_status_message(
          csm,
          "ccl_f2d.c: ccl_f2d_t_eval(): growth factor splines have not been precomputed!");
        return NAN;
      }
      gz = (
        ccl_growth_factor(csm, a, status) /
        ccl_growth_factor(csm, a_ev, status));
    }
    else if (f2d->extrap_linear_growth == ccl_f2d_customgrowth) // Use internal growth function
      gz = f2d->growth(a) / f2d->growth(a_ev);
    else // Use constant growth factor
      gz = f2d->growth_factor_0;

    fka_post *= pow(gz, f2d->growth_exponent);
  }

  return fka_post;
}

void ccl_f2d_t_free(ccl_f2d_t *f2d)
{
  if(f2d != NULL) {
    if(f2d->fka != NULL)
      gsl_spline2d_free(f2d->fka);
    if(f2d->fk != NULL)
      gsl_spline_free(f2d->fk);
    if(f2d->fa != NULL)
      gsl_spline_free(f2d->fa);
    free(f2d);
  }
}

max_stars_repo_name: borisbolliet/CCL

max_issues_repo_name: borisbolliet/CCL

max_forks_repo_name: borisbolliet/CCL


text:
/* specfunc/result.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_result.h>

int
gsl_sf_result_smash_e(const gsl_sf_result_e10 * re, gsl_sf_result * r)
{
  if(re->e10 == 0) {
    /* nothing to smash */
    r->val = re->val;
    r->err = re->err;
    return GSL_SUCCESS;
  }
  else {
    const double av = fabs(re->val);
    const double ae = fabs(re->err);

    if(   GSL_SQRT_DBL_MIN < av && av < GSL_SQRT_DBL_MAX
       && GSL_SQRT_DBL_MIN < ae && ae < GSL_SQRT_DBL_MAX
       && 0.49*GSL_LOG_DBL_MIN  < re->e10 && re->e10 < 0.49*GSL_LOG_DBL_MAX
       ) {
      const double scale = exp(re->e10 * M_LN10);
      r->val = re->val * scale;
      r->err = re->err * scale;
      return GSL_SUCCESS;
    }
    else {
      return gsl_sf_exp_mult_err_e(re->e10*M_LN10, 0.0, re->val, re->err, r);
    }
  }
/*
  int stat_v;
  int stat_e;

  if(re->val == 0.0) {
    r->val = 0.0;
    stat_v = GSL_SUCCESS;
  }
  else {
    gsl_sf_result r_val;
    const double s = GSL_SIGN(re->val);
    const double x_v = re->e10*M_LN10 + log(fabs(re->val));
    stat_v = gsl_sf_exp_e(x_v, &r_val);
    r->val = s * r_val.val;
  }

  if(re->err == 0.0) {
    r->err = 0.0;
    stat_e = GSL_SUCCESS;
  }
  else if(re->val != 0.0) {
    r->err = fabs(r->val * re->err/re->val);
    stat_e = GSL_SUCCESS;
  }
  else {
    gsl_sf_result r_err;
    const double x_e = re->e10*M_LN10 + log(fabs(re->err));
    stat_e = gsl_sf_exp_e(x_e, &r_err);
    r->err = r_err.val;
  }

  return GSL_ERROR_SELECT_2(stat_v, stat_e);
*/
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef TRANS_TYPE
#define TRANS_TYPE

#include <gsl/gsl_vector.h>

// corresponding python variable: element of 'list_n_original_diff'
typedef struct {
    gsl_vector *inc1;
    gsl_vector *inc2;
    gsl_vector *skp1;
    gsl_vector *skp2;
    int inclu_len;
    int skip_len;
    int flag;
    char* id;
} odiff;

typedef struct {
    int batch_size;
    void **datum;
} batch_datum;

// used to represent the reture value of
// MLE_marginal_iteration and MLE_marginal_iteration_constrain.
// corresponding python variable: '[current_sum,[psi1,psi2,beta_0,beta_1,var1,var2]]'
typedef struct {
    double sum;
    struct {
        gsl_vector *psi1;
        gsl_vector *psi2;
        double beta0;
        double beta1;
        double var1;
        double var2;
    } params;
} mle_result;

// we use a linked list to represent a python list.
// corresponding python variable: 'list_n_original_diff'
typedef struct dnode{
    odiff* data;
    struct dnode* next;
    struct dnode* end;
} diff_list_node;

typedef struct node {
    void* data;
    struct node* next;
    struct node* end;
} list_node;

// we use fortran l_bfgs_b routine to solve the numerical optimization problem.
// (i.e. f = 0.0 at the optimal solution.)
// only the minimal intersection of C's and Fortran's many data types can be relied on:
// the following are some essential mapping between c type and fortran type
// counterpart of f2c.h
typedef int integer;
typedef float real;
typedef double doublereal;
typedef long int logical;

#endif

max_stars_repo_name: chunjie-sam-liu/rmats-turbo

max_issues_repo_name: chunjie-sam-liu/rmats-turbo

max_forks_repo_name: chunjie-sam-liu/rmats-turbo


text:
#include <sstream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>

// pretty print
std::string ppMatrix(const gsl_matrix *m) {
    std::ostringstream str;
    int r, c;

    str << "";
    for (r = 0; r < (int) m->size1; r++) {
        for (c = 0; c < (int) m->size2; c++) {
            str << gsl_matrix_get(m, r, c);
            if (c < m->size2 - 1) {
                str << " ";
            }
        }
        str << "\n";
    }
    return str.str();
}

void
matrixInverseCheck() {
	// Define the dimension n of the matrix
	// and the signum s (for LU decomposition)
	int n = 2;
	int s;

	// Define all the used matrices
	gsl_matrix * m = gsl_matrix_alloc (n, n);
	gsl_matrix * inverse = gsl_matrix_alloc (n, n);
	gsl_permutation * perm = gsl_permutation_alloc (n);

	// Fill the matrix m
        gsl_matrix_set (m, 0, 0, 0.378589);
        gsl_matrix_set (m, 0, 1, 0.971711);
        gsl_matrix_set (m, 0, 2, 0.016087);
        gsl_matrix_set (m, 0, 3, 0.037668);
        gsl_matrix_set (m, 0, 4, 0.312398);

        gsl_matrix_set (m, 1, 0, 0.756377);
        gsl_matrix_set (m, 1, 1, 0.345708);
        gsl_matrix_set (m, 1, 2, 0.922947);
        gsl_matrix_set (m, 1, 3, 0.846671);
        gsl_matrix_set (m, 1, 4, 0.856103);

        gsl_matrix_set (m, 2, 0, 0.732510);
        gsl_matrix_set (m, 2, 1, 0.108942);
        gsl_matrix_set (m, 2, 2, 0.476969);
        gsl_matrix_set (m, 2, 3, 0.398254);
        gsl_matrix_set (m, 2, 4, 0.507045);

        gsl_matrix_set (m, 3, 0, 0.162608);
        gsl_matrix_set (m, 3, 1, 0.227770);
        gsl_matrix_set (m, 3, 2, 0.533074);
        gsl_matrix_set (m, 3, 3, 0.807075);
        gsl_matrix_set (m, 3, 4, 0.180335);

        gsl_matrix_set (m, 4, 0, 0.517006);
        gsl_matrix_set (m, 4, 1, 0.315992);
        gsl_matrix_set (m, 4, 2, 0.914848);
        gsl_matrix_set (m, 4, 3, 0.460825);
        gsl_matrix_set (m, 4, 4, 0.731980);

	// Make LU decomposition of matrix m
	gsl_linalg_LU_decomp (m, perm, &s);

	// Invert the matrix m
	gsl_linalg_LU_invert (m, perm, inverse);
}

void
matrixOpsCheck() {
    fprintf(stdout, "\n Matrix Operations ..............................................\n");
    int L = 3, M = 3;
    int i, j;
    gsl_matrix *m = gsl_matrix_alloc(L, M);

    for (i = 0; i < L; i++)
        for (j = 0; j < M; j++)
            gsl_matrix_set(m, i, j, 0.23 + 100 * i + j);

    fprintf(stdout, "m = %s\n", ppMatrix(m).c_str());

    gsl_matrix_free(m);

    fprintf(stdout, "\n");
}

int
main(void) {
    matrixOpsCheck();

    return 0;
}

max_stars_repo_name: anthonyohare/SwiftScientificLibrary

max_issues_repo_name: anthonyohare/SwiftScientificLibrary

max_forks_repo_name: anthonyohare/SwiftScientificLibrary


text:
/*
Copyright (c) 2011-2013, UT-Battelle, LLC
All rights reserved

[PsimagLite, Version 1.0.0]
[by G.A., Oak Ridge National Laboratory]

UT Battelle Open Source Software License 11242008

OPEN SOURCE LICENSE

Subject to the conditions of this License, each
contributor to this software hereby grants, free of
charge, to any person obtaining a copy of this software
and associated documentation files (the "Software"), a
perpetual, worldwide, non-exclusive, no-charge,
royalty-free, irrevocable copyright license to use, copy,
modify, merge, publish, distribute, and/or sublicense
copies of the Software.

1. Redistributions of Software must retain the above
copyright and license notices, this list of conditions,
and the following disclaimer.  Changes or modifications
to, or derivative works of, the Software should be noted
with comments and the contributor and organization's
name.

2. Neither the names of UT-Battelle, LLC or the
Department of Energy nor the names of the Software
contributors may be used to endorse or promote products
derived from this software without specific prior written
permission of UT-Battelle.

3. The software and the end-user documentation included
with the redistribution, with or without modification,
must include the following acknowledgment:

"This product includes software produced by UT-Battelle,
LLC under Contract No. DE-AC05-00OR22725  with the
Department of Energy."

*********************************************************
DISCLAIMER

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

NEITHER THE UNITED STATES GOVERNMENT, NOR THE UNITED
STATES DEPARTMENT OF ENERGY, NOR THE COPYRIGHT OWNER, NOR
ANY OF THEIR EMPLOYEES, REPRESENTS THAT THE USE OF ANY
INFORMATION, DATA, APPARATUS, PRODUCT, OR PROCESS
DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.

*********************************************************

*/
/** \ingroup PsimagLite */
/*@{*/

/*! \file GslWrapper.h
 *
 *  Wrapper for GSL functions and types
 */

#ifndef GSL_WRAPPER_H_
#define GSL_WRAPPER_H_

#ifdef USE_GSL
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_expint.h>
#endif

#include <stdexcept>
#include "Vector.h"

namespace PsimagLite {

#ifndef USE_GSL
class GslWrapper {
public:
	typedef int DummyType;
	typedef DummyType gsl_integration_workspace;
	typedef double (* GslWrapperFunctionType) (double, void * );
	typedef struct {
		double val;
		double err;
	} gsl_sf_result;

	typedef void (* gsl_error_handler_t) (const char *,const char *,int,int);

	struct gsl_function {
		GslWrapperFunctionType function;
		void * params;
	};

	gsl_error_handler_t  gsl_set_error_handler (gsl_error_handler_t) const
	{

		thereSnoGsl();
		gsl_error_handler_t *fx = new gsl_error_handler_t();

		return *fx;
	}

	gsl_integration_workspace * gsl_integration_workspace_alloc (SizeType) const
	{
		thereSnoGsl();
		int* x = new int;
		return x;
	}

	void gsl_integration_workspace_free (gsl_integration_workspace*) const
	{
		thereSnoGsl();
	}

	int gsl_integration_qag(const gsl_function*,
	                        double,
	                        double,
	                        double,
	                        double,
	                        size_t,
	                        int,
	                        gsl_integration_workspace*,
	                        double*,
	                        double*) const
	{
		thereSnoGsl();
		return 0;
	}

	int gsl_integration_qagiu(gsl_function*,
	                          double,
	                          double,
	                          double,
	                          size_t ,
	                          gsl_integration_workspace*,
	                          double*,
	                          double*) const
	{
		thereSnoGsl();
		return 0;
	}

	int gsl_integration_qagp (const gsl_function*,
	                          double*,
	                          SizeType,
	                          double,
	                          double,
	                          SizeType,
	                          gsl_integration_workspace*,
	                          double*,
	                          double*) const
	{
		thereSnoGsl();
		return 0;
	}

	int gsl_integration_qagi(gsl_function*,
	                         double,
	                         double,
	                         size_t,
	                         gsl_integration_workspace*,
	                         double*,
	                         double*) const
	{
		thereSnoGsl();
		return 0;
	}

	void printError(int) const
	{
		thereSnoGsl();
	}

	int gsl_sf_lngamma_complex_e(double,
	                             double,
	                             gsl_sf_result*,
	                             gsl_sf_result*) const
	{
		thereSnoGsl();
		return 0;
	}

	int gsl_sf_Ci_e(double, gsl_sf_result*)
	{
		thereSnoGsl();
		return 0;
	}

private:

	void thereSnoGsl() const
	{
		throw RuntimeError("You need to compile with the GSL\n");
	}

}; // class GslWrapper

#else

class GslWrapper {
public:

	typedef ::gsl_integration_workspace gsl_integration_workspace;
	typedef ::gsl_function gsl_function;
	typedef ::gsl_sf_result gsl_sf_result;

	void printError(int status) const
	{
		std::cerr<<"GslWrapper: error: "<<gsl_strerror(status)<<"\n";
	}

	gsl_error_handler_t * gsl_set_error_handler (gsl_error_handler_t * new_handler) const
	{
		return ::gsl_set_error_handler(new_handler);
	}

	gsl_integration_workspace * gsl_integration_workspace_alloc (SizeType n) const
	{
		return ::gsl_integration_workspace_alloc(n);
	}

	void gsl_integration_workspace_free (gsl_integration_workspace * w) const
	{
		return ::gsl_integration_workspace_free(w);
	}

	int gsl_integration_qagi(gsl_function* f,
	                         double epsabs,
	                         double epsrel,
	                         size_t limit,
	                         gsl_integration_workspace* workspace,
	                         double* result,
	                         double* abserr) const
	{
		return ::gsl_integration_qagi(f,epsabs,epsrel,limit,workspace,result,abserr);
	}

	int gsl_integration_qagiu(gsl_function* f,
	                          double a,
	                          double epsabs,
	                          double epsrel,
	                          size_t limit,
	                          gsl_integration_workspace* workspace,
	                          double* result,
	                          double* abserr) const
	{
		return ::gsl_integration_qagiu(f,a,epsabs,epsrel,limit,workspace,result,abserr);
	}

	int gsl_integration_qagp (const gsl_function * f,
	                          double * pts,
	                          SizeType npts,
	                          double epsabs,
	                          double epsrel,
	                          SizeType limit,
	                          gsl_integration_workspace * workspace,
	                          double * result,
	                          double * abserr) const
	{
		return ::gsl_integration_qagp(f,
		                              pts,
		                              npts,
		                              epsabs,
		                              epsrel,
		                              limit,
		                              workspace,
		                              result,
		                              abserr);
	}

	int gsl_integration_qag(const gsl_function * f,
	                        double a,
	                        double b,
	                        double epsabs,
	                        double epsrel,
	                        size_t limit,
	                        int key,
	                        gsl_integration_workspace* workspace,
	                        double* result,
	                        double* abserr) const
	{
		return ::gsl_integration_qag(f,
		                             a,
		                             b,
		                             epsabs,
		                             epsrel,
		                             limit,
		                             key,
		                             workspace,
		                             result,
		                             abserr);
	}

	int gsl_sf_lngamma_complex_e(double zr,
	                             double zi,
	                             gsl_sf_result* lnr,
	                             gsl_sf_result* arg) const
	{
		return ::gsl_sf_lngamma_complex_e(zr, zi, lnr, arg);
	}

	int gsl_sf_Ci_e(double x, gsl_sf_result* result) const
	{
		return ::gsl_sf_Ci_e(x,result);
	}

}; // class GslWrapper

#endif

} // namespace PsimagLite

/*@}*/
#endif // GSL_WRAPPER_H_