home *** CD-ROM | disk | FTP | other *** search
/ vis-ftp.cs.umass.edu / vis-ftp.cs.umass.edu.tar / vis-ftp.cs.umass.edu / pub / Software / ASCENDER / ascender.tar.Z / ascender.tar / Triangulate / ludcmp.c < prev    next >
C/C++ Source or Header  |  1995-04-13  |  1KB  |  60 lines
  1. #include <math.h>
  2.  
  3. #define TINY 1.0e-20;
  4.  
  5. void ludcmp(a,n,indx,d)
  6. int n,*indx;
  7. double **a,*d;
  8. {
  9.     int i,imax,j,k;
  10.     double big,dum,sum,temp;
  11.     double *vv,*dvector();
  12.     void nrerror(),free_dvector();
  13.  
  14.     vv=dvector(0,n-1);
  15.     *d=1.0;
  16.     for (i=0;i<n;i++) {
  17.         big=0.0;
  18.         for (j=0;j<n;j++)
  19.             if ((temp=fabs(a[i][j])) > big) big=temp;
  20.         if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
  21.         vv[i]=1.0/big;
  22.     }
  23.     for (j=0;j<n;j++) {
  24.         for (i=0;i<j;i++) {
  25.             sum=a[i][j];
  26.             for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
  27.             a[i][j]=sum;
  28.         }
  29.         big=0.0;
  30.         for (i=j;i<n;i++) {
  31.             sum=a[i][j];
  32.             for (k=0;k<j;k++)
  33.                 sum -= a[i][k]*a[k][j];
  34.             a[i][j]=sum;
  35.             if ( (dum=vv[i]*fabs(sum)) >= big) {
  36.                 big=dum;
  37.                 imax=i;
  38.             }
  39.         }
  40.         if (j != imax) {
  41.             for (k=0;k<n;k++) {
  42.                 dum=a[imax][k];
  43.                 a[imax][k]=a[j][k];
  44.                 a[j][k]=dum;
  45.             }
  46.             *d = -(*d);
  47.             vv[imax]=vv[j];
  48.         }
  49.         indx[j]=imax;
  50.         if (a[j][j] == 0.0) a[j][j]=TINY;
  51.         if (j != n) {
  52.             dum=1.0/(a[j][j]);
  53.             for (i=j+1;i<n;i++) a[i][j] *= dum;
  54.         }
  55.     }
  56.     free_dvector(vv,0,n-1);
  57. }
  58.  
  59. #undef TINY
  60.