#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>

int  main (int argc, char *argv[])    {
     int i,j,ndim;
     FILE *fp;
     float *X;
      
     fp=fopen(argv[1],"r");
     fscanf(fp,"%d",&ndim);
     X=(float *)malloc(ndim*ndim*sizeof(float)); 

     for(i=0; i < ndim ; i++){ 
       for(j=0; j < ndim; j++){
         fscanf(fp,"%f",&X[j+ndim*i]); 
      }
      fscanf(fp,"\n"); 
     }
     fclose(fp); 


     gsl_vector *S  = gsl_vector_alloc (ndim);
     gsl_vector *work  = gsl_vector_alloc (ndim);

     gsl_matrix *A  = gsl_matrix_alloc (ndim, ndim);
     gsl_matrix *V  = gsl_matrix_alloc (ndim, ndim);

    for(i=0; i < ndim; i++){
        for(j=0; j < ndim; j++)
           gsl_matrix_set (A, i, j, X[j+ndim*i]);
    }

     gsl_linalg_SV_decomp (A,V,S,work);


    printf("U\n");
    for(i=0; i < ndim ; i++){
       for(j=0; j < ndim; j++){
         printf("%6.2f",gsl_matrix_get (A, i, j));
      }
      printf("\n");
     }

   printf("V\n");
    for(i=0; i < ndim ; i++){
       for(j=0; j < ndim; j++){
         printf("%6.2f",gsl_matrix_get (V, i, j));
      }
      printf("\n");
     }

   printf("S\n");
   for(j=0; j < ndim; j++)
         printf("   %6.2f   ",gsl_vector_get (S,j));
   printf("\n");


     return 0;
 }