/* run: gcc -o eigtridi_test eigtridi_test.c eigtridi.c -lm */

#include "eigtridi.h"

int main (int argc, char* argv[])
{
  const int size=10;
  double diag[size];
  double sub[size];
  double ev[2*size];

  int i;
  int ierr;

  /* note: 
   * the eigenvalues of this matrix are the nodes for Gauss-Laguerre
   * Interpolation. That means: 
   * If you want to compute \int_0^{\infty} f(x) dx, you can use
   * \sum_{i=0}^n w_i f(x_i) exp(x_i), where x_i are the eigenvalues of
   * the example matrix an w_i are the squares of the first components 
   * of the computed eigenvectors
   */

  for ( i=0; i<size; ++i) {
    diag[i] = 2*i+1;
    sub[i]  = i;
  }
  for ( i=0; i<2*size; ++i) {
    ev[i]   = 0;
  }
  ev[0]=1;
  ev[size+1]=1;
  
  printf("INPUT\n\n");
  printf("diag: ");
  for ( i=0; i<size; ++i) {
    printf("%g, ", diag[i]);
  }
  printf("\nsub: ");
  for ( i=0; i<size; ++i) {
    printf("%g, ", sub[i]);
  }
  printf("\nev:  ");
  for ( i=0; i<2*size; ++i) {
    printf("%g, ", ev[i]);
  }
  printf("\n");

  tridi_eigenvalues(2, size, diag, sub, ev, 2, &ierr);

  printf("OUTPUT\n\n");
  printf("diag: ");
  for ( i=0; i<size; ++i) {
    printf("%g, ", diag[i]);
  }
  printf("\nsub: ");
  for ( i=0; i<size; ++i) {
    printf("%g, ", sub[i]);
  }
  printf("\nev:  ");
  for ( i=0; i<2*size; ++i) {
    printf("%g, ", ev[i]);
  }
  printf("\n");
  
}