Impala Documentation

template<class C>

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

Definition at line 481 of file Matrix.h.

References HxMatrix_EPS.

Referenced by CxMatrixTem< C >::i().

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