Impala Documentation

template<class C>
int Impala::Core::Matrix::Ludcmp	(	C *	a,
		int	n,
		short *	indx,
		double *	d
	)			`[inline]`
Definition at line 21 of file MatInverse.h.
References Matrix_EPS.
Referenced by MatInverse().
00022 {
00023     double *vv = new double [n];
00024 
00025     *d=1.0;
00026     int i, j, k;
00027     for ( i=0;i<n;i++)
00028     {
00029         double big=0.0;
00030         for (j=0;j<n;j++)
00031         {
00032             double temp = fabs(a[i*n+j]);
00033             if (temp > big)
00034                 big=temp;
00035         }
00036         if (big == 0.0)
00037             return(0); // Singular matrix 
00038         vv[i]=1.0/big;
00039     }
00040     for (j=0;j<n;j++)
00041     {
00042         for (i=0;i<j;i++)
00043         {
00044             double sum=a[i*n+j];
00045             for (k=0;k<i;k++)
00046                 sum -= a[i*n+k]*a[k*n+j];
00047             a[i*n+j]=sum;
00048         }
00049         double big=0.0;
00050         int imax = 0;
00051         for (i=j;i<n;i++)
00052         {
00053             double sum=a[i*n+j];
00054             for (k=0;k<j;k++)
00055                 sum -= a[i*n+k]*a[k*n+j];
00056             a[i*n+j]=sum;
00057 
00058             double dum = vv[i]*fabs(sum);
00059             if ( dum >= big)
00060             {
00061                 big=dum;
00062                 imax=i;
00063             }
00064         }
00065         if (j != imax)
00066         {
00067             for (k=0;k<n;k++)
00068             {
00069                 double dum=a[imax*n+k];
00070                 a[imax*n+k]=a[j*n+k];
00071                 a[j*n+k]=dum;
00072             }
00073             *d = -(*d);
00074             vv[imax]=vv[j];
00075         }
00076         indx[j]=imax;
00077         if (a[j*n+j] == 0.0)
00078             a[j*n+j]=Matrix_EPS;
00079         if (j != n-1)
00080         {
00081             double dum=1.0/a[j*n+j];
00082             for (i=j+1;i<n;i++)
00083                 a[i*n+j] *= dum;
00084         }
00085     }
00086     delete [] vv;
00087 
00088     return(1);
00089 }