home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_07_04 / v7n4062a.txt < prev    next >
Text File  |  1980-07-09  |  4KB  |  108 lines
  1. #define max(i,j)    ((i)>(j)?(i):(j))
  2. #define min(i,j)    ((i)<(j)?(i):(j))
  3. #ifdef __STDC__
  4. #include <math.h>    /* include this if available */
  5. #include <limits.h>
  6. #define FABSMASK    INT_MAX
  7. #else
  8. #define FABSMASK    0x7fffffff
  9. #endif
  10. #ifndef fabsf
  11. #ifndef fabs
  12.     double fabs();
  13. #endif
  14. #define fabsf(x)    fabs((double)x)
  15. #endif
  16. float lufact(a,n,indx)
  17.     float **a; /* a[1..n][1..n] as in Press, Flannery et al */
  18.     int n,indx[]; /* return indx[1..n],factored a[][] */
  19. {
  20.     static float *p,*p2,*p3,accest;
  21.     float *scale=(float *)&indx[1-sizeof(float)/sizeof(int)];
  22. /* put first scale[] then indx[] in caller's array */
  23.     static union intflt{float f;int i;} xmax,xx;
  24. /* smart compilers will use register variables */
  25.     static int i,j,k,ipiv,newmax;
  26.     static double d,atemp,atmp2;
  27.     for(i=1;i<=n;++i){ /* determine scale of each row */
  28.         xmax.f=fabsf(a[i][1]);
  29.         for(j=2;j<=n;++j){
  30. #ifdef ICMP    /* if comparison is to be done by int arithmetic,
  31. **    noting that only + numbers are possible */
  32.             xx.f=a[i][j];
  33. /* consistent style so that this loop will not dominate 
  34. ** time of the other search loop;
  35. ** pointer casting has been tried here,
  36. ** but does not agree with compiler optimizers */
  37.             xmax.i=max(xx.i&FABSMASK,xmax.i);
  38. #else
  39.             xmax.f=max(fabsf(a[i][j]),xmax.f);
  40. #endif
  41.         }
  42. /* except for this assignment,
  43. ** compilers could choose concurrent outer code
  44. ** (assign each loop iteration to independent process) */
  45.         scale[i]=1/xmax.f;
  46.     }
  47.     accest=1; /* initialize to "no loss of accuracy" */
  48.     for(j=2;j<=n;++j){ /* main factorization loop */
  49.         xmax.f=0.; /* search down column for max scaled pivot */
  50.         for(i=j-1;i<=n;++i){
  51. #ifdef ICMP
  52.             xx.f=a[i][j-1]*scale[i];
  53. /* give opportunity for pipelining by using ?
  54. ** instead of if(){} as preferred by certain compilers
  55. **  in 1988 */
  56.             xmax.i=(newmax=(xx.i&=FABSMASK)>xmax.i)?xx.i:xmax.i;
  57. #else
  58.             xmax.f=(newmax=(xx.f=fabsf(a[i][j-1]*scale[i]))>xmax.f)
  59.                 ?xx.f:xmax.f;
  60. #endif
  61.             ipiv=newmax?i:ipiv;
  62.         }
  63.         p=a[ipiv];
  64.         d=p[j-1]; /* this is the pivot element */
  65.         a[ipiv]=a[j-1]; /* swap row pointers */
  66.         a[j-1]=p;
  67.         accest=min(accest,xmax.f); /* track cancellation error */
  68.         scale[ipiv]=scale[j-1];
  69.         indx[j-1] = ipiv; /* record changes for caller */
  70. /* column j and row j-1 calculation combined, unrolled into pairs */
  71.         for(i=j;i<=n;i += 2){
  72.             double sum11=0.,sum12=0.,sum21=0.,sum22=0.;
  73.             p2=a[i];
  74.             p3=a[i+1]; /* may point beyond a[n] */
  75.             for(k=1;k<j-1;++k){
  76. /* use += in inner loop in case sumnn are not allocated to registers */
  77.                 sum11 += a[i][k]*(double)a[k][j];
  78.                 sum12 += a[i+1][k]*(double)a[k][j];
  79.                 sum21 += a[k][i]*(double)a[j-1][k];
  80.                 sum22 += a[k][i+1]*(double)a[j-1][k];
  81.             }
  82. #if defined(__STDC__) && (FLT_ROUNDS>=0) /* rounded normal precision best */
  83.             atemp=p2[j-1]/p[j-1];
  84.             atmp2=p3[j-1]/p[j-1];
  85. #else /* force into double and encourage compiler to preinvert */
  86.             atemp=p2[j-1]/d;
  87.             atmp2=p3[j-1]/d;
  88. #endif
  89. /* reassociation of additions is OK for mixed precision, not for full double */
  90.             p[i] -= sum21; /* completing upper factor */
  91. /* these results will be discarded in the case i==n */
  92.             sum22 = p[i+1]-sum22;
  93. /* note i>=j ; observe sequence */
  94.             sum12 = p3[j]-sum12-atmp2*p[j];
  95. /* storage assignments last to allow more parallel computation */
  96.             p2[j] -= sum11+atemp*p[j];
  97.             p2[j-1]=atemp;
  98.             if(i==n)break;
  99.             p[i+1]=sum22;
  100.             p3[j-1]=atmp2;
  101.             p3[j]=sum12;
  102.         }
  103.     }
  104.     accest = min(fabsf(a[n][n]*scale[n]),accest);
  105.     indx[n]=n;
  106.     return accest;
  107. }
  108.