bmks/polykernels1/trisolv.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <sys/time.h>
   5
   6 #ifdef SELFTEST
   7         #define NMAX 150
   8 #else
   9         #define NMAX 3000
  10         #define MEASURE_TIME 1
  11 #endif // SELFTEST
  12
  13 // #define S1(i,j,k) B[j][i]=B[j][i]-L[j][k]*B[k][i]
  14 // #define S2(i,j) B[j][i]=B[j][i]/L[j][j]
  15
  16 static double B[NMAX][NMAX], L[NMAX][NMAX];
  17 static double BB[NMAX][NMAX];
  18
  19 void trisolv(long N) {
  20
  21     long i,j,k;
  22
  23     for (i=0;i<=N-1;i++) {
  24
  25         for (j=0;j<=N-1;j++) {
  26
  27             for (k=0;k<=j-1;k++) {
  28
  29                 B[j][i]=B[j][i]-L[j][k]*B[k][i];  //S1 ;
  30
  31             }
  32
  33             B[j][i]=B[j][i]/L[j][j]; // S2 ;
  34
  35         } // for j
  36
  37     } // for i
  38
  39 } // end trisolv()
  40
  41 #ifdef SELFTEST
  42 void ref(long N) {
  43
  44     long i,j,k;
  45
  46     for (i=0;i<=N-1;i++) {
  47
  48         for (j=0;j<=N-1;j++) {
  49
  50             for (k=0;k<=j-1;k++) {
  51
  52                 BB[j][i]=BB[j][i]-L[j][k]*BB[k][i];  //S1 ;
  53
  54             }
  55
  56             BB[j][i]=BB[j][i]/L[j][j]; // S2 ;
  57
  58         } // for j
  59
  60     } // for i
  61 }
  62 #endif // SELFTEST
  63
  64 int main()
  65 {
  66         struct timeval start;
  67         struct timeval end;
  68     long N=NMAX;
  69         int i,j;
  70         int errors = 0;
  71
  72         for (i = 0; i < NMAX; i++) {
  73                 for (j = 0; j < NMAX; j++) {
  74                         B[i][j]  = BB[i][j] = i+j;
  75                     L[i][j] = (i+j+3.45) *i*j*0.5;
  76                 }
  77         }
  78
  79         gettimeofday(&start, NULL);
  80     trisolv(N);
  81         gettimeofday(&end, NULL);
  82
  83 #ifdef SELFTEST
  84         ref(N);
  85
  86         for (i = 0; i < NMAX; i++) {
  87                 for (j = 0; j < NMAX ; j++) {
  88                         if (fabs(B[i][j]/BB[i][j] -1) > 10e-8) {
  89                                 errors++;
  90 #ifdef ST_VERB
  91                                 printf("Error: expected BB[%d][%d] = %lf found B[%d][%d] = %lf\t=>\tabs(a/b -1) = %lf\n",
  92                                                         i, j, BB[i][j], i, j, B[i][j], fabs(BB[i][j]/B[i][j] -1) );
  93 #endif // ST_VERB
  94                         }
  95                 }
  96         }
  97
  98         printf("SELFTEST: TRISOLV ");
  99         if(errors != 0) printf("FAILED with %d errors!\n", errors);
 100         else printf("SUCCESS!\n");
 101 #endif
 102
 103 #ifdef MEASURE_TIME
 104         printf("Elapsed Time=%lf\n",
 105                         (end.tv_sec - start.tv_sec +
 106                         (double)(end.tv_usec - start.tv_usec)/1000000));
 107 #endif
 108     return (errors != 0);
 109 }