Impala Documentation

template<class C>

static int Impala::Core::Matrix::ludcmp	(	C *	a,
		int	n,
		short *	indx,
		double *	d
	)			`[static]`

Definition at line 489 of file MatrixTem.h.

References HxMatrix_EPS.

Referenced by Impala::Core::Matrix::MatrixTem< C >::i().

00490 {
00491     double *vv = new double [n];
00492 
00493     *d=1.0;
00494         int i, j, k;
00495     for ( i=0;i<n;i++) {
00496         double big=0.0;
00497         for (j=0;j<n;j++) {
00498             double temp = fabs(a[i*n+j]);
00499             if (temp > big)
00500                 big=temp;
00501         }
00502         if (big == 0.0)
00503             return(0); // Singular matrix 
00504         vv[i]=1.0/big;
00505     }
00506     for (j=0;j<n;j++) {
00507         for (i=0;i<j;i++) {
00508             double sum=a[i*n+j];
00509             for (k=0;k<i;k++)
00510                 sum -= a[i*n+k]*a[k*n+j];
00511             a[i*n+j]=sum;
00512         }
00513         double big=0.0;
00514         int imax = 0;
00515         for (i=j;i<n;i++) {
00516             double sum=a[i*n+j];
00517             for (k=0;k<j;k++)
00518                 sum -= a[i*n+k]*a[k*n+j];
00519             a[i*n+j]=sum;
00520 
00521             double dum = vv[i]*fabs(sum);
00522             if ( dum >= big) {
00523                 big=dum;
00524                 imax=i;
00525             }
00526         }
00527         if (j != imax) {
00528             for (k=0;k<n;k++) {
00529                 double dum=a[imax*n+k];
00530                 a[imax*n+k]=a[j*n+k];
00531                 a[j*n+k]=dum;
00532             }
00533             *d = -(*d);
00534             vv[imax]=vv[j];
00535         }
00536         indx[j]=imax;
00537         if (a[j*n+j] == 0.0)
00538             a[j*n+j]=HxMatrix_EPS;
00539         if (j != n-1) {
00540             double dum=1.0/a[j*n+j];
00541             for (i=j+1;i<n;i++)
00542                 a[i*n+j] *= dum;
00543         }
00544     }
00545     delete [] vv;
00546 
00547     return(1);
00548 }