home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource2 / sclib_1 / 1_4 / v7n4060a.txt < prev    next >
Encoding:
Text File  |  1995-11-01  |  1.8 KB  |  59 lines
  1. main(){ /* test driver for lufact() */
  2. #define NP 6 /* max size for single precision */
  3. #define NP2 NP*2 /* row pointer arrays doubled in size for guard */
  4. float ac[NP][NP],xlc[NP][NP],xuc[NP][NP],xc[NP][NP],
  5.     factr,lufact(),
  6.     *pa[NP2],*pxl[NP2],*pxu[NP2],*px[NP2],
  7. /* shifts to 1 based addressing */
  8.     **a= &pa[-1],**xl= &pxl[-1],**xu= &pxu[-1],**x= &px[-1];
  9. /* make sure indxc has room for float array [NP] */
  10. int indxc[NP2+1],jndxc[NP],*indx= &indxc[-1],*jndx= &jndxc[-1],i,j,k;
  11.     factr=3; /* scale up to make binary value exact */
  12.     for(i=5;i<NP2;i += 2)factr *= i;
  13.     printf("\n Original Scaled Hilbert Matrix");
  14.     for(i=1;i<=NP;++i){ /* set up row pointers base 1 */
  15.         a[i+NP]=a[i]= &ac[i-1][-1]; /* wrap overflow around */
  16.         xl[i+NP]=xl[i]= &xlc[i-1][-1];
  17.         xu[i+NP]=xu[i]= &xuc[i-1][-1];
  18.         x[i+NP]=x[i]= &xc[i-1][-1];
  19.         printf("\n");
  20.         for(j=1;j<=NP;++j)
  21.             printf("%13.3f",a[i][j]=factr/(i+j-1));
  22.     }
  23.     printf("\n\nMatrix Factorization Accuracy: %g",lufact(a,NP,indx));
  24.     printf("\n\n factored matrix");
  25.     for(i=1;i<=NP;++i){ /* copy factors and multiply them as check */
  26.         printf("\n indx[%d]=%d\n",i,indx[i]);
  27.         for(j=1;j<=NP;++j){
  28.             printf("%13g",ac[i-1][j-1]);
  29.             if(j>=i){
  30.                 xu[i][j]=a[i][j];
  31.                 xl[i][j]=(j==i);
  32.             }
  33.             else{
  34.                 xu[i][j]=0;
  35.                 xl[i][j]=a[i][j];
  36.             }
  37.         }
  38.     }
  39.     for(i=1;i<=NP;++i){
  40.         jndx[i]=i;
  41.         for(j=1;j<=NP;++j){
  42.             double sum=0;
  43.             for(k=1;k<=NP;++k)sum += xl[i][k]*xu[k][j];
  44.             x[i][j]=sum;
  45.         }
  46.     }
  47.     printf("\n\n Product of lower and upper matrices");
  48.     for(i=1;i<=NP;++i){
  49.         register int dum=jndx[indx[i]];
  50.         jndx[indx[i]]=jndx[i];
  51.         jndx[i]=dum;
  52.     }
  53.     for(i=1;i<=NP;++i){
  54.         for(j=1;jndx[j]!=i;++j);
  55.         printf("\n");
  56.         for(k=1;k<=NP;++k)printf("%13.6f",x[j][k]);
  57.     }
  58. }
  59.