home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 188_01 / eqd.c < prev    next >
Text File  |  1987-09-30  |  6KB  |  187 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        general matrix factorization;
  4. VERSION:    1.0;
  5.  
  6. DESCRIPTION║ááá    Sourcσááácodσááá fo≥áá general full matrix ì
  7.                matrix factorization and linear system solution
  8.  
  9.         Employ≤ Crou⌠ scaleΣ partia∞ pivotinτ fo≥ ì
  10.         ful∞ matrix«  Solutioε oµ linea≥ equation≤ ì
  11.         anΣ accurac∙ test≤ arσ demonstrated.
  12.  
  13. KEYWORDS:    Matrix linear equation solution;
  14. SYSTEM:        CP/M-80, V3.1;
  15. FILENAME:    EQD.C;
  16. WARNINGS║    require≤ compile≥ whicΦ support≤ doublσ ì
  17.         subscript≤ [][▌ anΣ floatinτ point.
  18.  
  19. AUTHORS:    Tim Prince;
  20. COMPILERS:    MIX C, v2.0.1; Aztec C 1.06B
  21. */
  22. typedef double OPTREAL; /* can work on float or double arrays*/
  23. main(){/* sample test of simultaneous linear solution */
  24. #define ndef 12 /* 12 for double, 6 for single precision test*/
  25. #define abs(d) (d>0?d:-d)
  26. #define max(d,a) d>a?d:a
  27. float scale[ndef],accest;
  28. OPTREAL a[ndef][ndef],b[ndef],aa[ndef][ndef],factr;
  29. double sum,det();
  30. char iperm[ndef];
  31. int maxdvr,i,j;
  32. factr=3;
  33. maxdvr=ndef*2-1;
  34. for(i=2;i<=maxdvr;i++)
  35.   factr*=i;/* maxdvr! */
  36. /* scale factor 1 for true Hilbert matrix but this allows
  37. ** exact integral data and is a more severe test */
  38. for(i=0;i<ndef;i++)
  39.   {/* set up to test accuracy of solution for all 1's */
  40.   for(sum=j=0;j<ndef;j++)
  41.     /* set up Hilbert matrix */
  42.     sum+=aa[j][i]=a[j][i]=factr/(i+j+1);
  43.   b[i]=sum;}; /* [i++] fails on several compilers */
  44. factor(a,iperm,scale,&accest);/* replace a by factors*/
  45. printf("\n accest =%e",accest); /* %g doesn't work in Mix C */
  46. #define DETRMNT 1 
  47. #ifdef DETRMNT
  48. sum=det(a,iperm,&i);
  49. printf("\n determinant = ldexp(%g,%d)",sum,i);
  50. #endif
  51. fwdbak(a,iperm,b,1);/* solve b system */
  52. for(sum=i=0;i<ndef;i++)
  53.   sum=max(abs(b[i]-1),sum);
  54. printf("\n solution error:%g",sum);
  55. }
  56. factor(a,iperm,scale,accest)
  57. OPTREAL a[ndef][ndef];
  58. char iperm[ndef];
  59. float scale[ndef],*accest;/* need not be accurate */
  60. {register double d;
  61. register float x,xmax,scali;/* need not be accurate */
  62. register double sum; /* long double preferred */
  63. register int i,j,k,ipiv,l;
  64. /* n: order of matrix
  65. a[ndef][ndef]: matrix to be overwritten by its triangular factors
  66. iperm: record of row exchanges
  67. accest: estimated relative accuracy of factorization
  68. accest =0. if singular */
  69. #define min(i,j) (i<j?i:j)
  70. /* matrix storage: (row,column) stored by columns in FORTRAN
  71. determine scale of each row */
  72. for(i=0;i<ndef;i++)
  73.   {
  74.   for(scali=j=0;j<ndef;j++)
  75.     scali=max(abs(a[j][i]),scali);
  76.   scale[i]=1/scali;};/* avoid more slow divides
  77. now factor by method in Jennings "Matrix Computation for
  78. #Engineers and Scientists" p. 114 */
  79. *accest=1;/* assume exact input*/
  80. for(j=1;j<ndef;j++)
  81.   {
  82.   xmax=0;
  83. /*  look for pivot */
  84.   for(i=l=j-1;i<ndef;i++)
  85.     {if((x=abs(a[l][i])*scale[i])>xmax){
  86. /* find max scaled pivot */
  87.       xmax=x;
  88.       ipiv=i;};};
  89.   /* move row ipiv to row j-1 */
  90.   if(xmax<*accest)*accest=xmax;
  91. /* indicates relative accuracy i.e. .1 for loss of 1
  92. ** digit due to cancellations */
  93.   if(l!=(iperm[l]=ipiv)){/* swap rows; test to save time */
  94.     x=scale[l];
  95.     scale[l]=scale[ipiv];
  96.     scale[ipiv]=x;
  97.     for(k=0;k<ndef;k++)
  98.       {d=a[k][l];
  99.       a[k][l]=a[k][ipiv];
  100.       a[k][ipiv]=d;};};
  101.   d=a[l][l]; /* d=1/a[l][l] OK for float OPTREAL or long
  102. ** double d */
  103.   for(i=1;i<ndef;i++)
  104.     { /* can split in 2 loops and eliminate explicit if else 
  105. ** the i<j part doesn't depend on the swapping code and may be
  106. ** executed in parallel.  In the i>=j part the i iterations are
  107. ** not "vector dependent" and may be executed in parallel. */
  108.     if(i>=j){
  109.       a[j-1][i]/=d;/* complete lower factor */
  110.       l=j;}
  111.     else l=i;
  112.     for(sum=k=0;k<l;k++)
  113.       sum+=a[k][i]*a[j][k];
  114.     a[j][i]-=sum;};/* replace a by factor */
  115. };
  116. *accest=min(abs(a[ndef-1][ndef-1]*scale[ndef-1]),*accest);
  117. }
  118. #ifdef DETRMNT
  119. double det(a,iperm,exp)
  120. OPTREAL a[ndef][ndef];
  121. char iperm[ndef];
  122. int *exp;
  123. { double frexp(),fract;
  124.   int i,pwr;
  125. /* calculate determinant as fraction * 2^exp to keep in range*/
  126.   fract=a[0][0];
  127.   *exp=0;
  128.   for(i=1;i<ndef;i++){
  129.     fract=frexp(fract*a[i][i],&pwr);
  130.     *exp+=pwr;
  131. /* each swap changes sign of determinant */
  132.     if(iperm[i-1]!=i-1)fract=-fract;}
  133.   return(fract);
  134. }
  135.   #ifndef AZTEC
  136. /* dblfmt describes floating point format
  137. ** machine dependencies in case frexp() and ldexp() are not
  138. ** supplied */
  139. struct dblfmt{
  140. /* MIX C like Microsoft format */
  141.   char mant[7]; /* full reverse byte order */
  142.   char expn;
  143. };
  144. double frexp(x,scale)
  145.   int *scale;
  146.   double x;
  147. {
  148.   #define BIAS 128 /* IEEE 1023 */
  149.   *scale=((struct dblfmt *)&x)->expn-BIAS;
  150.   ((struct dblfmt *)&x)->expn=BIAS;
  151.   return(x);
  152. }
  153.   #endif
  154. #endif
  155. fwdbak(a,iperm,b,m)
  156. int m;
  157. OPTREAL a[ndef][ndef],b[][ndef];
  158. char iperm[ndef];
  159. {register double x;
  160. register double sum;
  161. register int i,j,k,ipiv;
  162. /* b[m][ndef]: matrix of right side vectors to be overwritten by
  163.   solution 
  164. solve L U b" = b
  165. writing b' and b" over b 
  166. L stored in a[j][i]: i>j; L[j][j] =1 not stored */
  167. iperm[ndef-1]=ndef-1;
  168. for(k=0;k<m;k++)/* operate by columns, for efficiency if m=1 */
  169.   {
  170.   for(j=0;j<ndef;j++)
  171.     {
  172.     /* exchange as was done in factor */
  173.     x=b[k][ipiv=iperm[j]];
  174. /* eliminating x by setting sum here OK for float OPTREAL */
  175.     b[k][ipiv]=b[k][j];
  176. /* multiply by L inverse */
  177.     for(sum=i=0;i<j;i++)
  178.       sum+=a[i][j]*b[k][i];
  179.     b[k][j]=x-sum;};/* L inv b */
  180. /* multiply by U inverse */
  181.   for(j=ndef-1;j>=0;j--){
  182.     sum=0;
  183.     for(i=j+1;i<ndef;i++)
  184.     sum+=a[i][j]*b[k][i];
  185.     b[k][j]=(b[k][j]-sum)/a[j][j];};};
  186. }
  187.