#include<math.h>
#define TINY 1e-20;

void ludcmp(float **a, int n, int *indx, float *d, float *vv)
{
  int i, imax, j, k;
  float big, dum, sum, temp;
  
  *d = 1.0;
 
 for(i=0; i<=n; i++)
   {
     big = 0.0;
     for (j=0; j<=n; j++)
       if(( temp = fabs(a[i][j]))>big) big=temp;
     if(big == 0.0) printf("\n Singular matrix in routine ludcmp");
     vv[i] = 1.0/big;
   }
 
 for(j=0; j<=n; j++)
   {
     for(i=0; i<j; i++)
       {
	 sum = a[i][j];
	 for(k=0; k<i; k++)
	   sum -= a[i][k]*a[k][j];
	 a[i][j] = sum;
       }
     big = 0.0;
     for(i=j;i<=n; i++)
       {
	 sum = a[i][j];
	 for(k=0; k<j; k++)
	   sum -= a[i][k]*a[k][j];
	 a[i][j] = sum;
	 if((dum = vv[i]*fabs(sum)) >= big)
	   {
	     big = dum;
	     imax = i;
	   }
       }
     
     if(j != imax)
       {
	 for(k=0; k<=n; k++)
	   {
	     dum = a[imax][k];
	     a[imax][k] = a[j][k];
	     a[j][k]=dum;
	   }
	 *d = -(*d);
	 vv[imax] = vv[j];
       }
     indx[j] = imax;
     if(a[j][j] == 0.0)
       {
	 printf("zero a[j][j]\n");
	 a[j][j] = TINY;
       }
     if(j!=n)
       {
	 dum = 1.0/a[j][j];
	 for(i=j+1; i<=n; i++)
	   a[i][j] *= dum;
       }
   }
}

void lubksb(float **a, int n, int *indx, float *b, float *x, float *vv)
{
  int i, ii = -1, ip, j;
  float sum;

  for (i=0;i<=n;i++)
    vv[i] = -b[i];
  
  for (i=0; i<=n; i++)
    {
      ip = indx[i];
      sum = vv[ip];
      vv[ip] = vv[i];
      if(ii>-1)
      for (j=ii; j<=i-1; j++) sum -= a[i][j]*vv[j];
      else if (sum) ii=i;
      vv[i] = sum;
    }
  for(i=n; i>=0; i--)
    {
      sum = vv[i];
      for(j=i+1; j<=n; j++)
	sum -= a[i][j]*vv[j];
      x[i] = vv[i] = sum/a[i][i];
    }
}