#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;
}