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


/* Program to solve 
 * 
 * | 1 4 0 0 | | x1 |    | 8 |
 * | 3 4 1 0 | | x2 | =  | 4 |
 * | 0 2 3 4 | | x3 |    | 2 |
 * | 0 0 1 3 | | x4 |    | 4 |
 * 
 * using GNU Scientific library. To compile on cobra do 
 *
 * gcc tridiag.c -lgsl -lgslcblas
 *
 */


int main(int argc, char** argv)
  {
  
  double diag_data[] = { 1.0, 4.0, 3.0, 3.0 };
  double sup_data[] = {4.0, 1.0, 4.0};
  double sub_data[] = {3.0, 2.0, 1.0};

  double b_data[] = { 8.0, 4.0, 2.0, 4.0 };
 
  gsl_vector *x = gsl_vector_alloc(4);

  gsl_vector_view diag = gsl_vector_view_array(diag_data, 4); 
  gsl_vector_view sup  = gsl_vector_view_array(sup_data, 3); 
  gsl_vector_view sub  = gsl_vector_view_array(sub_data, 3); 
  gsl_vector_view b    = gsl_vector_view_array(b_data, 4); 

  gsl_linalg_solve_tridiag( &diag.vector, &sup.vector, &sub.vector, &b.vector, x );  

  printf("x = \n");
  gsl_vector_fprintf( stdout, x, "%g");


  gsl_vector_free(x);
  return 0;
}