#include <iostream>

#include <stdio.h>

#include <gsl/gsl_multifit.h>

#include <gsl/gsl_matrix.h>

#include <gsl/gsl_vector.h>

using namespace std;

double reg(int n, int p, double *mo, double *pg, double *l, double *ll)

{

double store[p];

gsl_multifit_linear_workspace *ws;

gsl_matrix *cov, *X;

gsl_vector *y, *c;

double chisq;

int i, j;

X = gsl_matrix_alloc(n, p);

y = gsl_vector_alloc(n);

c = gsl_vector_alloc(p);

cov = gsl_matrix_alloc(p, p);

for(i=0; i < n; i++)

{

gsl_matrix_set(X, i, 0, 1.0);

gsl_matrix_set(X, i, 1, mo[i]);

gsl_matrix_set(X, i, 2, pg[i]);

gsl_matrix_set(X, i, 3, l[i]);

gsl_vector_set(y, i, ll[i]);

}

ws = gsl_multifit_linear_alloc(n, p);

gsl_multifit_linear(X, y, c, cov, &chisq, ws);

/* store result ... */

for(i=0; i < p; i++)

{

store[i] = gsl_vector_get(c, i);

}

gsl_multifit_linear_free(ws);

gsl_matrix_free(X);

gsl_matrix_free(cov);

gsl_vector_free(y);

gsl_vector_free(c);

return *store; /* we do not "analyse" the result (cov matrix mainly)

to know if the fit is "good" */

}

int main ()

{

int pixel = 4, parameter=4;

gsl_vector* st;

st = gsl_vector_alloc (parameter);

gsl_vector_free (st);

double modis[]={4.31, 3.83, 3.13, 3.51};

double pg3[]={2.5, 3.94, 3.63, 2.08};

double prel[]={2.50, 4.39, 4.09, 2.95};

double nexl[]={3.4,3.8,3.3,2.7};

double

reg(pixel, parameter, *modis, *pg3, *prel, *nexl);

system("pause");

return 0;

}