text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_odeiv2.h>

#include <stdlib.h>
#include <string.h>
#include <stdbool.h>

#include <time.h>
#include <math.h>

#include <GL/glut.h>
#include <GL/gl.h>

#include <unistd.h>
#include <signal.h>
#include <sys/wait.h>
#include <sys/types.h>

#define SLICE	360

typedef struct _line_vertex line_vertex;
struct _line_vertex
{
	float x1;
	float y1;

	float x2;
	float y2;
};

float text_x = 64.0f, text_y = 58.0f;
const float TRANS_VAL = 100.0f;

uint32_t mouse_loc;
int x, y;

bool clicked = false;
int pos_x, pos_y;
int mov_x, mov_y;
int click_x, click_y;
int special_key;

bool run = false;

line_vertex vertices[9] = {
	{ 50.0f, 50.0f, 50.0f, 30.0f },
	{ 50.0f, 50.0f, 55.0f, 50.0f },
	{ 75.0f, 50.0f, 80.0f, 50.0f },
	{ 80.0f, 50.0f, 80.0f, 30.0f },
	{ 55.0f, 55.0f, 55.0f, 45.0f },
	{ 55.0f, 55.0f, 75.0f, 55.0f },
	{ 75.0f, 55.0f, 75.0f, 45.0f },
	{ 55.0f, 45.0f, 75.0f, 45.0f }
};

line_vertex pick_vertices[9] = {
	{ 45.0f, 25.0f, 45.0f, 65.0f },
	{ 45.0f, 25.0f, 85.0f, 25.0f },
	{ 45.0f, 65.0f, 85.0f, 65.0f },
	{ 85.0f, 25.0f, 85.0f, 65.0f }
};

line_vertex power_vertices[5] = {
	{ 40.0f, 50.0f, 40.0f, 60.0f },
	{ 50.0f, 50.0f, 30.0f, 50.0f },
	{ 45.0f, 45.0f, 35.0f, 45.0f },
	{ 40.0f, 45.0f, 40.0f, 35.0f }
};

line_vertex ground_vertices[5] = {
	{ 40.0f, 50.0f, 40.0f, 60.0f },
	{ 50.0f, 50.0f, 30.0f, 50.0f },
	{ 47.5f, 47.5f, 32.5f, 47.5f },
	{ 45.0f, 45.0f, 35.0f, 45.0f }
};

int *jac;

int rhs (double t, const double y[], double f[], void *params_ptr);
//int jacobian (double t, const double y[], double *dfdy, double dfdt[], void *params_ptr);

void calc_differential_eq(void);
void draw_differential_eq_plot(void);
void draw_lattice(void);
void draw_lattice2(void);
void draw_volt_plot(void);

void reshape(int w, int h);

typedef struct _xy_plot xy_plot;
struct _xy_plot
{
	double x;
	double y;
};

//xy_plot plot_res[1001];
xy_plot plot_res[100000001];
//double volt_res[100000001];

float resistance;
float inductance;
float voltage;

int sub_window;
int sub_window2;

void sim_reshape(int w, int h)
{
    GLfloat n_range = 100.0f;

    if(h == 0)
        h = 1;

    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

	// glOrtho(표기하고자 하는 최소 x값, 최대 x값, 화면상 나타는 최소 y, 최대 y, z는 관계없음)
	// 즉 내가 보고자 하는 x, y, z의 범위를 지정할 수 있음
	glOrtho(-0.000005, 0.000005, -0.7, 0.7, -1, 1);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
}

