void ludcmp(double a[], int n, int indx[], int d)
{
	const int nmax=100;
	const double tiny=1e-20;
	double vv[100];
	double sum,dum,aamax;
	int i,j,k,imax;
	d=1;
	for (i=1; i<=n; i++)
	{
		aamax=0.0;
		for (j=1; j<=n; j++)
		{
			if (fabs(a[(i-1)*n+j])>aamax)
			{
				aamax=fabs(a[(i-1)*n+j]);
			}
		}
		if (aamax==0.0)
		{
			cout<<"singular matrix"<<endl;
		}
		vv[i]=1.0/aamax;
	}
	for (j=1; j<=n; j++)
	{
        if (j > 1)
		{
            for (i=1; i<=j-1; i++)
			{
                sum=a[(i-1)*n+j];
                if (i>1)
				{
                    for (int k=1; k<=i-1; k++)
					{
                        sum=sum-a[(i-1)*n+k]*a[(k-1)*n+j];
                    }
                    a[(i-1)*n+j]=sum;
                }
            }
        }
        aamax=0.0;
        for (i=j; i<=n; i++)
		{
            sum=a[(i-1)*n+j];
            if (j>1)
			{
                for (k=1; k<=j-1; k++)
				{
                    sum=sum-a[(i-1)*n+k]*a[(k-1)*n+j];
				}
                a[(i-1)*n+j]=sum;
            }
            dum=vv[i]*fabs(sum);
            if (dum>=aamax)
			{
                imax=i;
                aamax=dum;
            }
        }
        if (j!=imax)
		{
            for (k=1; k<=n; k++)
			{
                dum=a[(imax-1)*n+k];
                a[(imax-1)*n+k]=a[(j-1)*n+k];
                a[(j-1)*n+k]=dum;
            }
            d=-d;
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (j!=n)
		{
            if (a[(j-1)*n+j]==0.0)
			{
				a[(j-1)*n+j]=tiny;
			}
            dum=1.0/a[(j-1)*n+j];
            for (i=j+1; i<=n; i++)
			{
                a[(i-1)*n+j]=a[(i-1)*n+j]*dum;
            }
		}
    }
    if (a[(n-1)*n+n]==0.0)
	{
		a[(n-1)*n+n]=tiny;
	}
}