/* 
   Matrix Inversion using 
   LU Decomposition from
   Numerical Recipes in C
   Chapter 2
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nrutil.h"

#define TINY 1.0e-20




void ludcmp(float **a, int n, int *indx, float *d)
{
 int i,imax,j,k;
 float   big,dum,sum,temp;
 float   *vv;

 vv = (float*)malloc((unsigned)(n*sizeof(float)));
 if (!vv) 
   {
    fprintf(stderr,"Error Allocating Vector Memory\n");
    exit(1);
   }
 *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)
       {
        fprintf(stderr,"Singular Matrix in Routine LUDCMP\n");
     for (j=0;j<n;j++) printf(" %f ",a[i][j]); printf("/n");
     exit(1);
       }
     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) a[j][j] = TINY;
     if (j != n-1)
       {
     dum = 1.0 / a[j][j];
     for (i=j+1;i<n;i++) a[i][j] *= dum;
       }
    }
 free(vv);
}