home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The C Users' Group Library 1994 August
/
wc-cdrom-cusersgrouplibrary-1994-08.iso
/
listings
/
v_07_04
/
v7n4060a.txt
< prev
next >
Wrap
Text File
|
1980-07-09
|
2KB
|
59 lines
main(){ /* test driver for lufact() */
#define NP 6 /* max size for single precision */
#define NP2 NP*2 /* row pointer arrays doubled in size for guard */
float ac[NP][NP],xlc[NP][NP],xuc[NP][NP],xc[NP][NP],
factr,lufact(),
*pa[NP2],*pxl[NP2],*pxu[NP2],*px[NP2],
/* shifts to 1 based addressing */
**a= &pa[-1],**xl= &pxl[-1],**xu= &pxu[-1],**x= &px[-1];
/* make sure indxc has room for float array [NP] */
int indxc[NP2+1],jndxc[NP],*indx= &indxc[-1],*jndx= &jndxc[-1],i,j,k;
factr=3; /* scale up to make binary value exact */
for(i=5;i<NP2;i += 2)factr *= i;
printf("\n Original Scaled Hilbert Matrix");
for(i=1;i<=NP;++i){ /* set up row pointers base 1 */
a[i+NP]=a[i]= &ac[i-1][-1]; /* wrap overflow around */
xl[i+NP]=xl[i]= &xlc[i-1][-1];
xu[i+NP]=xu[i]= &xuc[i-1][-1];
x[i+NP]=x[i]= &xc[i-1][-1];
printf("\n");
for(j=1;j<=NP;++j)
printf("%13.3f",a[i][j]=factr/(i+j-1));
}
printf("\n\nMatrix Factorization Accuracy: %g",lufact(a,NP,indx));
printf("\n\n factored matrix");
for(i=1;i<=NP;++i){ /* copy factors and multiply them as check */
printf("\n indx[%d]=%d\n",i,indx[i]);
for(j=1;j<=NP;++j){
printf("%13g",ac[i-1][j-1]);
if(j>=i){
xu[i][j]=a[i][j];
xl[i][j]=(j==i);
}
else{
xu[i][j]=0;
xl[i][j]=a[i][j];
}
}
}
for(i=1;i<=NP;++i){
jndx[i]=i;
for(j=1;j<=NP;++j){
double sum=0;
for(k=1;k<=NP;++k)sum += xl[i][k]*xu[k][j];
x[i][j]=sum;
}
}
printf("\n\n Product of lower and upper matrices");
for(i=1;i<=NP;++i){
register int dum=jndx[indx[i]];
jndx[indx[i]]=jndx[i];
jndx[i]=dum;
}
for(i=1;i<=NP;++i){
for(j=1;jndx[j]!=i;++j);
printf("\n");
for(k=1;k<=NP;++k)printf("%13.6f",x[j][k]);
}
}