#include "find_interpolation_parameters.h"

int find_interpolation_parameters_SVD( const size_t dim, 
  const size_t num,
  double * point_data, 
  double * value_data,
  double * parameters )
{
  gsl_matrix * AU   = gsl_matrix_alloc( num,     dim + 1 );
  gsl_matrix * V    = gsl_matrix_alloc( dim + 1, dim + 1 );
  gsl_vector * S    = gsl_vector_alloc( dim + 1 );
  gsl_vector * work = gsl_vector_alloc( dim + 1 );
  gsl_vector * b    = gsl_vector_alloc( num );
  gsl_vector * x    = gsl_vector_alloc( dim + 1 );

  size_t i;
  size_t j;

  for ( i = 0; i < num; ++i )
    {
    for ( j = 0; j < dim; ++j )
      {
      gsl_matrix_set( AU, i, j, point_data[ i * dim + j ] );
      }
    gsl_matrix_set( AU, i, dim, 1.0 );
    }

  for ( i = 0; i < num; ++i )
    {
    gsl_vector_set( b, i, value_data[i] );
    }

  gsl_linalg_SV_decomp( AU, V, S, work );
  gsl_linalg_SV_solve(  AU, V, S, b, x );

  for ( i = 0; i < dim + 1; ++i )
    {
    parameters[i] = gsl_vector_get( x, i );
    }

  gsl_matrix_free( AU );
  gsl_matrix_free( V );
  gsl_vector_free( S );
  gsl_vector_free( work );
  gsl_vector_free( b );
  gsl_vector_free( x );

  return EXIT_SUCCESS;
}