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

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



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

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