home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / rad386 / radiosit / src / gauss.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-21  |  1.6 KB  |  64 lines
  1. /*
  2.     Gauss Siedel Method of iterative Simultaneous Equation Solver.
  3.     -------------------------------------------------------------
  4.         Interface :
  5.                 gauss_iter(n,a,x,b)
  6.                 ----------
  7.                 where
  8.                         n - is the dimension
  9.                         a - Array of pointers
  10.                             pointing to array of float.
  11.                         b,x - arrays of float
  12.                         tolerance - A scalar float value
  13.                                 much smaller than 0.
  14.     
  15.     WARNING:
  16.         Matrix must be Diagonal Dominant.
  17.         Otherwise converging is not guaranteed.
  18. */
  19. #include <stdio.h>
  20. #include <math.h>
  21. extern int verbose_flag;
  22. double **matrix(),*vector();
  23. int gauss_iter(n,a,x,b,tolerance)
  24. int n;
  25. double **a,*x,b[];
  26. double tolerance;
  27. {
  28.     double newx;
  29.     double maxdelta;
  30.     int count=0;
  31.     if (verbose_flag)fprintf(stderr,"Iter Solve : ");
  32.     do{
  33.         int i;
  34.         maxdelta=0.0;
  35.         for(i=0; i < n; i++){
  36.             int j;
  37.             float delta;
  38.             newx = b[i];
  39.             for (j = 0; j < n; j++)
  40.                  if (j!=i)newx=newx-a[i][j]*x[j];
  41.             newx /= a[i][i];
  42.             delta = (float)fabs(newx - x[i]);
  43.             x[i]=newx;
  44.             if (delta > maxdelta) maxdelta=delta;
  45.         }
  46.         if (verbose_flag)fprintf(stderr,".");
  47.         count++;
  48.         if (count > n) break;
  49.     }while (maxdelta > tolerance);
  50.     if (verbose_flag)fprintf(stderr,"\n");
  51. }
  52. int diagonally_dominant(n,a)
  53. int n;
  54. double **a;
  55. {
  56.     int i,j;
  57.     double sum;
  58.     for(i=0;i<n;i++){
  59.         for(j = 0,sum = -fabs(a[i][i]);j<n;j++) sum += fabs(a[i][j]);
  60.         if(fabs(a[i][i])<sum)return(0);
  61.     }
  62.     return(1);
  63. }
  64.