void simulate(void)
{
    glClearColor(0.0, 0.0, 0.0, 1.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glLoadIdentity();

    glColor3f(1, 0, 0);

    glBegin(GL_LINE_LOOP);
    	glVertex3f(100.0, 0.0, 0.0);
        glVertex3f(-100.0, 0.0, 0.0);
    glEnd();

    glColor3f(0.0, 1.0, 0.0);

    glBegin(GL_LINE_LOOP);
        glVertex3f(0.0, 100.0, 0.0);
        glVertex3f(0.0, -100.0, 0.0);
    glEnd();

    draw_differential_eq_plot();
    draw_lattice();

    glutSwapBuffers();
}

void sim_reshape2(int w, int h)
{
    GLfloat n_range = 100.0f;

    if(h == 0)
        h = 1;

    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

	glOrtho(-0.000005, 0.000005, -2000000, 500000, -1, 1);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
}

void simulate2(void)
{
    glClearColor(0.0, 0.0, 0.0, 1.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glLoadIdentity();

    glColor3f(1, 0, 0);

    glBegin(GL_LINE_LOOP);
    	glVertex3f(100.0, 0.0, 0.0);
        glVertex3f(-100.0, 0.0, 0.0);
    glEnd();

    glColor3f(0.0, 1.0, 0.0);

    glBegin(GL_LINE_LOOP);
        glVertex3f(0.0, 500000.0, 0.0);
        glVertex3f(0.0, -2000000.0, 0.0);
    glEnd();

    draw_volt_plot();
    draw_lattice2();

    glutSwapBuffers();
}

void keyboard_handler(unsigned char key, int x, int y)
{
	int pid;
	int i, status;

	// fork -> execve
	// OpenGL -> GPU
#if 0
	char *argv[5] = {"print_simulation"};
	char *envp[] = {0};
	char volt[64] = {0};
	char coil[64] = {0};
	char resist[64] = {0};
#endif

	switch(key)
	{
		case 's':
			printf("Simulation Preparation\n");
			printf("저항값 입력: ");
			scanf("%f", &resistance);
			printf("인덕턴스값 입력: ");
			scanf("%f", &inductance);
			printf("입력 전원값 입력: ");
			scanf("%f", &voltage);

			printf("R = %f, L = %f, V = %f\n", resistance, inductance, voltage);

			calc_differential_eq();
			run = true;
			break;
		case 'r':
			printf("Run Simulation\n");

#if 0
			pid = fork();

			if (pid > 0)
			{
				printf("wait for simulation finished\n");
				wait(&status);
			}
			else if (pid == 0)
			{
				execve("./print_simulation", argv, envp);
			}
			else
			{
				perror("fork() ");
				exit(-1);
			}
#endif

			sub_window = glutCreateWindow("Current i(t)");

			glutDisplayFunc(simulate);
			glutReshapeFunc(sim_reshape);

			sub_window2 = glutCreateWindow("Voltage V(t)");

			glutDisplayFunc(simulate2);
			glutReshapeFunc(sim_reshape2);
			break;
		case 't':
			printf("Termination\n");
			glutDestroyWindow(sub_window);
			glutDestroyWindow(sub_window2);
			break;
		case 27:
			exit(1);
			break;
	}
}

void on_mouse(int button, int state, int x, int y)
{
	special_key = glutGetModifiers();

	if (button == GLUT_LEFT_BUTTON & state == GLUT_DOWN)
	{
		if (clicked == false)
		{
			printf("Click - Down\n");
			click_x = x;
			click_y = y;
			clicked = true;

			printf("click_x = %d, click_y = %d\n", click_x, click_y);
		}
	}

	if (button == GLUT_LEFT_BUTTON & state == GLUT_UP)
	{
		printf("Click - Up\n");
		clicked = false;
	}

	int window_width = glutGet(GLUT_WINDOW_WIDTH);
	int window_height = glutGet(GLUT_WINDOW_HEIGHT);

	GLbyte color[4];
	GLfloat depth;
	GLuint index;

	glReadPixels(x, window_height - y - 1, 1, 1, GL_RGBA, GL_UNSIGNED_BYTE, color);
	glReadPixels(x, window_height - y - 1, 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &depth);
	glReadPixels(x, window_height - y - 1, 1, 1, GL_STENCIL_INDEX, GL_UNSIGNED_INT, &index);

	printf("Clicked on pixel %d, %d, color %02hhx%02hhx%02hhx%02hhx, depth %f, stencil index %u\n",
			x, y, color[0], color[1], color[2], color[3], depth, index);

}

void drag_mouse(int x, int y)
{
	pos_x = x;
	pos_y = y;

	special_key = glutGetModifiers();

	if (clicked == true)
	{
		printf("clicked!\n");

		mov_x = (click_x - pos_x) / 3;
		mov_y = (click_y - pos_y) / 2.5;

		click_x = pos_x;
		click_y = pos_y;

		printf("mov_x = %d, mov_y = %d\n", mov_x, mov_y);

		int i;

		for (i = 0; i < 8; i++)
		{
			vertices[i].x1 -= mov_x;
			vertices[i].x2 -= mov_x;
			vertices[i].y1 += mov_y;
			vertices[i].y2 += mov_y;
		}

		for (i = 0; i < 4; i++)
		{
			pick_vertices[i].x1 -= mov_x;
			pick_vertices[i].x2 -= mov_x;
			pick_vertices[i].y1 += mov_y;
			pick_vertices[i].y2 += mov_y;
		}

		for (i = 0; i < 4; i++)
		{
			power_vertices[i].x1 -= mov_x;
			power_vertices[i].x2 -= mov_x;
			power_vertices[i].y1 += mov_y;
			power_vertices[i].y2 += mov_y;
		}

		for (i = 0; i < 4; i++)
		{
			ground_vertices[i].x1 -= mov_x;
			ground_vertices[i].x2 -= mov_x;
			ground_vertices[i].y1 += mov_y;
			ground_vertices[i].y2 += mov_y;
		}

		text_x -= mov_x;
		text_y += mov_y;

		glutPostRedisplay();
	}
}

void drawString (char *s)
{
	unsigned int i;

	for (i = 0; i < strlen (s); i++)
		glutBitmapCharacter (GLUT_BITMAP_HELVETICA_10, s[i]);
}

void drawStringBig (char *s)
{
	unsigned int i;

	for (i = 0; i < strlen (s); i++)
		glutBitmapCharacter (GLUT_BITMAP_HELVETICA_18, s[i]);
}

void draw_outline(void)
{
	int i;

	glColor3f(0, 1, 0);

	glBegin(GL_LINES);
		for (i = 0; i < 8; i++)
		{
			glVertex2f(pick_vertices[i].x1, pick_vertices[i].y1);
			glVertex2f(pick_vertices[i].x2, pick_vertices[i].y2);
		}
	glEnd();
}

void draw_inductance_outline(void)
{
	int i;

	glColor3f(0, 1, 0);

	glBegin(GL_LINES);
		for (i = 0; i < 4; i++)
		{
			glVertex2f(pick_vertices[i].x1 + 50.0f, pick_vertices[i].y1);
			glVertex2f(pick_vertices[i].x2 + 50.0f, pick_vertices[i].y2);
		}
	glEnd();
}

void draw_resistance(void)
{
	int i;
	static char label[100];

	glColor3f(1, 0, 1);

	glBegin(GL_LINES);
		for (i = 0; i < 8; i++)
		{
			glVertex2f(vertices[i].x1, vertices[i].y1);
			glVertex2f(vertices[i].x2, vertices[i].y2);
		}
	glEnd();

	sprintf (label, "R");
	glRasterPos2f (text_x, text_y);
	drawStringBig (label);
}

void draw_inductance(void)
{
	int i;
	static char label[100];

	glColor3f(1, 0, 1);

	glBegin(GL_LINES);
		for (i = 0; i < 8; i++)
		{
			glVertex2f(vertices[i].x1 + 50.0f, vertices[i].y1);
			glVertex2f(vertices[i].x2 + 50.0f, vertices[i].y2);
		}
	glEnd();

	sprintf (label, "L");
	glRasterPos2f (text_x + 50.0f, text_y);
	drawStringBig (label);
}

void draw_power_outline(void)
{
	int i;

	glColor3f(0, 1, 0);

	glBegin(GL_LINES);
		for (i = 0; i < 4; i++)
		{
			glVertex2f(pick_vertices[i].x1 - 75.0f, pick_vertices[i].y1 + 20.0f);
			glVertex2f(pick_vertices[i].x2 - 75.0f, pick_vertices[i].y2 + 20.0f);
		}
	glEnd();
}

void draw_power(void)
{
	int i;
	static char label[100];

	glColor3f(1, 0, 1);

	glBegin(GL_LINES);
		for (i = 0; i < 4; i++)
		{
			glVertex2f(power_vertices[i].x1 - 50.0f, power_vertices[i].y1 + 15.0f);
			glVertex2f(power_vertices[i].x2 - 50.0f, power_vertices[i].y2 + 15.0f);
		}
	glEnd();

	sprintf (label, "V");
	glRasterPos2f (text_x - 70.0f, text_y + 20.0f);
	drawStringBig (label);
}

void draw_ground_outline(void)
{
	int i;

	glColor3f(0, 1, 0);

	glBegin(GL_LINES);
		for (i = 0; i < 4; i++)
		{
			glVertex2f(pick_vertices[i].x1 - 75.0f, pick_vertices[i].y1 - 20.0f);
			glVertex2f(pick_vertices[i].x2 - 75.0f, pick_vertices[i].y2 - 20.0f);
		}
	glEnd();
}

void draw_ground(void)
{
	int i;
	static char label[100];

	glColor3f(1, 0, 1);

	glBegin(GL_LINES);
		for (i = 0; i < 4; i++)
		{
			glVertex2f(ground_vertices[i].x1 - 50.0f, ground_vertices[i].y1 - 30.0f);
			glVertex2f(ground_vertices[i].x2 - 50.0f, ground_vertices[i].y2 - 30.0f);
		}
	glEnd();

	sprintf (label, "GND");
	glRasterPos2f (text_x - 70.0f, text_y - 30.0f);
	drawStringBig (label);
}

void display(void)
{
	glClearColor(0.0, 0.0, 0.0, 1.0);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();

	glColor3f(1, 0, 0);

	glBegin(GL_LINE_LOOP);
	glVertex3f(100.0, 0.0, 0.0);
        glVertex3f(-100.0, 0.0, 0.0);
        glEnd();

        glColor3f(0.0, 1.0, 0.0);

        glBegin(GL_LINE_LOOP);
        glVertex3f(0.0, 100.0, 0.0);
        glVertex3f(0.0, -100.0, 0.0);
    glEnd();

	draw_resistance();
	draw_outline();

	draw_inductance();
	draw_inductance_outline();

	draw_power();
	draw_power_outline();

	draw_ground();
	draw_ground_outline();

#if 0
	if (run)
	{
		draw_differential_eq_plot();
	}

	draw_lattice();
#endif

	glutSwapBuffers();
}

void reshape(int w, int h)
{
    GLfloat n_range = 100.0f;

    if(h == 0)
        h = 1;

    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

#if 0
    if(w <= h)
        glOrtho(-n_range, n_range, -n_range * h / w, n_range * h / w, -n_range, n_range);
    else
        glOrtho(-n_range * w / h, n_range * w / h, -n_range, n_range, -n_range, n_range);
#endif

	//glOrtho(-10, 10, -0.5, 0.5, -1, 1);
	glOrtho(-100, 100, -100, 100, 1.0, -1.0);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
}

int calc_vertices_num(void)
{
	double tmin = 0.0, tmax = 10.0, delta_t = 0.01;
	return (tmax - tmin) / delta_t;
}

void calc_differential_eq(void)
{
    int dim = 1;

	gsl_odeiv2_system sys = {rhs, NULL, dim, NULL};
	gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);

	int cnt = 0;

	double i;
    double x0 = 0.0,  xf = 10.0;
    double x = x0;
    //double y[1] = { 0.5 };
    double y[1] = { 0.0 };
	double delta = 0.0000001;

	for (i = 0; i <= xf; i += delta)
    {
        double xi = x0 + i * (xf-x0) / xf;
        int status = gsl_odeiv2_driver_apply (d, &x, xi, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error, return value=%d\n", status);
            break;
        }

		plot_res[cnt].x = x;
		plot_res[cnt].y = y[0];

		//volt_res[cnt] = y[0] * (resistance / inductance);

		if (cnt < 1001)
		{
       		printf ("%.8e %.8e\n", plot_res[cnt].x, plot_res[cnt].y);
		}

		cnt++;
    }

	printf("Finish Calc\n");

    gsl_odeiv2_driver_free (d);

#if 0
    const gsl_odeiv_step_type *type_ptr = gsl_odeiv_step_rkf45;

    gsl_odeiv_step *step_ptr = gsl_odeiv_step_alloc (type_ptr, dimension);
    gsl_odeiv_control *control_ptr = gsl_odeiv_control_y_new (eps_abs, eps_rel);
    gsl_odeiv_evolve *evolve_ptr = gsl_odeiv_evolve_alloc (dimension);

    gsl_odeiv_system my_system;	/* structure with the rhs function, etc. */

    double mu = 10;		/* parameter for the diffeq */
    double y[2];			/* current solution vector */

    double t, t_next;		/* current and next independent variable */
    double tmin, tmax, delta_t;	/* range of t and step size for output */

    double h = 1e-6;		/* starting step size for ode solver */

	int cnt = 0;
	int vertex_num = calc_vertices_num();

    my_system.function = rhs;	/* the right-hand-side functions dy[i]/dt */
    my_system.jacobian = NULL;
    my_system.dimension = dimension;	/* number of diffeq's */
    my_system.params = NULL;

    tmin = 0.;			/* starting t value */
    tmax = 10.;			/* final t value */
    delta_t = 0.01;

    y[1] = 0.5;			/* initial x value */

    t = tmin;             /* initialize t */

    for (t_next = tmin + delta_t; t_next <= tmax; t_next += delta_t)
    {
        while (t < t_next)	/* evolve from t to t_next */
        {
            gsl_odeiv_evolve_apply (evolve_ptr, control_ptr, step_ptr,
                                    &my_system, &t, t_next, &h, y);
        }

		plot_res[cnt].x = t;
		plot_res[cnt].y = y[0];
        printf ("%.5e %.5e %.5e\n", t, plot_res[cnt].x, plot_res[cnt].y); /* print at t=t_next */

		cnt++;
    }

    gsl_odeiv_evolve_free (evolve_ptr);
    gsl_odeiv_control_free (control_ptr);
    gsl_odeiv_step_free (step_ptr);
#endif
}

void draw_lattice(void)
{
	int i;
	float lattice = 0.000001;
	char buf[64] = {0};

	glColor3f(1, 1, 0);
	glBegin(GL_LINES);
		for (i = 1; i <= 10; i++)
		{
			glVertex2f(i * 0.000001, -0.01);
			glVertex2f(i * 0.000001, 0.01);
		}
	glEnd();

	for (i = 1; i <= 10; i++)
	{
		sprintf(buf, "%f", lattice * i);
	    glRasterPos2f (lattice * i - 0.0000005, -0.05);
	    drawStringBig (buf);
	}

	glBegin(GL_LINES);
		for (i = 1; i <= 10; i++)
		{
			glVertex2f(-0.0000002, 0.1 * i);
			glVertex2f(0.0000002, 0.1 * i);
		}
	glEnd();

	for (i = 1; i <= 10; i++)
	{
		sprintf(buf, "%.1f", 0.1 * i);
	    glRasterPos2f (-0.0000004, 0.1 * i - 0.05);
	    drawStringBig (buf);
	}
}

void draw_differential_eq_plot(void)
{
	static char label[100];

	float x = 0, x2 = 0, y2, cx, cy;
	float tmp;
	int cache = 0;
	int i;

	glColor3f(0, 1, 1);
	glBegin(GL_LINES);
	//for(; ; t += step)
	for (i = 0; i < 100000000; i++)
	{
		x2 = plot_res[i].x;
		y2 = plot_res[i].y;

		if(cache)
		{
			glVertex2f(cx, cy);	// 이전값
			glVertex2f(x2, y2);	// 현재값
		}

		cache = 1;
		cx = x2;
		cy = y2;
		//printf("t = %f, y2 = %f\n", x2, y2);
	}
	glEnd();

	sprintf (label, "i(t)");
	// -0.000005
    glRasterPos2f (-0.0000045, 0.0000045);
    drawStringBig (label);
}

void draw_lattice2(void)
{
	int i;
	float lattice = 0.000001;
	char buf[64] = {0};

	glColor3f(1, 0, 0);
	glBegin(GL_LINES);
		for (i = 1; i <= 10; i++)
		{
			glVertex2f(i * 0.000001, -20000);
			glVertex2f(i * 0.000001, 20000);
		}
	glEnd();

	for (i = 1; i <= 10; i++)
	{
		sprintf(buf, "%f", lattice * i);
	    glRasterPos2f (lattice * i - 0.0000005, 50000);
	    drawStringBig (buf);
	}

	glBegin(GL_LINES);
		for (i = 1; i <= 10; i++)
		{
			glVertex2f(-0.0000002, -200000 * i);
			glVertex2f(0.0000002, -200000 * i);
		}
	glEnd();

	for (i = 1; i <= 10; i++)
	{
		sprintf(buf, "%d", -200000 * i);
	    glRasterPos2f (-0.0000015, -200000 * i);
	    drawStringBig (buf);
	}

}

void draw_volt_plot(void)
{
	static char label[100];

	float x = 0, x2 = 0, y2, cx, cy;
	float tmp;
	int cache = 0;
	int i;

	glColor3f(1, 1, 0);
	glBegin(GL_LINES);
	//for(; ; t += step)
	for (i = 0; i < 100000000; i++)
	{
		x2 = plot_res[i].x;
		y2 = plot_res[i].y * (-resistance / inductance);

		if(cache)
		{
			glVertex2f(cx, cy);	// 이전값
			glVertex2f(x2, y2);	// 현재값
		}

		cache = 1;
		cx = x2;
		cy = y2;
		//printf("t = %f, y2 = %f\n", x2, y2);
	}
	glEnd();

	sprintf (label, "V(t)");
	// -0.000005
    glRasterPos2f (-0.0000045, 0.000004);
    drawStringBig (label);
}

int main (int argc, char **argv)
{
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE);
	glutInitWindowSize(800, 800);
	glutInitWindowPosition(0, 0);
	glutCreateWindow("Digital Signal Processing");

	printf("저항 1000, 인덕터: 0.0003, 전압: 5V\n");

	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutMouseFunc(on_mouse);
	glutMotionFunc(drag_mouse);
	// 새롭게 추가된 요것은 키보드 입력을 처리하는 콜백 등록 함수
	glutKeyboardFunc(keyboard_handler);
	glutMainLoop();

    return 0;
}

