/*    ATK IV - Numerical Programming (2006)      */
/*    Excercise 8.1 -  Gaussian Elimination      */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#include<nrutil.h>

void gaussel(float **a, float *b, float *x, int n);
void gausspivot(float **a, float *b, float *x, int n);

int main(void){
  float **a, **c;
  float *b, *x, *e;
  int *indx;
  int n;
  float d;

  FILE *file;
  int i,j;

  file=fopen("matrix_01.dat","rt");
  fscanf(file,"%d\n",&n); 
   
  a=matrix(1,n,1,n+1);
  b=vector(1,n);
  x=vector(1,n);
  indx = ivector(1,n);
  c=matrix(1,n,1,n);
  e=vector(1,n);

  for(i=1;i<=n;i++){
    for(j=1;j<=n;j++){
      fscanf(file,"%f",&a[i][j]);
    }    
  }

  for(i=1;i<=n;i++){
    fscanf(file,"%f",&b[i]);
  }
  fclose(file);
  for(i=1;i<=n;i++){
    e[i]=b[i];
    for(j=1;j<=n;j++){
      c[i][j]=a[i][j];
    }
  }

  gausspivot(a,b,x,n);

  ludcmp(c,n,indx,&d);
  lubksb(c,n,indx,e);

  for(i=1;i<=n;i++){
    printf("%f %f\n",e[i],x[i]);
  }
  
  free_ivector(indx,1,n);
  free_vector(e,1,n);
  free_matrix(c,1,n,1,n);
  free_matrix(a,1,n,1,n);
  free_vector(b,1,n);
  free_vector(x,1,n);
  return 0;
}

void gaussel(float **a, float *b, float *x, int n){
  int i,j,k;

  for(i=1;i<=n;i++){
    a[i][n+1]=b[i];
  }

  for(i=1;i<=n;i++){
    for(j=i+1;j<=n+1;j++){
      a[i][j]=a[i][j]/a[i][i];
    }

    for(j=i+1;j<=n;j++){
      for(k=i+1;k<=n+1;k++){
	a[j][k]=a[j][k]-a[j][i]*a[i][k];
      }
    }
  }

  x[n]=a[n][n+1];

  for(i=n-1;i>=1;i--){
    x[i]=a[i][n+1];
    for(j=i+1;j<=n;j++){
      x[i]=x[i]-a[i][j]*x[j];
    }
  }
}

void gausspivot(float **a, float *b, float *x, int n){
  int i,j,k,ind;
  float dummy, tiny=1.0e-20;

  for(i=1;i<=n;i++){
    a[i][n+1]=b[i];
  }

  for(i=1;i<=n;i++){
    ind=i;
    for(j=i+1;j<=n;j++){
      if(fabs(a[j][i])>fabs(a[ind][i])) ind=j;
    }

    if(ind != i){
      for(k=i;k<=n+1;k++){
	dummy = a[i][k];
	a[i][k] = a[ind][k];
	a[ind][k] = dummy;
      }
    }
  for(j=i+1;j<=n+1;j++){
    a[i][j]=a[i][j]/a[i][i];
  }

  for(j=i+1;j<=n;j++){
    for(k=i+1;k<=n+1;k++){
      a[j][k]=a[j][k]-a[j][i]*a[i][k];
    }
  }
}

  x[n]=a[n][n+1];

  for(i=n-1;i>=1;i--){
    x[i]=a[i][n+1];
    for(j=i+1;j<=n;j++){
      x[i]=x[i]-a[i][j]*x[j];
    }
  }
}