/*
  GSL version: 1.9
  
  Compile and run:
   gcc `gsl-config --cflags` svdbug.c `gsl-config --libs` -o svdbug && ./svdbug
 */
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

int main(int argc, char** argv) {
  int m = 462, n = 421;
  gsl_matrix* A = gsl_matrix_alloc(m, n);

  // Read the matrix A from the file svdbug.c, show the first 10 lines and report
  // the number of zeros and ones.
  FILE* fp = fopen("svdbug.csv", "r");
  int i, j, n_ones = 0, n_zeros = 0;
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      double a;
      fscanf(fp, "%lf", &a);
      gsl_matrix_set(A, i, j, a);
      if (i < 10) { printf("%g", a); }
      if (a == 0) { n_zeros++; }
      else if (a == 1) { n_ones++; }
    }
    if (i < 10) { printf("\n"); }
  }
  fclose(fp);
  printf("m = %d, n = %d, %d zeros, %d ones, %d elements\n", m, n, n_zeros,
      n_ones, m * n);
  
  gsl_matrix* V = gsl_matrix_alloc(n, n);
  gsl_vector* S = gsl_vector_alloc(n);
  gsl_vector* work = gsl_vector_alloc(n);
  gsl_linalg_SV_decomp(A, V, S, work);
  printf("first 10 singular values:\n");
  for (i = 0; i < 10; i++) {
    if (i > 0) { printf(", "); }
    printf("%g", gsl_vector_get(S, i));
  }
  return 0;
}