/*************************** rhs ****************************/
/*
   x' = v                  ==>  dy[0]/dt = f[0] = y[1]
   v' = -x + \mu v (1-x^2) ==>  dy[1]/dt = f[1] = -y[0] + mu*y[1]*(1-y[0]*y[0])

   x''(t) + x'(t) + x(t) = 0
   x''(t) = -x'(t) - x(t)

   x' = v            ===> dy[0]/dt = f[0] = y[1]
   v' = -v - x       ===> dy[1]/dt = f[1] = -y[1] - y[0]
*/

// https://www.wolframalpha.com/input/?i=y%27%27+%2B+y%27+%2B+y+%3D+0%2C+y%280%29+%3D+1%2C+y%27%280%29+%3D+0, y(0) = 1
// https://www.wolframalpha.com/input/?i=%28sqrt%283%29+*+sin%28sqrt%283%29+%2F+2%29+%2B+3+*+cos%28sqrt%283%29+%2F+2%29%29+%2F+%283+*+sqrt%28e%29%29, y(1)
// (sqrt(3) * sin(sqrt(3) / 2) + 3 * cos(sqrt(3) / 2)) / (3 * sqrt(e)), y(1) = 0.6597001...
int
rhs (double t, const double y[], double f[], void *params_ptr)
{
#if 0
	f[0] = y[1];
	f[1] = -y[1] - 9 * y[0];
#endif

	f[0] =  -(resistance / inductance) * y[0] + voltage / inductance;

    return GSL_SUCCESS;		/* GSL_SUCCESS defined in gsl/errno.h as 0 */
}

#if 0

주석: 현재 아래 상황은 방전 상황입니다.
      또한 현재 코드는 충전 상황입니다.
	  두 가지 버전으로 작성이 필요함.

V = iR + L * di /dt
초기조건 V_L(t) | (t=0) ===> V_L(0) = 0V

V * (1 / L) = i(t) * (R / L) + i`(t)

P(x) = R / L

mu = C * exp^(int P(x)dx)

mu = C * exp^(R/L)t

d (mu * i(t)) / dt = (R/L) exp^(R/L)t * i(t) + exp^(R/L)t * i`(t)

int d(mu * i(t)) / dt = mu * i(t) = c

exp^(R/L)t * i(t) = C

i(t) = C * exp^(-R/L)t

i(t) | (t=0) = 0.5A ===> C = 0.5

i(t) = 0.5 * exp^(-R/L)t
V(t) = L di / dt
     = L * 0.5 * (-R/L) * exp^(-R/L)t
	-(resistance / inductance)
#endif


text:
/*
** get transformation matrix for co-registration
** Translation, rotation and one scaling parameter (7 degrees of freedom)
**
** G.Lohmann, 2015
*/
#include <viaio/VImage.h>
#include <viaio/Vlib.h>
#include <viaio/mu.h>
#include <viaio/os.h>

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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>

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

#define PI 3.14159265358979323846

double **MatrixAlloc(int n,int m)
{
  int i,j;
  double **A = (double **)VCalloc(n,sizeof(double));
  for(i=0; i<n; i++) {
    A[i]=(double *)VCalloc(m,sizeof(double));
  }
  for (i=0; i<n; i++) {
    for (j=0; j<m; j++) {
      A[i][j] = 0.0;
    }
  }
  return A;
}

void MatrixFree(double **A,int n,int m)
{
  int i;
  for(i=0; i<n; i++) VFree(A[i]);
  VFree(A);
}

/* matrix mult of 3x3 rotation matrices */
void MatrixMult3x3(double **A, double **B, double **C)
{
  int i,j,k,n=3;
  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
      double sum = 0;
      for (k=0; k<n; k++) {
        sum += A[i][k] * B[k][j];
      }
      C[i][j] = sum;
    }
  }
}

/* matrix ini */
void MatrixIni3x3(double **A)
{
  int i,j,n=3;
  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
      A[i][j] = 0.0;
    }
    for (i=0; i<n; i++) A[i][i] = 1.0;
  }
}

/* 3x3 matrix times vector */
void MatrixVector3x3(double **A, double *x,double *y)
{
  int i,j,n=3;
  for (i=0; i<n; i++) {
    double sum = 0;
    for (j=0; j<n; j++) {
      sum += A[i][j]*x[j];
    }
    y[i] = sum;
  }
}

void MatrixCopy3x3(double **A,double **B)
{
  int i,j,n=3;
  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
      B[i][j] = A[i][j];
    }
  }
}

void MatrixScale3x3(double **R,double sx,double sy,double sz)
{
  int i;
  for (i=0; i<3; i++) {
    R[i][0] *= sx;
    R[i][1] *= sy;
    R[i][2] *= sz;
  }
}

