#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

#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
                        );
                }
        }

        printf( "M =\n" );

        for ( i = 0; i < N; ++i ) {
                for ( j = 0; j < N; ++j ) {
                        printf( "%14.5f ", gsl_matrix_get( M, i, j ) );
                }

                printf( "\n" );
        }

        printf( "\n" );

        /**********************************
         * Find the decomposition M = PLU *
         **********************************/

        gsl_linalg_LU_decomp( M, Pt, &signum );

        printf( "LU =\n" );

        for ( i = 0; i < N; ++i ) {
                for ( j = 0; j < N; ++j ) {
                        printf( "%14.5f ", gsl_matrix_get( M, i, j ) );
                }

                printf( "\n" );
        }

        printf( "\n" );

        printf( "Pt =\n" );

        for ( i = 0; i < N; ++i ) {
                printf( "      %2d\n", gsl_permutation_get( Pt, i ) );
        }

        printf( "\n" );

        /*********************************************
         * Use the PLU decomposition to solve Mx = b *
         *********************************************/

        gsl_linalg_LU_solve( M, Pt, b, x );

        printf( "x =\n" );

        for ( i = 0; i < N; ++i ) {
                printf( "%14.10f\n", gsl_vector_get( x, i ) );
        }

        printf( "\n" );

        /*******************************
         * 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;
}