void lubksb(double a[], int n, int indx[], double b[])
{
	int i,j,ll;
	double sum;
    int ii = 0;
    for (i = 1; i<=n; i++)
	{
        ll = indx[i];
        sum = b[ll];
        b[ll] = b[i];
        if (ii != 0 )
		{
            for (j = ii; j<=i-1; j++)
			{
                sum = sum - a[(i-1)*n+j] * b[j];
            }
		}
        else
		{
			if (sum != 0.0)
			{
				ii = i;
			}
        }
        b[i] = sum;
    }
    for (i = n; i>=1; i--)
	{
        sum = b[i];
        if (i < n)
		{
            for (j = i + 1; j<=n; j++)
			{
                sum = sum - a[(i-1)*n+j] * b[j];
            }
        }
        b[i] = sum / a[(i-1)*n+i];
    }
}