#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit.h>

int main() {
  gsl_multifit_linear_workspace *work; 
  gsl_matrix *M, *cov;
  gsl_vector *vweight, *vimage, *b;
  double chisq;
  int status;
  size_t i, j, rank;
  FILE *f;

  M = gsl_matrix_alloc(92, 91);
  vweight = gsl_vector_alloc(92);
  vimage = gsl_vector_alloc(92);
  b = gsl_vector_alloc(91);
  cov = gsl_matrix_alloc(91, 91);
  work = gsl_multifit_linear_alloc(M->size1, M->size2);
  printf("Creating random data...\n");
  fflush(stdout);
  for (i=0; i<M->size1; i++) 
    for (j=0; j<M->size2; j++) 
      gsl_matrix_set(M, i, j, drand48());
  for (i=0; i<vweight->size; i++) 
    gsl_vector_set(vweight, i, drand48());
  for (i=0; i<vimage->size; i++) 
    gsl_vector_set(vimage, i, drand48());
  printf("Starting decomposition...\n");
  fflush(stdout);
  status = gsl_multifit_wlinear_svd(M, vweight, vimage, 
				    1e-6, &rank, b, cov,
				    &chisq, work);
  printf("Done.\n");
  fflush(stdout);
  printf("Reading data from file...\n");
  f = fopen("gsl-test.dat", "rb");
  gsl_matrix_fread(f, M);
  gsl_vector_fread(f, vweight);
  gsl_vector_fread(f, vimage);
  printf("Starting decomposition...\n");
  fflush(stdout);
  status = gsl_multifit_wlinear_svd(M, vweight, vimage, 
				    1e-6, &rank, b, cov,
				    &chisq, work);
  printf("Done.\n");
  fflush(stdout);
  
  return 0;
}