#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_linalg.h"

#include <stdio.h>
#define M 5
#define N 3

int main(void)
{
	double A[M][N] = {{1.0, 1.0, 1.0},
			  {1.0, 2.0, 4.0},
			  {1.0, 3.0, 9.0},
			  {1.0, 4.0, 16.0},
			  {1.0, 5.0, 25.0}};
	double b[M] = {1.0, 2.3, 4.6, 3.1, 1.2};
	gsl_vector *vb;
	gsl_matrix *mA;
	int rc;

	vb = gsl_vector_alloc(M);
	mA = gsl_matrix_alloc(M, N);	
	for (int i = 0; i < M; ++i) {
		gsl_vector_set(vb, i, b[i]);
		for (int j = 0; j < N; ++j)
			gsl_matrix_set(mA, i, j, A[i][j]);
	}
	
	rc = gsl_linalg_HH_svx(mA, vb);
	printf("Solution Vector x:");
	for (int i = 0; i < N; ++i)
		printf("\t%0.6f", gsl_vector_get(vb, i));
	printf("\n");
	return rc;
}