#include <math.h>
#include "proto.h"

static void pivot(double **a,double *b,int n,int k)
{
  int i,mi;
  double *tmp,t,max;;

  max=fabs(a[k][k]);
  mi=k;
  for (i=k+1; i<n; i++)
    if (fabs(a[i][k]) > max ) { 
      max=fabs(a[i][k]);
      mi=i;
    }
  if (max == 0.0) 
    error_exit("Matrix singular",__FILE__,__LINE__);

  if (mi != k) {
    tmp=a[k];a[k]=a[mi];a[mi]=tmp;
    t=b[k];b[k]=b[mi];b[mi]=t;
  }
  return;
}

double *gauss(double **a,double *b,int n)
{
  int i,k;

  for(k=0; k<n-1; k++) {
    pivot(a,b,n,k);
    scale(a[k]+k+1,1.0/a[k][k],n-(k+1));
    b[k] /= a[k][k];
    for(i=k+1; i<n; i++) {
      xpay(a[i]+k+1,-a[i][k],a[k]+k+1,n-(k+1));
      b[i] -= a[i][k]*b[k];   
    }
  }
  if ( a[n-1][n-1] == 0.0) 
    error_exit("Matrix singular",__FILE__,__LINE__);
  b[n-1] /= a[n-1][n-1];

  for(i=n-2; i>=0; i--)
    b[i] -= xyT(a[i]+i+1,b+i+1,n-(i+1));

  return b;

}