/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*       Solving Nonlinear System of Equations */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-10-04                    */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multiroots.h>

/* Dimension of Matrix and Vectors */
#define DIM 3

/* Limit of Iterative Times */
#define MAXTIMES 100

/* func(x) = 0 */
/* x1 + x2 + x3 - 6 = 0 */
/* x1 * x2 * x3 - 6 = 0 */
/* x1^2 + x2^2 + x3^2 - 14 = 0 */
int func_m(const gsl_vector *x, void *param, gsl_vector *ret)
{
	double x1, x2, x3;

	x1 = gsl_vector_get(x, 0);
	x2 = gsl_vector_get(x, 1);
	x3 = gsl_vector_get(x, 2);

	gsl_vector_set(ret, 0, x1 + x2 + x3 - 6.0);
	gsl_vector_set(ret, 1, x1 * x2 * x3 - 6.0);
	gsl_vector_set(ret, 2, x1*x1 + x2*x2 + x3*x3 - 14.0);

	return GSL_SUCCESS;
}

/* J = [      1,       1,       1 ] */
/*     [x2 * x3, x1 * x3, x1 * x2 ] */
/*     [2 * x1 ,  2 * x2, 2 * x3  ] */
int dfunc_m(const gsl_vector *x, void *param, gsl_matrix *ret)
{
	double x1, x2, x3;

	x1 = gsl_vector_get(x, 0);
	x2 = gsl_vector_get(x, 1);
	x3 = gsl_vector_get(x, 2);

	gsl_matrix_set(ret, 0, 0, 1.0);
	gsl_matrix_set(ret, 0, 1, 1.0);
	gsl_matrix_set(ret, 0, 2, 1.0);

	gsl_matrix_set(ret, 1, 0, x2 * x3);
	gsl_matrix_set(ret, 1, 1, x1 * x3);
	gsl_matrix_set(ret, 1, 2, x1 * x2);

	gsl_matrix_set(ret, 2, 0, 2.0 * x1);
	gsl_matrix_set(ret, 2, 1, 2.0 * x2);
	gsl_matrix_set(ret, 2, 2, 2.0 * x3);

	return GSL_SUCCESS;
}

/* simultaneous eval of f and df */
int fdfunc_m(const gsl_vector *x, void *param, gsl_vector *ret_vec, gsl_matrix *ret_mat)
{
	func_m(x, param, ret_vec);
	dfunc_m(x, param, ret_mat);

	return GSL_SUCCESS;
}

int main(void)
{
	int i, times, status;
	gsl_multiroot_function f;
	gsl_multiroot_fsolver *workspace_f;
	gsl_multiroot_function_fdf fdf;
	gsl_multiroot_fdfsolver *workspace;
	gsl_vector *x;


	/* allocate a, x, b */
	x = gsl_vector_alloc(DIM);

/* Fsolver */

	/* Define Solver */
	workspace_f = gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrid, DIM);

	printf("F solver: %s\n", gsl_multiroot_fsolver_name(workspace_f));

	f.f = &func_m;
	f.n = DIM;
	f.params = 0;

	/* set initial value */
	for(i = 0; i < DIM; i++)
		gsl_vector_set(x, i, 1.0);

	/* set solver */
	gsl_multiroot_fsolver_set(workspace_f, &f, x);

	/* main loop */
	for(times = 0; times < MAXTIMES; times++)
	{
		status = gsl_multiroot_fsolver_iterate(workspace_f);

		printf("%d times: ", times);
		for(i = 0; i < DIM; i++)
			printf("%10.3e ", gsl_vector_get(workspace_f->x, i));
		printf("\n");

		if((status == GSL_EBADFUNC) || (status == GSL_ENOPROG))
		{
			printf("Status: %s\n", gsl_strerror(status));
			break;
		}
	}

	/* print answer */
	for(i = 0; i < DIM; i++)
		printf("%3d %25.17e\n", i, gsl_vector_get(x, i));

	gsl_multiroot_fsolver_free(workspace_f);

/* FDFsolver */

	/* Define Solver */
	workspace = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybridsj, DIM);

	printf("FDF solver: %s\n", gsl_multiroot_fdfsolver_name(workspace));

	fdf.f = &func_m;
	fdf.df = &dfunc_m;
	fdf.fdf = &fdfunc_m;
	fdf.n = DIM;
	fdf.params = 0;

	/* set initial value */
	for(i = 0; i < DIM; i++)
		gsl_vector_set(x, i, 1.0);

	/* set solver */
	gsl_multiroot_fdfsolver_set(workspace, &fdf, x);

	/* main loop */
	for(times = 0; times < MAXTIMES; times++)
	{
		status = gsl_multiroot_fdfsolver_iterate(workspace);

		printf("%d times: ", times);
		for(i = 0; i < DIM; i++)
			printf("%10.3e ", gsl_vector_get(workspace->x, i));
		printf("\n");

		if((status == GSL_EBADFUNC) || (status == GSL_ENOPROG))
		{
			printf("Status: %s\n", gsl_strerror(status));
			break;
		}
	}

	/* print answer */
	for(i = 0; i < DIM; i++)
		printf("%3d %25.17e\n", i, gsl_vector_get(workspace->x, i));

	gsl_multiroot_fdfsolver_free(workspace);

	/* free x */
	gsl_vector_free(x);

	return 0;
}