ceph/src/boost/libs/numeric/odeint/performance/c_lorenz.c

   1 #include <stdio.h>
   2 #include <time.h>
   3 #include <math.h>
   4
   5 void lorenz(const double *x, double *restrict y) {
   6     y[0] = 10.0 * (x[1] - x[0]);
   7     y[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
   8     y[2] = x[0] * x[1] - (8.0 / 3.0) * x[2];
   9 }
  10
  11 int main(int argc, const char *argv[])
  12 {
  13     const int nb_steps = 20000000;
  14     const double h = 1.0e-10;
  15     const double h2 = 0.5 * h;
  16     const double nb_loops = 21;
  17     double x[3];
  18     double y[3];
  19     double f1[3];
  20     double f2[3];
  21     double f3[3];
  22     double f4[3];
  23     double min_time = 1E6;
  24     clock_t begin, end;
  25     double time_spent;
  26
  27     for (int j = 0; j < nb_loops; j++) {
  28         x[0] = 8.5;
  29         x[1] = 3.1;
  30         x[2] = 1.2;
  31         begin = clock();
  32         for (int k = 0; k < nb_steps; k++) {
  33             lorenz(x, f1);
  34             for (int i = 0; i < 3; i++) {
  35                 y[i] = x[i] + h2 * f1[i];
  36             }
  37             lorenz(y, f2);
  38             for (int i = 0; i < 3; i++) {
  39                 y[i] = x[i] + h2 * f2[i];
  40             }
  41             lorenz(y, f3);
  42             for (int i = 0; i < 3; i++) {
  43                 y[i] = x[i] + h * f3[i];
  44             }
  45             lorenz(y, f4);
  46             for (int i = 0; i < 3; i++) {
  47                 x[i] = x[i] + h * (f1[i] + 2 * (f2[i] + f3[i]) + f4[i]) / 6.0;
  48             }
  49         }
  50         end = clock();
  51         min_time = fmin(min_time, (double)(end-begin)/CLOCKS_PER_SEC);
  52         printf("Result: %f\t runtime: %f\n", x[0], (double)(end-begin)/CLOCKS_PER_SEC);
  53     }
  54     printf("Minimal Runtime: %f\n", min_time);
  55
  56     return 0;
  57 }