home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 188_01 / eqvb.c < prev    next >
Text File  |  1987-09-30  |  4KB  |  116 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        Unsymmetric variable band matrix factorization
  4. VERSION:    1.0;
  5.  
  6. DESCRIPTION║ááá    Factorization and linear equation solution withì
  7.                variable banded unsymmetric matrix
  8.  
  9. KEYWORDS:    Unymmetric variable band matrix solution
  10. SYSTEM:        CP/M-80, V3.1;
  11. FILENAME:    EQVB.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 au[ndef*(ndef-1)],al[ndef*(ndef+1)],b[ndef],factr;
  23. double sum,d;
  24. int maxdvr,i,j,iu[ndef],il[ndef],icu,icl;
  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. icu=icl=-1;
  32. for(i=0;i<ndef;i++)
  33.   {/* set up to test accuracy of solution for all 1's */
  34.   for(sum=j=0;j<ndef;j++){
  35.     /* set up Hilbert matrix, symmetric variable bandwidth */
  36.     d=factr/(i+j+1);
  37.     /* rows & columns are full length in this example
  38. ** upper triangle stored by columns, lower triangle
  39. ** by rows including diagonal element */
  40.     if(j<i){
  41.       icu++;
  42.       au[icu]=d;} /* a[i][j] */
  43.     if(j<=i){
  44.       icl++;
  45.       al[icl]=d;} /* a[j][i] */
  46.     sum+=d;}
  47.   iu[i]=icu; /* index to element above diagonal (may not exist) */
  48.   il[i]=icl; /* index to diagonal element */
  49.   b[i]=sum;}; /* [i++] fails on several compilers */
  50. vbfac(au,iu,al,il,ndef);/* replace a by factors*/
  51. vbfwb(au,iu,al,il,ndef,b);/* 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. vbfac(au,iu,al,il,n)
  57. OPTREAL au[],al[];
  58. int n,iu[],il[];
  59. {
  60. register double sum; /* long double preferred */
  61. register int i,k,ic,iuc,ilr;
  62. int j,iu1;
  63. /* n: order of matrix
  64. a[]: matrix to be overwritten by its triangular factors
  65. optionally split into au (upper triangle)
  66. and al (lower triangle including main diagonal) for clarity
  67. al[il]: diagonal elements: end of row vector
  68. au[iu]: element of column just above al[il]: end of column
  69. matrix storage: variable length row & column (unsymmetric) */
  70. for(j=1;j<n;j++){ /* factors L U 
  71. ** start by calculating U from 1st row of column j */
  72.   for(i=(ic=iu1=iu[j-1]+1)-iu[j]-1+j;i<j;ic++,i++){
  73.    /* set pointers to start dot product L row i, U col j
  74.    ** first assuming row i equal or longer than col j
  75.    ** then correcting if "reentrant" */
  76.     ilr=(iuc=iu1)-ic+il[i];
  77.     sum=0;
  78.     if(ic>iu1){
  79.       if((k=ilr-il[i-1])<=0){ 
  80.     iuc+=k=1-k; /* skip over implicit zeros in row i */
  81.     ilr+=k;}
  82.       for(;iuc<ic;iuc++,ilr++)
  83.     sum += au[iuc] * al[ilr];}
  84.     au[ic]=(au[ic]-sum)/al[ilr];};/* replace by factor U */
  85. /* now do L from 2nd column of row j */
  86.   for(i=(ic=(iu1=il[j-1]+1)+1)-il[j]+j;i<=j;ic++,i++){
  87.     if((k=(iuc=(ilr=iu1)-ic+iu[i]+1)-iu[i-1])<=0){
  88.       iuc+=k=1-k; /* skip reentrancy in column i */
  89.       ilr+=k;}
  90.     for(sum=0;ilr<ic;iuc++,ilr++)
  91.       sum += au[iuc] * al[ilr];
  92.     al[ic] -= sum;} /* replace with L */
  93. }
  94. }
  95. vbfwb(au,iu,al,il,n,b)
  96. int n,iu[],il[];
  97. OPTREAL au[],al[],b[];
  98. {
  99. register double sum;
  100. register int i,j,ic,ilr;
  101. /* b[n]: right side vector to be overwritten by
  102.   solution 
  103. solve L U b" = b
  104. writing b' and b" over b 
  105. U[j][j] =1 not stored */
  106. for(ilr=-1,j=0;j<n;j++){ /* b' = Linv b */
  107. /* multiply by L inverse */
  108.   for(sum=0,i=++ilr+j-il[j];i<j;ilr++,i++)
  109.     sum += al[ilr] * b[i];
  110.   b[j] = (b[j]-sum)/al[ilr];};/* L inv b */
  111. /* multiply by Lt inverse */
  112. for(ilr=iu[n-1];--j;)
  113.   for(i=j-1;ilr>iu[j-1];ilr--,i--)
  114.     b[i] -= b[j] * au[ilr];
  115. }
  116.