void printmat(VImage src)
{
  int i,j;
  fprintf(stderr,"printmat:\n");
  for (i=0; i<VImageNRows(src); i++) {
    for (j=0; j<VImageNColumns(src); j++) {
      float u = VGetPixel(src,0,i,j);
      fprintf(stderr," %8.3f",u);
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");
}

void xprintmat(double **A,int n,int m)
{
  int i,j;
  fprintf(stderr,"printmat:\n");
  for (i=0; i<n; i++) {
    for (j=0; j<m; j++) {
      fprintf(stderr," %8.3f",A[i][j]);
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");
}

void VCheckImage(VImage src)
{
  VAttrList out_list = VCreateAttrList();
  VAppendAttr(out_list,"image",NULL,VImageRepn,src);
  FILE *out_file = fopen("test.v","w");
  VWriteFile (out_file, out_list);
}

double clamp(double v)
{
  v = fmod(v, 2.0*PI);
  return v>=0 ? v : 2.0*PI+v;
}

void EulerAngles(double **R,double *ea)
{
  int i;

  double u = -R[0][2];
  if (u < -1.0) u = u + ceil(u);
  if (u >  1.0) u = u - floor(u);
  if (u < -1.0) u = -1.0;
  if (u >  1.0) u =  1.0;
  ea[1] = asin(u);

  double xcos = cos(ea[1]);
  if (ABS(xcos)<0.001*PI/180.0) {   /* Gimball lock? */
    ea[0] = atan2(R[1][0], R[1][1]);
    ea[2] = 0;
  }
  else {
    ea[0] = atan2(R[1][2], R[2][2]);
    ea[2] = atan2(R[0][1], R[0][0]);
  }
  for (i=0; i<3; i++) ea[i] = fmod(ea[i],2.0*PI);
  /*   for (i=0; i<3; i++) ea[i] = clamp(ea[i]); */
}

void Euler2RotationMatrix(double **R,double *ea)
{
  double c1 = cos(ea[0]);
  double c2 = cos(ea[1]);
  double c3 = cos(ea[2]);
  double s1 = sin(ea[0]);
  double s2 = sin(ea[1]);
  double s3 = sin(ea[2]);

  R[0][0] = c2*c1;
  R[1][0] = -c3*s1+s3*s2*c1;
  R[2][0] = s3*s1+c3*s2*c1;
  R[0][1] = c2*s1;
  R[1][1] = c3*c1+s3*s2*s1;
  R[2][1] = -s3*c1+c3*s2*s1;
  R[0][2] = -s2;
  R[1][2] = s3*c2;
  R[2][2] = c3*c2;

  double qfac = -1.0;
  R[0][2] *= qfac;
  R[1][2] *= qfac;
  R[2][2] *= qfac;
}

void Normalize(VImage src)
{
  int b,r,c;
  double tiny = 1.0e-6;
  if (VPixelRepn(src) == VUByteRepn || VPixelRepn(src) == VShortRepn
      || VPixelRepn(src) == VSByteRepn)tiny = 1.0;

  double sum1=0,sum2=0,nx=0;
  for (b=0; b<VImageNBands(src); b++) {
    for (r=0; r<VImageNRows(src); r++) {
      for (c=0; c<VImageNColumns(src); c++) {
	double u = VGetPixel(src,b,r,c);
	if (u >= tiny) {
	  sum1 += u;
	  sum2 += u*u;
	  nx++;
	}
      }
    }
  }
  double mean = sum1/nx;
  double sigma = sqrt((sum2 - nx * mean * mean) / (nx - 1.0));

  double xmax = VPixelMaxValue(src);
  double xmin = VPixelMinValue(src);
  double mean0 = (xmax + xmin)/2.0;
  double sig0 = xmax*0.2;
  fprintf(stderr," normalize, mean,sd= %f %f\n",mean0,sig0);

  for (b=0; b<VImageNBands(src); b++) {
    for (r=0; r<VImageNRows(src); r++) {
      for (c=0; c<VImageNColumns(src); c++) {
	double u = VGetPixel(src,b,r,c);
	if (u >= tiny) {
	  u = (u-mean)/sigma;
	  u = u*sig0 + mean0;
	  if (u > xmax) u = xmax;
	  if (u < xmin) u = xmin;
	  VSetPixel(src,b,r,c,u);
	}
      }
    }
  }
}

/* set up transformation info */
void ReadTransformation(VImage transimage,double **trans,double *shift,
			double *ref_pixdim,double *ref_dim,double *src_pixdim,double *src_dim)
{
  int i,j;
  for (i=0; i<3; i++) {
    for (j=0; j<3; j++) {
      trans[i][j] = VGetPixel(transimage,0,i,j);
    }
  }
  for (i=0; i<3; i++)
    shift[i] = VGetPixel(transimage,0,i,3);

  for (i=0; i<3; i++)
    ref_pixdim[i] = VGetPixel(transimage,0,i,4);

  for (i=0; i<3; i++)
    ref_dim[i] = VGetPixel(transimage,0,i,5);

  for (i=0; i<3; i++)
    src_pixdim[i] = VGetPixel(transimage,0,i,6);

  for (i=0; i<3; i++)
    src_dim[i] = VGetPixel(transimage,0,i,7);
}

void VSetResolution(VImage transimage,double **trans,double newres,double *shift,double *scale)
{
  /* resolution of reference image */
  double resx0 = VGetPixel(transimage,0,0,6);
  double resy0 = VGetPixel(transimage,0,1,6);
  double resz0 = VGetPixel(transimage,0,2,6);

  /* change resolution */
  scale[0] = resx0/newres;
  scale[1] = resy0/newres;
  scale[2] = resz0/newres;

  int i;
  MatrixScale3x3(trans,scale[0],scale[1],scale[2]);
  for (i=0; i<3; i++) shift[i] *= scale[i];
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/* specfunc/bessel_k.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 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 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"
#include "check.h"

#include "bessel.h"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* [Abramowitz+Stegun, 10.2.4 + 10.2.6]
 * with lmax=15, precision ~ 15D for x < 3
 *
 * assumes l >= 1
 */
static int bessel_kl_scaled_small_x(int l, const double x, gsl_sf_result * result)
{
  gsl_sf_result num_fact;
  double den  = gsl_sf_pow_int(x, l+1);
  int stat_df = gsl_sf_doublefact_e((unsigned int) (2*l-1), &num_fact);

  if(stat_df != GSL_SUCCESS || den == 0.0) {
    OVERFLOW_ERROR(result);
  }
  else {
    const int lmax = 50;
    gsl_sf_result ipos_term;
    double ineg_term;
    double sgn = (GSL_IS_ODD(l) ? -1.0 : 1.0);
    double ex  = exp(x);
    double t = 0.5*x*x;
    double sum = 1.0;
    double t_coeff = 1.0;
    double t_power = 1.0;
    double delta;
    int stat_il;
    int i;

    for(i=1; i<lmax; i++) {
      t_coeff /= i*(2*(i-l) - 1);
      t_power *= t;
      delta = t_power*t_coeff;
      sum += delta;
      if(fabs(delta/sum) < GSL_DBL_EPSILON) break;
    }

    stat_il = gsl_sf_bessel_il_scaled_e(l, x, &ipos_term);
    ineg_term =  sgn * num_fact.val/den * sum;
    result->val = -sgn * 0.5*M_PI * (ex*ipos_term.val - ineg_term);
    result->val *= ex;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_il;
  }
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_k0_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else {
    result->val = M_PI/(2.0*x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_k1_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < (M_SQRTPI+1.0)/(M_SQRT2*GSL_SQRT_DBL_MAX)) {
    OVERFLOW_ERROR(result);
  }
  else {
    result->val = M_PI/(2.0*x) * (1.0 + 1.0/x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_k2_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0/GSL_ROOT3_DBL_MAX) {
    OVERFLOW_ERROR(result);
  }
  else {
    result->val = M_PI/(2.0*x) * (1.0 + 3.0/x * (1.0 + 1.0/x));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_kl_scaled_e(int l, const double x, gsl_sf_result * result)
{
  if(l < 0 || x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(l == 0) {
    return gsl_sf_bessel_k0_scaled_e(x, result);
  }
  else if(l == 1) {
    return gsl_sf_bessel_k1_scaled_e(x, result);
  }
  else if(l == 2) {
    return gsl_sf_bessel_k2_scaled_e(x, result);
  }
  else if(x < 3.0) {
    return bessel_kl_scaled_small_x(l, x, result);
  }
  else if(GSL_ROOT3_DBL_EPSILON * x > (l*l + l + 1)) {
    int status = gsl_sf_bessel_Knu_scaled_asympx_e(l + 0.5, x, result);
    double pre = sqrt((0.5*M_PI)/x);
    result->val *= pre;
    result->err *= pre;
    return status;
  }
  else if(GSL_MIN(0.29/(l*l+1.0), 0.5/(l*l+1.0+x*x)) < GSL_ROOT3_DBL_EPSILON) {
    int status = gsl_sf_bessel_Knu_scaled_asymp_unif_e(l + 0.5, x, result);
    double pre = sqrt((0.5*M_PI)/x);
    result->val *= pre;
    result->err *= pre;
    return status;
  }
  else {
    /* recurse upward */
    gsl_sf_result r_bk;
    gsl_sf_result r_bkm;
    int stat_1 = gsl_sf_bessel_k1_scaled_e(x, &r_bk);
    int stat_0 = gsl_sf_bessel_k0_scaled_e(x, &r_bkm);
    double bkp;
    double bk  = r_bk.val;
    double bkm = r_bkm.val;
    int j;
    for(j=1; j<l; j++) {
      bkp = (2*j+1)/x*bk + bkm;
      bkm = bk;
      bk  = bkp;
    }
    result->val  = bk;
    result->err  = fabs(bk) * (fabs(r_bk.err/r_bk.val) + fabs(r_bkm.err/r_bkm.val));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_1, stat_0);
  }
}

int
gsl_sf_bessel_kl_scaled_array(const int lmax, const double x, double * result_array)
{
  if(lmax < 0 || x <= 0.0) {
    GSL_ERROR("domain error", GSL_EDOM);
  } else if (lmax == 0) {
    gsl_sf_result result;
    int stat = gsl_sf_bessel_k0_scaled_e(x, &result);
    result_array[0] = result.val;
    return stat;
  } else {
    int ell;
    double kellp1, kell, kellm1;
    gsl_sf_result r_kell;
    gsl_sf_result r_kellm1;
    gsl_sf_bessel_k1_scaled_e(x, &r_kell);
    gsl_sf_bessel_k0_scaled_e(x, &r_kellm1);
    kell   = r_kell.val;
    kellm1 = r_kellm1.val;
    result_array[0] = kellm1;
    result_array[1] = kell;
    for(ell = 1; ell < lmax; ell++) {
      kellp1 = (2*ell+1)/x * kell + kellm1;
      result_array[ell+1] = kellp1;
      kellm1 = kell;
      kell   = kellp1;
    }
    return GSL_SUCCESS;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_k0_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_k0_scaled_e(x, &result));
}

double gsl_sf_bessel_k1_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_k1_scaled_e(x, &result));
}

double gsl_sf_bessel_k2_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_k2_scaled_e(x, &result));
}

double gsl_sf_bessel_kl_scaled(const int l, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_kl_scaled_e(l, x, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_asum (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int N = 1;
   float X[] = { 0.239f };
   int incX = -1;
   float expected = 0.0f;
   float f;
   f = cblas_sasum(N, X, incX);
   gsl_test_rel(f, expected, flteps, "sasum(case 40)");
  };

  {
   int N = 1;
   double X[] = { -0.413 };
   int incX = -1;
   double expected = 0;
   double f;
   f = cblas_dasum(N, X, incX);
   gsl_test_rel(f, expected, dbleps, "dasum(case 41)");
  };

  {
   int N = 1;
   float X[] = { 0.1f, 0.017f };
   int incX = -1;
   float expected = 0.0f;
   float f;
   f = cblas_scasum(N, X, incX);
   gsl_test_rel(f, expected, flteps, "scasum(case 42)");
  };

  {
   int N = 1;
   double X[] = { -0.651, 0.079 };
   int incX = -1;
   double expected = 0;
   double f;
   f = cblas_dzasum(N, X, incX);
   gsl_test_rel(f, expected, dbleps, "dzasum(case 43)");
  };

  {
   int N = 2;
   float X[] = { 0.899f, -0.72f };
   int incX = 1;
   float expected = 1.619f;
   float f;
   f = cblas_sasum(N, X, incX);
   gsl_test_rel(f, expected, flteps, "sasum(case 44)");
  };

  {
   int N = 2;
   double X[] = { 0.271, -0.012 };
   int incX = 1;
   double expected = 0.283;
   double f;
   f = cblas_dasum(N, X, incX);
   gsl_test_rel(f, expected, dbleps, "dasum(case 45)");
  };

  {
   int N = 2;
   float X[] = { -0.567f, -0.645f, 0.098f, 0.256f };
   int incX = 1;
   float expected = 1.566f;
   float f;
   f = cblas_scasum(N, X, incX);
   gsl_test_rel(f, expected, flteps, "scasum(case 46)");
  };

  {
   int N = 2;
   double X[] = { -0.046, -0.671, -0.323, 0.785 };
   int incX = 1;
   double expected = 1.825;
   double f;
   f = cblas_dzasum(N, X, incX);
   gsl_test_rel(f, expected, dbleps, "dzasum(case 47)");
  };

  {
   int N = 2;
   float X[] = { 0.169f, 0.833f };
   int incX = -1;
   float expected = 0.0f;
   float f;
   f = cblas_sasum(N, X, incX);
   gsl_test_rel(f, expected, flteps, "sasum(case 48)");
  };

  {
   int N = 2;
   double X[] = { -0.586, -0.486 };
   int incX = -1;
   double expected = 0;
   double f;
   f = cblas_dasum(N, X, incX);
   gsl_test_rel(f, expected, dbleps, "dasum(case 49)");
  };

  {
   int N = 2;
   float X[] = { -0.314f, -0.318f, -0.835f, -0.807f };
   int incX = -1;
   float expected = 0.0f;
   float f;
   f = cblas_scasum(N, X, incX);
   gsl_test_rel(f, expected, flteps, "scasum(case 50)");
  };

  {
   int N = 2;
   double X[] = { -0.927, 0.152, -0.554, -0.844 };
   int incX = -1;
   double expected = 0;
   double f;
   f = cblas_dzasum(N, X, incX);
   gsl_test_rel(f, expected, dbleps, "dzasum(case 51)");
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


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

#define EXPECT_EQ(x, y) do {                                            \
        if (fabs(x - y) >= 1e-8)                                        \
            fprintf(stderr, "%d:FAILED: %f != %f\n",                    \
                    __LINE__, (double) x, (double) y);                  \
    } while (0)

int main() {
    double a[4] = { 1,  2,
                    3,  4};
    double b[4] = {-1, -2,
                   -3, -4};
    double c[4];

    /* test matrix multiplication */
    bls_dgemm('t', 't', 2, 2, 2, 1, a, 2, b, 2, 1, c, 2);
    EXPECT_EQ(c[0], -7);      EXPECT_EQ(c[2], -10);
    EXPECT_EQ(c[1], -15);     EXPECT_EQ(c[3], -22);

    /* test absolute sum */
    EXPECT_EQ(bls_dasum(4, a, 1), 10);

    return 0;
}

max_stars_repo_name: Rufflewind/blas-shim

max_issues_repo_name: Rufflewind/blas-shim

max_forks_repo_name: Rufflewind/blas-shim


text:
/**
 *
 * @file test_dotu.c
 *
 * @author Samuel  D. Relton
 * @author Pedro   V. Lara
 * @author Mawussi Zounon
 * @date 2016-06-15
 * @precisions normal z -> c
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include<time.h>
#ifdef BBLAS_WITH_MKL
    #include <mkl_cblas.h>
    #include <mkl_lapacke.h>
#else
    #include <cblas.h>
    #include <lapacke.h>
#endif
#include "bblas_z.h"

#undef REAL
#define COMPLEX

int main(int argc, char *argv[])
{
  //Set parameters to initilize random value generation
  __time_t t;
  srand((unsigned) time(&t));
  int IONE     = 1;
  int ISEED[4] ={0,0,0,1};

  // Generation of batch_count value between 100 and 1000
  int batch_min = 100;
  int batch_max = 10000;
  int batch_count = rand() % (batch_max - batch_min + 1) + batch_min;

  // Generation of x and y vector sizes
  int * n = (int *) malloc(batch_count*sizeof(int));
  int n_max = 1000;
  int n_min = 10;
  for(int iter=0; iter < batch_count; iter++)
    n[iter] = rand() % (n_max - n_min + 1) + n_min;

  // common calling parameters
  int * incx = (int*) malloc(batch_count*sizeof(int));
  int * incy = (int*) malloc(batch_count*sizeof(int));
  int * info = (int *) malloc(batch_count*sizeof(int));
  enum BBLAS_OPTS batch_opts = BBLAS_VARIABLE;

  //Generation of x and y vectors  and the result vector
  // for using variable api
  BBLAS_Complex64_t ** x, ** y, *result;
  x = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
  y = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
  result = (BBLAS_Complex64_t *) malloc(batch_count*sizeof(BBLAS_Complex64_t));

  for(int iter=0; iter < batch_count; iter++){
    x[iter] = (BBLAS_Complex64_t*) malloc(n[iter]*sizeof(BBLAS_Complex64_t));
    y[iter] = (BBLAS_Complex64_t*) malloc(n[iter]*sizeof(BBLAS_Complex64_t));

#if defined(BBLAS_WITH_MKL)
    LAPACKE_zlarnv_work(IONE, ISEED, n[iter], (MKL_Complex16*) x[iter]);
    LAPACKE_zlarnv_work(IONE, ISEED, n[iter], (MKL_Complex16*) y[iter]);
#else
    LAPACKE_zlarnv_work(IONE, ISEED, n[iter], x[iter]);
    LAPACKE_zlarnv_work(IONE, ISEED, n[iter], y[iter]);
#endif
    // set values of incx and incy
    incx[iter] = 1;
    incy[iter] = 1;
  }

  //Calling variable api
  printf("Calling variable api: batchv_zdotu_sub \n");
  batchv_zdotu_sub(
		  (const int *) n,
		  (BBLAS_Complex64_t const * const *) x,
		  (const int *) incx,
		  (BBLAS_Complex64_t const * const *) y,
		  (const int *) incy,
		  result,
		  (const int) batch_count, info);

  //Calling tag api to solve variable problem
  printf("Calling tag api to solve variable problem: batch_zdotu_sub\n");

    batch_zdotu_sub(
		  (const int *) n,
		  (BBLAS_Complex64_t const * const *) x,
		  (const int *) incx,
		  (BBLAS_Complex64_t const * const *) y,
		  (const int *) incy,
		  result,
		  (const int) batch_count,
		  (const enum BBLAS_OPTS) batch_opts,
		  info);

    //Free x and y vectors
    printf("Free x and y vectors after variable case\n");

    for(int iter=0; iter < batch_count; iter++){
      free(x[iter]);
      free(y[iter]);
    }
    // free n incx and incy and info
    free(n);
    free(incx);
    free(incy);
    free(info);

    // Calling fixe api
    //Generation of x and y vectors, result will be overwritten
    int n_fixe = rand() % (n_max - n_min + 1) + n_min;
    int incx_fixe = 1;
    int incy_fixe = 1;
    int * info_fixe = (int*)malloc(sizeof(int));

    for(int iter=0; iter < batch_count; iter++){
      x[iter] = (BBLAS_Complex64_t*) malloc(n_fixe*sizeof(BBLAS_Complex64_t));
      y[iter] = (BBLAS_Complex64_t*) malloc(n_fixe*sizeof(BBLAS_Complex64_t));

#if defined(BBLAS_WITH_MKL)
      LAPACKE_zlarnv_work(IONE, ISEED, n_fixe, (MKL_Complex16*) x[iter]);
      LAPACKE_zlarnv_work(IONE, ISEED, n_fixe, (MKL_Complex16*) y[iter]);
#else
      LAPACKE_zlarnv_work(IONE, ISEED, n_fixe, x[iter]);
      LAPACKE_zlarnv_work(IONE, ISEED, n_fixe, y[iter]);
#endif
    }
  //Calling fixe api
    printf("Calling fixe api: batchf_zdotu_sub \n");
  batchf_zdotu_sub(
		  (const int ) n_fixe,
		  (BBLAS_Complex64_t const * const *) x,
		  (const int) incx_fixe,
		  (BBLAS_Complex64_t const * const *) y,
		  (const int) incy_fixe,
		  result,
		  (const int) batch_count, info_fixe);

  //Calling tag api to solve fixed batch  problem
  printf("Calling tag api to solve fixed batch  problem: batch_zdotu_sub \n");
  batch_opts = BBLAS_FIXED;

  batch_zdotu_sub(
		  (const int *) &n_fixe,
		  (BBLAS_Complex64_t const * const *) x,
		  (const int *) &incx_fixe,
		  (BBLAS_Complex64_t const * const *) y,
		  (const int *) &incy_fixe,
		  result,
		  (const int) batch_count,
		  (const enum BBLAS_OPTS) batch_opts,
		  info_fixe);
  //Free x and y vectors
  printf("Free memory after the call to fixed batch case\n");

  for(int iter=0; iter < batch_count; iter++){
    free(x[iter]);
    free(y[iter]);
  }

  //Creation of variable specific to group api
  int group_min = 1;
  int group_max = 10;
  int group_count = rand() % (group_max - group_min + 1) + group_min;
  int * group_size = (int *) malloc(group_count*sizeof(int));
  int tmp_size = batch_count/group_count;

  //Set value for group_size
  for(int iter=0; iter < group_count-1; iter++)
    group_size[iter] = tmp_size;

  //Set last value of group_size
  group_size[group_count-1] = batch_count - (tmp_size*(group_count-1));

  //Allocate memory for vectors

  int *n_group    = (int *) malloc(group_count*sizeof(int));
  int *incx_group = (int *) malloc(group_count*sizeof(int));
  int *incy_group = (int *) malloc(group_count*sizeof(int));
  int *info_group = (int *) malloc(group_count*sizeof(int));
  int offset = 0;
  //Set vector entries

  for(int group_iter = 0; group_iter < group_count ; group_iter++){
    n_group[group_iter] = rand() % (n_max - n_min + 1) + n_min;
    incx_group[group_iter] = 1;
    incy_group[group_iter] = 1;

    for (int iter = 0; iter < group_size[group_iter]; iter++){
      x[offset+iter] = (BBLAS_Complex64_t*) malloc(n_group[group_iter]*sizeof(BBLAS_Complex64_t));
      y[offset+iter] = (BBLAS_Complex64_t*) malloc(n_group[group_iter]*sizeof(BBLAS_Complex64_t));

#if defined(BBLAS_WITH_MKL)
      LAPACKE_zlarnv_work(IONE, ISEED, n_group[group_iter], (MKL_Complex16*) x[offset+iter]);
      LAPACKE_zlarnv_work(IONE, ISEED, n_group[group_iter], (MKL_Complex16*) y[offset+iter]);
#else
      LAPACKE_zlarnv_work(IONE, ISEED, n_group[group_iter], x[offset+iter]);
      LAPACKE_zlarnv_work(IONE, ISEED, n_group[group_iter], y[offset+iter]);
#endif
    }
    offset += group_size[group_iter];
  }

  //Calling group api
  printf("Calling group api:batchg_zdotu_sub\n");
   batchg_zdotu_sub(
		    (const int *)n_group,
		    (BBLAS_Complex64_t const * const *)x,
		    (const int *)incx_group,
		    (BBLAS_Complex64_t const * const *)y,
		    (const int *)incy_group,
		    result,
		    (const int)group_count,
		    (const int *) group_size,
		    info_group);

   //Free allocated memories
  for(int iter=0; iter < batch_count; iter++){
    free(x[iter]);
    free(y[iter]);
  }
  free(result);
  free(x);
  free(y);
  free(n_group);
  free(incx_group);
  free(incy_group);
  free(group_size);
  free(info_group);
  return 0;
}

max_stars_repo_name: mawussi/BBLAS-group


text:
/* specfunc/exp.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 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 */

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

#include "error.h"

/* Evaluate the continued fraction for exprel.
 * [Abramowitz+Stegun, 4.2.41]
 */
static
int
exprel_n_CF(const double N, const double x, gsl_sf_result * result)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 5000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = 1.0;
  double b1 = 1.0;
  double a2 = -x;
  double b2 = N+1;
  double an, bn;

  double fn;

  double An = b1*Anm1 + a1*Anm2;   /* A1 */
  double Bn = b1*Bnm1 + a1*Bnm2;   /* B1 */

  /* One explicit step, before we get to the main pattern. */
  n++;
  Anm2 = Anm1;
  Bnm2 = Bnm1;
  Anm1 = An;
  Bnm1 = Bn;
  An = b2*Anm1 + a2*Anm2;   /* A2 */
  Bn = b2*Bnm1 + a2*Bnm2;   /* B2 */

  fn = An/Bn;

  while(n < maxiter) {
    double old_fn;
    double del;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    an = ( GSL_IS_ODD(n) ? ((n-1)/2)*x : -(N+(n/2)-1)*x );
    bn = N + n - 1;
    An = bn*Anm1 + an*Anm2;
    Bn = bn*Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An /= RECUR_BIG;
      Bn /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;

    if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
  }

  result->val = fn;
  result->err = 4.0*(n+1.0)*GSL_DBL_EPSILON*fabs(fn);

  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_exp_e(const double x, gsl_sf_result * result)
{
  if(x > GSL_LOG_DBL_MAX) {
    OVERFLOW_ERROR(result);
  }
  else if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else {
    result->val = exp(x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result)
{
  if(x > INT_MAX-1) {
    OVERFLOW_ERROR_E10(result);
  }
  else if(x < INT_MIN+1) {
    UNDERFLOW_ERROR_E10(result);
  }
  else {
    const int N = (x > GSL_LOG_DBL_MAX || x < GSL_LOG_DBL_MIN) ? (int) floor(x/M_LN10) : 0;
    result->val = exp(x-N*M_LN10);
    result->err = 2.0 * (fabs(x)+1.0) * GSL_DBL_EPSILON * fabs(result->val);
    result->e10 = N;
    return GSL_SUCCESS;
  }
}

int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result)
{
  const double ay  = fabs(y);

  if(y == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
    ) {
    const double ex = exp(x);
    result->val = y * ex;
    result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double ly  = log(ay);
    const double lnr = x + ly;

    if(lnr > GSL_LOG_DBL_MAX - 0.01) {
      OVERFLOW_ERROR(result);
    }
    else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
      UNDERFLOW_ERROR(result);
    }
    else {
      const double sy   = GSL_SIGN(y);
      const double M    = floor(x);
      const double N    = floor(ly);
      const double a    = x  - M;
      const double b    = ly - N;
      const double berr = 2.0 * GSL_DBL_EPSILON * (fabs(ly) + fabs(N));
      result->val  = sy * exp(M+N) * exp(a+b);
      result->err  = berr * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * (M + N + 1.0) * fabs(result->val);
      return GSL_SUCCESS;
    }
  }
}

int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result)
{
  const double ay  = fabs(y);

  if(y == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
    ) {
    const double ex = exp(x);
    result->val = y * ex;
    result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else {
    const double ly  = log(ay);
    const double l10_val = (x + ly)/M_LN10;

    if(l10_val > INT_MAX-1) {
      OVERFLOW_ERROR_E10(result);
    }
    else if(l10_val < INT_MIN+1) {
      UNDERFLOW_ERROR_E10(result);
    }
    else {
      const double sy  = GSL_SIGN(y);
      const int    N   = (int) floor(l10_val);
      const double arg_val = (l10_val - N) * M_LN10;
      const double arg_err = 2.0 * GSL_DBL_EPSILON * (fabs(x) + fabs(ly) + M_LN10*fabs(N));

      result->val  = sy * exp(arg_val);
      result->err  = arg_err * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      result->e10 = N;

      return GSL_SUCCESS;
    }
  }
}

int gsl_sf_exp_mult_err_e(const double x, const double dx,
                             const double y, const double dy,
                             gsl_sf_result * result)
{
  const double ay  = fabs(y);

  if(y == 0.0) {
    result->val = 0.0;
    result->err = fabs(dy * exp(x));
    return GSL_SUCCESS;
  }
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
    ) {
    double ex = exp(x);
    result->val  = y * ex;
    result->err  = ex * (fabs(dy) + fabs(y*dx));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double ly  = log(ay);
    const double lnr = x + ly;

    if(lnr > GSL_LOG_DBL_MAX - 0.01) {
      OVERFLOW_ERROR(result);
    }
    else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
      UNDERFLOW_ERROR(result);
    }
    else {
      const double sy  = GSL_SIGN(y);
      const double M   = floor(x);
      const double N   = floor(ly);
      const double a   = x  - M;
      const double b   = ly - N;
      const double eMN = exp(M+N);
      const double eab = exp(a+b);
      result->val  = sy * eMN * eab;
      result->err  = eMN * eab * 2.0*GSL_DBL_EPSILON;
      result->err += eMN * eab * fabs(dy/y);
      result->err += eMN * eab * fabs(dx);
      return GSL_SUCCESS;
    }
  }
}

int gsl_sf_exp_mult_err_e10_e(const double x, const double dx,
                             const double y, const double dy,
                             gsl_sf_result_e10 * result)
{
  const double ay  = fabs(y);

  if(y == 0.0) {
    result->val = 0.0;
    result->err = fabs(dy * exp(x));
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
    ) {
    const double ex = exp(x);
    result->val  = y * ex;
    result->err  = ex * (fabs(dy) + fabs(y*dx));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else {
    const double ly  = log(ay);
    const double l10_val = (x + ly)/M_LN10;

    if(l10_val > INT_MAX-1) {
      OVERFLOW_ERROR_E10(result);
    }
    else if(l10_val < INT_MIN+1) {
      UNDERFLOW_ERROR_E10(result);
    }
    else {
      const double sy  = GSL_SIGN(y);
      const int    N   = (int) floor(l10_val);
      const double arg_val = (l10_val - N) * M_LN10;
      const double arg_err = dy/fabs(y) + dx + 2.0*GSL_DBL_EPSILON*fabs(arg_val);

      result->val  = sy * exp(arg_val);
      result->err  = arg_err * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      result->e10 = N;

      return GSL_SUCCESS;
    }
  }
}

int gsl_sf_expm1_e(const double x, gsl_sf_result * result)
{
  const double cut = 0.002;

  if(x < GSL_LOG_DBL_MIN) {
    result->val = -1.0;
    result->err = GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < -cut) {
    result->val = exp(x) - 1.0;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < cut) {
    result->val = x * (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_LOG_DBL_MAX) {
    result->val = exp(x) - 1.0;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

int gsl_sf_exprel_e(const double x, gsl_sf_result * result)
{
  const double cut = 0.002;

  if(x < GSL_LOG_DBL_MIN) {
    result->val = -1.0/x;
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < -cut) {
    result->val = (exp(x) - 1.0)/x;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < cut) {
    result->val = (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_LOG_DBL_MAX) {
    result->val = (exp(x) - 1.0)/x;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

int gsl_sf_exprel_2_e(double x, gsl_sf_result * result)
{
  const double cut = 0.002;

  if(x < GSL_LOG_DBL_MIN) {
    result->val = -2.0/x*(1.0 + 1.0/x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < -cut) {
    result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < cut) {
    result->val = (1.0 + 1.0/3.0*x*(1.0 + 0.25*x*(1.0 + 0.2*x*(1.0 + 1.0/6.0*x))));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_LOG_DBL_MAX) {
    result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

int
gsl_sf_exprel_n_CF_e(const double N, const double x, gsl_sf_result * result)
{
  return exprel_n_CF(N, x, result);
}

int
gsl_sf_exprel_n_e(const int N, const double x, gsl_sf_result * result)
{
  if(N < 0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(fabs(x) < GSL_ROOT3_DBL_EPSILON * N) {
    result->val = 1.0 + x/(N+1) * (1.0 + x/(N+2));
    result->err = 2.0 * GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(N == 0) {
    return gsl_sf_exp_e(x, result);
  }
  else if(N == 1) {
    return gsl_sf_exprel_e(x, result);
  }
  else if(N == 2) {
    return gsl_sf_exprel_2_e(x, result);
  }
  else {
    if(x > N && (-x + N*(1.0 + log(x/N)) < GSL_LOG_DBL_EPSILON)) {
      /* x is much larger than n.
       * Ignore polynomial part, so
       * exprel_N(x) ~= e^x N!/x^N
       */
      gsl_sf_result lnf_N;
      double lnr_val;
      double lnr_err;
      double lnterm;
      gsl_sf_lnfact_e(N, &lnf_N);
      lnterm = N*log(x);
      lnr_val  = x + lnf_N.val - lnterm;
      lnr_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(lnterm));
      lnr_err += lnf_N.err;
      return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
    }
    else if(x > N) {
      /* Write the identity
       *   exprel_n(x) = e^x n! / x^n (1 - Gamma[n,x]/Gamma[n])
       * then use the asymptotic expansion
       * Gamma[n,x] ~ x^(n-1) e^(-x) (1 + (n-1)/x + (n-1)(n-2)/x^2 + ...)
       */
      double ln_x = log(x);
      gsl_sf_result lnf_N;
      double lg_N;
      double lnpre_val;
      double lnpre_err;
      gsl_sf_lnfact_e(N, &lnf_N);    /* log(N!)       */
      lg_N  = lnf_N.val - log(N);       /* log(Gamma(N)) */
      lnpre_val  = x + lnf_N.val - N*ln_x;
      lnpre_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(N*ln_x));
      lnpre_err += lnf_N.err;
      if(lnpre_val < GSL_LOG_DBL_MAX - 5.0) {
        int stat_eG;
        gsl_sf_result bigG_ratio;
        gsl_sf_result pre;
        int stat_ex = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &pre);
        double ln_bigG_ratio_pre = -x + (N-1)*ln_x - lg_N;
        double bigGsum = 1.0;
        double term = 1.0;
        int k;
        for(k=1; k<N; k++) {
          term *= (N-k)/x;
          bigGsum += term;
        }
        stat_eG = gsl_sf_exp_mult_e(ln_bigG_ratio_pre, bigGsum, &bigG_ratio);
        if(stat_eG == GSL_SUCCESS) {
          result->val  = pre.val * (1.0 - bigG_ratio.val);
          result->err  = pre.val * (2.0*GSL_DBL_EPSILON + bigG_ratio.err);
          result->err += pre.err * fabs(1.0 - bigG_ratio.val);
          result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
          return stat_ex;
        }
        else {
          result->val = 0.0;
          result->err = 0.0;
          return stat_eG;
        }
      }
      else {
        OVERFLOW_ERROR(result);
      }
    }
    else if(x > -10.0*N) {
      return exprel_n_CF(N, x, result);
    }
    else {
      /* x -> -Inf asymptotic:
       * exprel_n(x) ~ e^x n!/x^n - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
       *             ~ - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
       */
      double sum  = 1.0;
      double term = 1.0;
      int k;
      for(k=1; k<N; k++) {
        term *= (N-k)/x;
        sum  += term;
      }
      result->val = -N/x * sum;
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
  }
}

int
gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result)
{
  const double adx = fabs(dx);

  /* CHECK_POINTER(result) */

  if(x + adx > GSL_LOG_DBL_MAX) {
    OVERFLOW_ERROR(result);
  }
  else if(x - adx < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else {
    const double ex  = exp(x);
    const double edx = exp(adx);
    result->val  = ex;
    result->err  = ex * GSL_MAX_DBL(GSL_DBL_EPSILON, edx - 1.0/edx);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result)
{
  const double adx = fabs(dx);

  /* CHECK_POINTER(result) */

  if(x + adx > INT_MAX - 1) {
    OVERFLOW_ERROR_E10(result);
  }
  else if(x - adx < INT_MIN + 1) {
    UNDERFLOW_ERROR_E10(result);
  }
  else {
    const int    N  = (int)floor(x/M_LN10);
    const double ex = exp(x-N*M_LN10);
    result->val = ex;
    result->err = ex * (2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) + adx);
    result->e10 = N;
    return GSL_SUCCESS;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_exp(const double x)
{
  EVAL_RESULT(gsl_sf_exp_e(x, &result));
}

double gsl_sf_exp_mult(const double x, const double y)
{
  EVAL_RESULT(gsl_sf_exp_mult_e(x, y, &result));
}

double gsl_sf_expm1(const double x)
{
  EVAL_RESULT(gsl_sf_expm1_e(x, &result));
}

double gsl_sf_exprel(const double x)
{
  EVAL_RESULT(gsl_sf_exprel_e(x, &result));
}

double gsl_sf_exprel_2(const double x)
{
  EVAL_RESULT(gsl_sf_exprel_2_e(x, &result));
}

double gsl_sf_exprel_n(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_exprel_n_e(n, x, &result));
}


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * 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.
 */

#ifndef SIM_PHYSICS_AERO_H_
#define SIM_PHYSICS_AERO_H_

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <stdint.h>

#include <memory>
#include <string>
#include <vector>

#include "common/c_math/force_moment.h"
#include "common/c_math/linalg.h"
#include "common/c_math/vec3.h"
#include "common/macros.h"
#include "control/system_types.h"
#include "sim/physics/aero_types.h"
#include "sim/physics/avl_aero_database.h"
#include "sim/physics/delta_coeff_aero_database.h"
#include "sim/physics/dvl_aero_database.h"
#include "sim/sim_types.h"

class Aero {
  friend class AeroTest;

 public:
  explicit Aero(const AeroSimParams &aero_sim_params);
  ~Aero() {}

  // Calculates the force-moment coefficient for a given (alpha, beta,
  // omega_hat, Reynolds number, flap deflection) operating point,
  // taking into account information from multiple databases.
  void CalcForceMomentCoeff(double alpha, double beta, const Vec3 &omega_hat,
                            const Vec &flaps, double reynolds_number,
                            ForceMoment *force_moment_coeff,
                            double thrust_coeff) const {
    CalcForceMomentCoeff(alpha, beta, omega_hat, flaps, reynolds_number,
                         force_moment_coeff, thrust_coeff, nullptr);
  }

  // See function above.  Optionally returns the breakdown of the
  // force-moment coefficient into stability and control derivatives
  // if raw_aero_coeffs is not null.
  void CalcForceMomentCoeff(double alpha, double beta, const Vec3 &omega_hat,
                            const Vec &flaps, double reynolds_number,
                            ForceMoment *force_moment_coeff,
                            double thrust_coeff,
                            RawAeroCoeffs *raw_aero_coeffs) const;

  void CalcLowIncidenceCoeffs(double alpha, double beta, const Vec3 &omega_hat,
                              const Vec &flaps, double reynolds_number,
                              ForceMoment *force_moment_coeff,
                              double thrust_coeff) const;

  void CalcHighIncidenceCoeffs(double alpha, double beta, const Vec3 &omega_hat,
                               const Vec &flaps, double reynolds_number,
                               ForceMoment *force_moment_coeff,
                               double thrust_coeff) const;

  double CalcHingeCM(FlapLabel label, double alpha, double delta) const;

  double CalcEmpiricalPitchingMoment(double alpha) const;

 private:
  // Returns true if the parameters to the aerodynamic model pass
  // basic consistency checks.
  bool ValidateParams(const AeroSimParams &params) const;

  // Linearly interpolates the force-moment coefficients between two
  // aerodynamics databases based on Reynolds number.  Returns the
  // weight used to mix the databases of the lower Reynolds database.
  template <typename Database, typename DatabaseCoeffs>
  double InterpolateAeroDatabases(
      double alpha, double beta, const Vec3 &omega_hat, const Vec &flaps,
      double reynolds_number,
      const std::vector<std::unique_ptr<Database>> &databases,
      DatabaseCoeffs *lower_coeffs, DatabaseCoeffs *upper_coeffs,
      ForceMoment *cfm, double thrust_coeff) const;
  const AeroSimParams &aero_sim_params_;

  // Set of databases to use when the flaps are deflected in the
  // linear range.  These databases are typically used during
  // crosswind flight.
  std::vector<std::unique_ptr<AvlAeroDatabase>> small_deflection_databases_;

  // Set of databases to use when the flaps may be deflected a large
  // amount, outside the linear range.  These databases are typically
  // used during hover and transition when the elevator may be
  // feathered.
  std::vector<std::unique_ptr<DvlAeroDatabase>> large_deflection_databases_;

  // Set of databases to use when the spoilers are deflected to
  // max range.  These databases are typically used during
  // crosswind flight, during high speed events.
  std::unique_ptr<DeltaCoeffAeroDatabase> spoiler_offset_database_;

  DISALLOW_COPY_AND_ASSIGN(Aero);
};

#endif  // SIM_PHYSICS_AERO_H_

ext: h


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

#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics.h>

#define N (4096*4*16)

int main (void)
{
    long countsq;
    double x, r2;
    //6 dimensional space
    int dim = 6;

    gsl_rng *r;

    r = gsl_rng_alloc (gsl_rng_taus2);
    gsl_rng_set (r, 1UL);

    long i, j, nn;

    printf ("# Iter        Npoints   Volume     AbsErr\n");
    nn = N;
    for (j = 0; j < 16; j++)
    {
        countsq = 0L;
        for (i = 0; i < nn; i++)
        {
            // 4-dimensional distance from the origin, squared
            r2 = 0.;
            for (int k = 0; k < dim; k++)
            {
                x = 2 * gsl_rng_uniform (r) - 1.;
                r2 += x * x;
            }

            // are we inside the unit ball
            if (r2 <= 1.)
            {
                countsq++;
            }
        }
        // ratio of the volume of the 4-dimensional unit ball and
        // and 4-dimensional cube of side length 2.
        double vol = ((double) countsq) * pow (2, (double) dim) / nn;
        double exact = M_PI * M_PI / 2;

        printf ("%3ld  %16ld   %f   %f\n", j, nn, vol, fabs (exact - vol));
        nn *= 2;
    }

    gsl_rng_free (r);

    return 0;
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <dc1394/dc1394.h>
#include <math.h>
#include <fitsio.h>
#include <gsl/gsl_statistics.h>
#include <time.h>
#include <sys/time.h>
#include <string.h>
#include <inttypes.h>

#define REG_CAMERA_ABS_MIN                  0x000U
#define REG_CAMERA_ABS_MAX                  0x004U
#define REG_CAMERA_ABS_VALUE                0x008U

int grab_frame(dc1394camera_t *c, unsigned char *buf, int nbytes) {
  dc1394video_frame_t *frame=NULL;
  dc1394error_t err;

  err = dc1394_capture_dequeue(c, DC1394_CAPTURE_POLICY_WAIT, &frame);
  if (err != DC1394_SUCCESS) {
    dc1394_log_error("Unable to capture.");
    dc1394_capture_stop(c);
    dc1394_camera_free(c);
    exit(1);
  }

  memcpy(buf, frame->image, nbytes);
  dc1394_capture_enqueue(c, frame);
  return 1;
}

int main(int argc, char *argv[]) {

  long nelements, naxes[2];
  dc1394_t * dc;
  dc1394camera_t * camera;
  dc1394camera_list_t * list;
  dc1394error_t err;
  dc1394video_mode_t mode;
  unsigned int max_height, max_width, winleft, wintop;

  unsigned char *buffer, *buffer2;
  double *average, *diff;
  int j, f, status, nimages, xsize, ysize;
  double mean, var, exp;
  struct timeval start_time, end_time;
  time_t start_sec, end_sec;
  suseconds_t start_usec, end_usec;
  float elapsed_time, fps;

  stderr = freopen("noise_test.log", "w", stderr);

  srand48((unsigned)time(NULL));

  nimages = 0;
  status = 0;
  xsize = 320;
  ysize = 240;
  naxes[0] = xsize;
  naxes[1] = ysize;

  nelements = naxes[0]*naxes[1];

  mode = DC1394_VIDEO_MODE_FORMAT7_1;
  dc = dc1394_new();
  if (!dc)
    return -1;
  err = dc1394_camera_enumerate(dc, &list);
  DC1394_ERR_RTN(err, "Failed to enumerate cameras.");

  if (list->num == 0) {
    dc1394_log_error("No cameras found.");
    return -1;
  }

  camera = dc1394_camera_new(dc, list->ids[0].guid);
  if (!camera) {
    dc1394_log_error("Failed to initialize camera with guid %"PRIx64".",
                     list->ids[0].guid);
    return -1;
  }
  dc1394_camera_free_list(list);

  printf("Using camera with GUID %"PRIx64"\n", camera->guid);

  // need to use legacy firewire400 mode for now.  800 not quite reliable.
  dc1394_video_set_iso_speed(camera, DC1394_ISO_SPEED_400);

  // configure camera for format7
  err = dc1394_video_set_mode(camera, mode);
  DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't choose video mode.");

  err = dc1394_format7_get_max_image_size(camera, mode, &max_width, &max_height);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot get max image size.");

  winleft = 0;
  wintop = 0;

  err = dc1394_format7_set_roi(camera,
                               mode,
                               DC1394_COLOR_CODING_MONO8,
                               DC1394_USE_MAX_AVAIL,
                               //				 packet,
                               winleft, wintop, // left, top
                               xsize, ysize);
  DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't set ROI.");

  err=dc1394_video_set_framerate(camera, DC1394_FRAMERATE_MAX);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set framerate");

  err = dc1394_feature_set_mode (camera, DC1394_FEATURE_SHUTTER, DC1394_FEATURE_MODE_MANUAL);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set shutter to manual");

  err = dc1394_feature_set_value (camera, DC1394_FEATURE_SHUTTER, 61);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set shutter");

  err = dc1394_feature_set_mode (camera, DC1394_FEATURE_GAIN, DC1394_FEATURE_MODE_MANUAL);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set gain to manual");

  err = dc1394_feature_set_value (camera, DC1394_FEATURE_GAIN, 48);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set gain");

  // set brightness manually.  use relative value in range 0 to 1023.
  err = dc1394_feature_set_mode(camera, DC1394_FEATURE_BRIGHTNESS, DC1394_FEATURE_MODE_MANUAL);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set brightness to manual");
  err = dc1394_feature_set_value(camera, DC1394_FEATURE_BRIGHTNESS, 500);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set brightness");
  printf ("I: brightness is %d\n", 500);

  err = dc1394_capture_setup(camera, 16, DC1394_CAPTURE_FLAGS_DEFAULT);
  DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Error capturing.");

  // start the camera up
  err = dc1394_video_set_transmission(camera, DC1394_ON);
  if (err != DC1394_SUCCESS) {
    dc1394_log_error("Unable to start camera iso transmission.");
    dc1394_capture_stop(camera);
    dc1394_camera_free(camera);
    return -1;
  }

  printf("Camera successfully initialized.\n");

  if (!(buffer = malloc(nelements*sizeof(char)))) {
    printf("Couldn't Allocate Image Buffer\n");
    exit(-1);
  }

  if (!(buffer2 = malloc(nelements*sizeof(char)))) {
    printf("Couldn't Allocate Image2 Buffer\n");
    exit(-1);
  }

  if (!(average = calloc(nelements, sizeof(double)))) {
    printf("Couldn't Allocate Average Image Buffer\n");
    exit(-1);
  }

  if (!(diff = calloc(nelements, sizeof(double)))) {
    printf("Couldn't Allocate Diff Image Buffer\n");
    exit(-1);
  }

  gettimeofday(&start_time, NULL);

  f = 0;

  err = dc1394_feature_set_absolute_control(camera, DC1394_FEATURE_SHUTTER, DC1394_TRUE);
  DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set shutter to absolute mode");
  exp = 8.0*1.0e-5;
  err = dc1394_feature_set_absolute_value(camera, DC1394_FEATURE_SHUTTER, exp);
  grab_frame(camera, buffer, nelements*sizeof(char));
  grab_frame(camera, buffer, nelements*sizeof(char));

  for (f=8; f<320; f++) {
    exp = f*1.0e-5;
    err = dc1394_feature_set_absolute_value(camera, DC1394_FEATURE_SHUTTER, exp);
    DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set shutter");

    grab_frame(camera, buffer, nelements*sizeof(char));
    grab_frame(camera, buffer2, nelements*sizeof(char));
    for (j=0; j<nelements; j++) {
      average[j] = (buffer[j]+buffer2[j])/2.0;
    }
    for (j=0; j<nelements; j++) {
      diff[j] = buffer[j] - buffer2[j];
    }

    var = gsl_stats_variance(diff, 1, nelements);
    mean = gsl_stats_mean(average, 1, nelements);

    var = var/2.0;

    printf("%e %10.2f %10.3f\n", exp, mean, var);

  }

  gettimeofday(&end_time, NULL);
  printf("End capture.\n");

  start_sec = start_time.tv_sec;
  start_usec = start_time.tv_usec;
  end_sec = end_time.tv_sec;
  end_usec = end_time.tv_usec;

  elapsed_time = (float)((end_sec + 1.0e-6*end_usec) - (start_sec + 1.0e-6*start_usec));
  fps = nimages/elapsed_time;
  printf("Elapsed time = %g seconds.\n", elapsed_time);
  printf("Framerate = %g fps.\n", fps);

  /*-----------------------------------------------------------------------
   *  stop data transmission
   *-----------------------------------------------------------------------*/
  err=dc1394_video_set_transmission(camera,DC1394_OFF);
  DC1394_ERR_RTN(err,"couldn't stop the camera?");

  return (status);
}

max_stars_repo_name: marissakotze/timDIMM

max_issues_repo_name: marissakotze/timDIMM

max_forks_repo_name: marissakotze/timDIMM


text:
// MIT License Copyright (c) 2020 Jarrett Wendt

#pragma once

#include <gsl/gsl>

#include "Macros.h"
#include "SharedPtr.h"
#include "Entity.h"

namespace Library
{
	class Engine final
	{
		friend Entity;

		/** the parentmost Entity */
		static inline SharedPtr<Entity> world{ nullptr };

		static inline std::string programName{};
		static inline std::string pythonSourceDirectory{};
		static inline const std::string initFileName{ "init.py" };
		static inline const std::string shutdownFileName{ "del.py" };

		static inline FILE* initFilePtr{ nullptr };

	public:
		STATIC_CLASS(Engine)
		using Args = gsl::span<const char*>;

		/**
		 * O(1)
		 *
		 * @returns		the parentmost Entity
		 */
		static Entity& World() noexcept;

		/**
		 * to be called by WinMain()
		 */
		static void Main(const Args& args);

		/**
		 * setting this to false will terminate the simulation
		 *
		 * @returns		whether or not the simulation is active
		 */
		static bool& IsActive() noexcept;

		/**
		 * runs before the first Update()
		 */
		static void Init();

		/**
		 * the main engine loop
		 */
		static void Update();

		/**
		 * runs after the last Update()
		 */
		static void Terminate();

	private:
		/**
		 * wraps all of the logic for parsing the command line args.
		 */
		static void ParseArgs(const Args& args);

		/**
		 * Invoked by WinMain's WndProc.
		 *
		 * @param hWnd		HWND
		 * @param uMsg		UINT
		 * @param wParam	WPARAM
		 * @param lParam	LPARAM
		 */
#ifdef _WIN64
		static long long __stdcall WndProc(void* hWnd, unsigned int uMsg, unsigned long long wParam, long long lParam);
#else
		static long __stdcall WndProc(void* hWnd, unsigned int uMsg, unsigned wParam, long lParam);
#endif
	};
}

max_stars_repo_name: JarrettWendt/FIEAEngine

max_issues_repo_name: JarrettWendt/FIEAEngine

max_forks_repo_name: JarrettWendt/FIEAEngine


text:
#ifdef GSLODE1

////////////////////////////// ODE V1 //////////////////////////////////////////

#include <gsl/gsl_odeiv.h>

typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;

int odefunc (double t, const double y[], double f[], void *params) {
    Tode * P = (Tode*) params;
    (P->f)(t,P->n,y,P->n,f);
    return GSL_SUCCESS;
}

int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
     Tode * P = ((Tode*) params);
     (P->j)(t,P->n,y,P->n,P->n,dfdy);
     int j;
     for (j=0; j< P->n; j++)
        dfdt[j] = 0.0;
     return GSL_SUCCESS;
}

int ode(int method, double h, double eps_abs, double eps_rel,
        int f(double, int, const double*, int, double*),
        int jac(double, int, const double*, int, int, double*),
        KRVEC(xi), KRVEC(ts), RMAT(sol)) {

    const gsl_odeiv_step_type * T;

    switch(method) {
        case 0 : {T = gsl_odeiv_step_rk2; break; }
        case 1 : {T = gsl_odeiv_step_rk4; break; }
        case 2 : {T = gsl_odeiv_step_rkf45; break; }
        case 3 : {T = gsl_odeiv_step_rkck; break; }
        case 4 : {T = gsl_odeiv_step_rk8pd; break; }
        case 5 : {T = gsl_odeiv_step_rk2imp; break; }
        case 6 : {T = gsl_odeiv_step_rk4imp; break; }
        case 7 : {T = gsl_odeiv_step_bsimp; break; }
        case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); }
        case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); }
        case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); }
        default: ERROR(BAD_CODE);
    }

    gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
    gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);

    Tode P;
    P.f = f;
    P.j = jac;
    P.n = xin;

    gsl_odeiv_system sys = {odefunc, odejac, xin, &P};

    double t = tsp[0];

    double* y = (double*)calloc(xin,sizeof(double));
    int i,j;
    for(i=0; i< xin; i++) {
        y[i] = xip[i];
        solp[i] = xip[i];
    }

       for (i = 1; i < tsn ; i++)
         {
           double ti = tsp[i];
           while (t < ti)
             {
               gsl_odeiv_evolve_apply (e, c, s,
                                       &sys,
                                       &t, ti, &h,
                                       y);
               // if (h < hmin) h = hmin;
             }
           for(j=0; j<xin; j++) {
               solp[i*xin + j] = y[j];
           }
         }

    free(y);
    gsl_odeiv_evolve_free (e);
    gsl_odeiv_control_free (c);
    gsl_odeiv_step_free (s);
    return 0;
}

#else

///////////////////// ODE V2 ///////////////////////////////////////////////////

#include <gsl/gsl_odeiv2.h>

typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;

int odefunc (double t, const double y[], double f[], void *params) {
    Tode * P = (Tode*) params;
    (P->f)(t,P->n,y,P->n,f);
    return GSL_SUCCESS;
}

int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
     Tode * P = ((Tode*) params);
     (P->j)(t,P->n,y,P->n,P->n,dfdy);
     int j;
     for (j=0; j< P->n; j++)
        dfdt[j] = 0.0;
     return GSL_SUCCESS;
}

int ode(int method, double h, double eps_abs, double eps_rel,
        int f(double, int, const double*, int, double*),
        int jac(double, int, const double*, int, int, double*),
        KRVEC(xi), KRVEC(ts), RMAT(sol)) {

    const gsl_odeiv2_step_type * T;

    switch(method) {
        case 0 : {T = gsl_odeiv2_step_rk2; break; }
        case 1 : {T = gsl_odeiv2_step_rk4; break; }
        case 2 : {T = gsl_odeiv2_step_rkf45; break; }
        case 3 : {T = gsl_odeiv2_step_rkck; break; }
        case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
        case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
        case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
        case 7 : {T = gsl_odeiv2_step_bsimp; break; }
        case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
        case 9 : {T = gsl_odeiv2_step_msadams; break; }
        case 10: {T = gsl_odeiv2_step_msbdf; break; }
        default: ERROR(BAD_CODE);
    }

    Tode P;
    P.f = f;
    P.j = jac;
    P.n = xin;

    gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};

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

    double t = tsp[0];

    double* y = (double*)calloc(xin,sizeof(double));
    int i,j;
    int status;
    for(i=0; i< xin; i++) {
        y[i] = xip[i];
        solp[i] = xip[i];
    }

       for (i = 1; i < tsn ; i++)
         {
           double ti = tsp[i];

           status = gsl_odeiv2_driver_apply (d, &t, ti, y);

           if (status != GSL_SUCCESS) {
         	  printf ("error in ode, return value=%d\n", status);
         	  break;
        	}

//           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);

           for(j=0; j<xin; j++) {
               solp[i*xin + j] = y[j];
           }
         }

    free(y);
    gsl_odeiv2_driver_free (d);

    return status;
}

#endif

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell