#include #include #include #include #include #include #define N 7 /************************************************* * Find the coefficients of the polynomial * interpolating the N points * (1, sin(1)), (2, sin(2)), ..., (N, sin(N)) * The equivalent Matlab code is format long N = 7; M = vander( 1:N ); b = sin( 1:N )'; [L U Pt] = lu( M ) y = L \ (Pt*b); x = U \ y * You will notice that the * permutation matrix differs in * the last two rows due to the * application of scaling * * Acknowledgment: T. Kouya, * http://na-inet.jp/na/gslsample/linear_system.c *************************************************/ int main() { int i, j, signum; gsl_matrix *M; gsl_vector *b, *x; gsl_permutation *Pt; /**************************** * Allocate memory * * - This is equivalent to * * calling 'new' in C++ * ****************************/ M = gsl_matrix_alloc( N, N ); b = gsl_vector_alloc( N ); x = gsl_vector_alloc( N ); Pt = gsl_permutation_alloc( N ); /************************************* * Create the vector T * * b = (sin(1) sin(2) ... sin(N)) * *************************************/ for ( i = 0; i < N; ++i ) { gsl_vector_set( b, i, sin( (double)(i + 1) ) ); } /************************************************ * Create the Vandermonde matrix for the points * * T * * x = (1 2 3 ... N) * ************************************************/ for ( i = 0; i < N; ++i ) { gsl_matrix_set( M, i, N - 1, 1.0 ); } for ( j = N - 2; j >= 0; --j ) { for ( i = 0; i < N; ++i ) { gsl_matrix_set( M, i, j, (double)(i + 1)*gsl_matrix_get( M, i, j ); } } /********************************************* * Find the decomposition M = PLU * * Use the PLU decomposition to solve Mx = b * *********************************************/ gsl_linalg_LU_decomp( M, Pt, &signum ); gsl_linalg_LU_solve( M, Pt, b, x ); /******************************* * Free allocated memory * * - This is equivalent to * * calling 'delete' in C++ * *******************************/ gsl_permutation_free( Pt ); gsl_matrix_free( M ); gsl_vector_free( b ); gsl_vector_free( x ); return 0; }