home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 188_01 / eqsym.c < prev    next >
Text File  |  1987-09-30  |  3KB  |  102 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        Symmetric matrix factorization
  4. VERSION:    1.0;
  5.  
  6. DESCRIPTION║ááá    Factorization and linear equation solution withì
  7.                full or variable banded symmetric matrix
  8.  
  9. KEYWORDS:    Symmetric matrix linear equation solution
  10. SYSTEM:        CP/M-80, V3.1;
  11. FILENAME:    EQSYM.C;
  12. WARNINGS║áá    requires floating pointì
  13.  
  14. AUTHORS:    Tim Prince;
  15. COMPILERS:    MIX C, v2.0.1; Eco C 3.45; Aztec C 1.06B
  16. */
  17. typedef double OPTREAL; /* can work on float or double arrays*/
  18. main(){/* sample test of simultaneous linear solution */
  19. #define ndef 12 /* 12 for double, 6 for single precision test*/
  20. #define abs(d) (d>0?d:-d)
  21. #define max(d,a) d>a?d:a
  22. OPTREAL a[ndef*(ndef+1)],b[ndef],factr;
  23. double sum,d;
  24. int maxdvr,i,j,id[ndef],ic;
  25. factr=3;
  26. maxdvr=ndef*2-1;
  27. for(i=2;i<=maxdvr;i++)
  28.   factr*=i;/* maxdvr! */
  29. /* scale factor 1 for true Hilbert matrix but this allows
  30. ** exact integral data and is a more severe test */
  31. ic=-1;
  32. for(i=0;i<ndef;i++)
  33.   {/* set up to test accuracy of solution for all 1's */
  34. /* multiple assignment in for doesn't work in Manx Aztec C */
  35.   for(sum=j=0;j<ndef;j++){
  36.     /* set up Hilbert matrix, symmetric variable bandwidth */
  37.     d=factr/(i+j+1);
  38.     if(j<=i){ /* columns are full length */
  39.       ic++;
  40.       a[ic]=d;}
  41.     sum+=d;}
  42.   id[i]=ic;
  43.   b[i]=sum;}; /* [i++] fails on several compilers */
  44. symfac(a,id,ndef);/* replace a by factors*/
  45. symfwb(a,id,ndef,b);/* solve b system */
  46. for(sum=i=0;i<ndef;i++)
  47.   sum=max(abs(b[i]-1),sum);
  48. printf("\n solution error:%g",sum);
  49. }
  50. symfac(a,id,n)
  51. OPTREAL a[];
  52. int n,id[];
  53. {
  54. register double u,sum; /* long double preferred */
  55. register int i,k,ic,iu,il;
  56. int j,k1,iu1;
  57. /* n: order of matrix
  58. a[]: matrix to be overwritten by its triangular factors
  59. a[id[]]: diagonal elements: end of column vector
  60. matrix storage: variable length columns (symmetric)
  61. for(j=1;j<n;j++){ /* symmetric factors L D Lt */
  62. /* start by calculating U=D Lt from 2nd row of column j */
  63.   for(i=k1=(ic=(iu1=id[j-1]+1)+1)-id[j]+j;i<j;ic++,i++){
  64.    /* set pointers to start dot product L row i, U col j
  65.    ** first assuming row i equal or longer than col j
  66.    ** then correcting if "reentrant" */
  67.     if((k=(il=(iu=iu1)-ic+id[i])-id[i-1])<=0){
  68.       iu+=k=1-k; /* skip over implicit zeros in col i */
  69.       il+=k;}
  70.     for(sum=0;iu<ic;iu++,il++)
  71.       sum += a[iu] * a[il];
  72.     a[ic]-=sum;};/* replace by factor U */
  73.   for(k=k1-1,iu=iu1,sum=0;k<j;iu++,k++){
  74.     a[iu]=(u=a[iu])/(a[id[k]]); /* replace with L = Dinv Ut */
  75.     sum += a[iu] * u ;} /* Lt U */
  76.   a[ic] -= sum; /* replace with D */
  77. }
  78. }
  79. symfwb(a,id,n,b)
  80. int n,id[];
  81. OPTREAL a[],b[];
  82. {
  83. register double sum;
  84. register int i,j,ic,il;
  85. /* b[n]: right side vector to be overwritten by
  86.   solution 
  87. solve L D Lt b" = b
  88. writing b' and b" over b 
  89. L[j][j] =1 not stored */
  90. for(j=1;j<n;j++){ /* b' = Linv b */
  91. /* multiply by L inverse */
  92.   for(sum=0,i=j+(il=id[j-1]+1)-id[j];i<j;il++,i++)
  93.       sum += a[il] * b[i];
  94.     b[j] -= sum;};/* L inv b */
  95. for(j=0;j<n;j++) /* b" = Dinv b' */
  96.   b[j] /= a[id[j]];
  97. /* multiply by Lt inverse */
  98. for(j=n-1;j;j--)
  99.   for(i=j-1,il=id[j]-1;il>id[j-1];il--,i--)
  100.     b[i] -= b[j] * a[il];
  101. }
  102.