Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Benchmarking/sources/HeuristicLab.Algorithms.Benchmarks/3.3/Linpack.cs @ 6920

Last change on this file since 6920 was 6920, checked in by ascheibe, 12 years ago

#1659 added benchmarking algorithms

File size: 12.5 KB
Line 
1using System;
2using System.Diagnostics;
3using System.Threading;
4using System.Threading.Tasks;
5using HeuristicLab.Common;
6using HeuristicLab.Core;
7using HeuristicLab.Data;
8using HeuristicLab.Optimization;
9using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
10
11namespace HeuristicLab.Algorithms.Benchmarks {
12  [Item("Linpack Algorithm", "A Linpack benchmark algorithm.")]
13  [Creatable("Benchmarks")]
14  [StorableClass]
15  public class Linpack : Algorithm {
16    private DateTime lastUpdateTime;
17
18    [Storable]
19    private ResultCollection results;
20
21    #region Benchmark Fields
22
23    private const int DEFAULT_PSIZE = 1500;
24
25    private double eps_result = 0.0;
26    private double mflops_result = 0.0;
27    private double residn_result = 0.0;
28    private double time_result = 0.0;
29    private double total = 0.0;
30
31    private Stopwatch sw = new Stopwatch();
32
33    #endregion
34
35    #region Properties
36
37    public override ResultCollection Results {
38      get { return results; }
39    }
40
41    #endregion
42
43    #region Costructors
44
45    public Linpack()
46      : base() {
47      results = new ResultCollection();
48    }
49
50    private Linpack(Linpack original, Cloner cloner)
51      : base(original, cloner) {
52      results = new ResultCollection();
53    }
54
55    #endregion
56
57    public override IDeepCloneable Clone(Cloner cloner) {
58      return new Linpack(this, cloner);
59    }
60
61    public override void Prepare() {
62      results.Clear();
63      OnPrepared();
64    }
65
66    public override void Start() {
67      var cancellationTokenSource = new CancellationTokenSource();
68      OnStarted();
69      Task task = Task.Factory.StartNew(Run, cancellationTokenSource.Token, cancellationTokenSource.Token);
70      task.ContinueWith(t => {
71        try {
72          t.Wait();
73        }
74        catch (AggregateException ex) {
75          try {
76            ex.Flatten().Handle(x => x is OperationCanceledException);
77          }
78          catch (AggregateException remaining) {
79            if (remaining.InnerExceptions.Count == 1) OnExceptionOccurred(remaining.InnerExceptions[0]);
80            else OnExceptionOccurred(remaining);
81          }
82        }
83        cancellationTokenSource.Dispose();
84        cancellationTokenSource = null;
85        OnStopped();
86      });
87    }
88
89    private void Run(object state) {
90      CancellationToken cancellationToken = (CancellationToken)state;
91      lastUpdateTime = DateTime.Now;
92      System.Timers.Timer timer = new System.Timers.Timer(250);
93      timer.AutoReset = true;
94      timer.Elapsed += new System.Timers.ElapsedEventHandler(timer_Elapsed);
95      timer.Start();
96      try {
97        RunBenchmark();
98      }
99      finally {
100        timer.Elapsed -= new System.Timers.ElapsedEventHandler(timer_Elapsed);
101        timer.Stop();
102        ExecutionTime += DateTime.Now - lastUpdateTime;
103      }
104
105      cancellationToken.ThrowIfCancellationRequested();
106    }
107
108    private void timer_Elapsed(object sender, System.Timers.ElapsedEventArgs e) {
109      System.Timers.Timer timer = (System.Timers.Timer)sender;
110      timer.Enabled = false;
111      DateTime now = DateTime.Now;
112      ExecutionTime += now - lastUpdateTime;
113      lastUpdateTime = now;
114      timer.Enabled = true;
115    }
116
117    #region Linpack Benchmark
118
119    private void RunBenchmark() {
120      int n = DEFAULT_PSIZE;
121      int ldaa = DEFAULT_PSIZE;
122      int lda = DEFAULT_PSIZE + 1;
123
124      double[][] a = new double[ldaa][];
125      double[] b = new double[ldaa];
126      double[] x = new double[ldaa];
127
128      double ops;
129      double norma;
130      double normx;
131      double resid;
132      int i;
133      int info;
134      int[] ipvt = new int[ldaa];
135
136      for (i = 0; i < ldaa; i++) {
137        a[i] = new double[lda];
138      }
139
140      ops = (2.0e0 * (((double)n) * n * n)) / 3.0 + 2.0 * (n * n);
141
142      norma = mathGen(a, lda, n, b);
143
144      sw.Reset();
145      sw.Start();
146
147      info = dgefa(a, lda, n, ipvt);
148
149      dgesl(a, lda, n, ipvt, b, 0);
150
151      sw.Stop();
152      total = sw.Elapsed.TotalMilliseconds / 1000;
153
154      for (i = 0; i < n; i++) {
155        x[i] = b[i];
156      }
157
158      norma = mathGen(a, lda, n, b);
159
160      for (i = 0; i < n; i++) {
161        b[i] = -b[i];
162      }
163
164      dmxpy(n, b, n, lda, x, a);
165
166      resid = 0.0;
167      normx = 0.0;
168
169      for (i = 0; i < n; i++) {
170        resid = (resid > abs(b[i])) ? resid : abs(b[i]);
171        normx = (normx > abs(x[i])) ? normx : abs(x[i]);
172      }
173
174      eps_result = epslon((double)1.0);
175
176      residn_result = resid / (n * norma * normx * eps_result);
177      residn_result += 0.005; // for rounding
178      residn_result = (int)(residn_result * 100);
179      residn_result /= 100;
180
181      time_result = total;
182      time_result += 0.005; // for rounding
183      time_result = (int)(time_result * 100);
184      time_result /= 100;
185
186      mflops_result = ops / (1.0e6 * total);
187      mflops_result += 0.0005; // for rounding
188      mflops_result = (int)(mflops_result * 1000);
189      mflops_result /= 1000;
190
191      //System.Console.WriteLine("Mflops/s: " + mflops_result + "  Time: " + time_result + " secs" + "  Norm Res: " + residn_result + "  Precision: " + eps_result);
192
193      Results.Add(new Result("Mflops/s", new DoubleValue(mflops_result)));
194      //Results.Add(new Result("ca. Mflops/s", new DoubleValue(mflops_result * Environment.ProcessorCount)));
195    }
196
197    private double abs(double d) {
198      return (d >= 0) ? d : -d;
199    }
200
201    private double mathGen(double[][] a, int lda, int n, double[] b) {
202      Random gen;
203      double norma;
204      int init, i, j;
205
206      init = 1325;
207      norma = 0.0;
208
209      gen = new Random(init);
210
211      // Next two for() statements switched.  Solver wants
212      // matrix in column order. --dmd 3/3/97
213
214      for (i = 0; i < n; i++) {
215        for (j = 0; j < n; j++) {
216          a[j][i] = gen.NextDouble() - .5;
217          norma = (a[j][i] > norma) ? a[j][i] : norma;
218        }
219      }
220
221      for (i = 0; i < n; i++) {
222        b[i] = 0.0;
223      }
224
225      for (j = 0; j < n; j++) {
226        for (i = 0; i < n; i++) {
227          b[i] += a[j][i];
228        }
229      }
230
231      return norma;
232    }
233
234    private int dgefa(double[][] a, int lda, int n, int[] ipvt) {
235      double[] col_k, col_j;
236      double t;
237      int j, k, kp1, l, nm1;
238      int info;
239
240      // gaussian elimination with partial pivoting
241
242      info = 0;
243      nm1 = n - 1;
244      if (nm1 >= 0) {
245        for (k = 0; k < nm1; k++) {
246          col_k = a[k];
247          kp1 = k + 1;
248
249          // find l = pivot index
250
251          l = idamax(n - k, col_k, k, 1) + k;
252          ipvt[k] = l;
253
254          // zero pivot implies this column already triangularized
255
256          if (col_k[l] != 0) {
257            // interchange if necessary
258
259            if (l != k) {
260              t = col_k[l];
261              col_k[l] = col_k[k];
262              col_k[k] = t;
263            }
264
265            // compute multipliers
266
267            t = -1.0 / col_k[k];
268            dscal(n - (kp1), t, col_k, kp1, 1);
269
270            // row elimination with column indexing
271
272            for (j = kp1; j < n; j++) {
273              col_j = a[j];
274              t = col_j[l];
275              if (l != k) {
276                col_j[l] = col_j[k];
277                col_j[k] = t;
278              }
279              daxpy(n - (kp1), t, col_k, kp1, 1,
280                col_j, kp1, 1);
281            }
282          } else {
283            info = k;
284          }
285        }
286      }
287
288      ipvt[n - 1] = n - 1;
289      if (a[(n - 1)][(n - 1)] == 0) info = n - 1;
290
291      return info;
292    }
293
294    private void dgesl(double[][] a, int lda, int n, int[] ipvt, double[] b, int job) {
295      double t;
296      int k, kb, l, nm1, kp1;
297
298      nm1 = n - 1;
299      if (job == 0) {
300        // job = 0 , solve  a * x = b.  first solve  l*y = b
301
302        if (nm1 >= 1) {
303          for (k = 0; k < nm1; k++) {
304            l = ipvt[k];
305            t = b[l];
306            if (l != k) {
307              b[l] = b[k];
308              b[k] = t;
309            }
310            kp1 = k + 1;
311            daxpy(n - (kp1), t, a[k], kp1, 1, b, kp1, 1);
312          }
313        }
314
315        // now solve  u*x = y
316
317        for (kb = 0; kb < n; kb++) {
318          k = n - (kb + 1);
319          b[k] /= a[k][k];
320          t = -b[k];
321          daxpy(k, t, a[k], 0, 1, b, 0, 1);
322        }
323      } else {
324        // job = nonzero, solve  trans(a) * x = b.  first solve  trans(u)*y = b
325
326        for (k = 0; k < n; k++) {
327          t = ddot(k, a[k], 0, 1, b, 0, 1);
328          b[k] = (b[k] - t) / a[k][k];
329        }
330
331        // now solve trans(l)*x = y
332
333        if (nm1 >= 1) {
334          //for (kb = 1; kb < nm1; kb++) {
335          for (kb = 0; kb < nm1; kb++) {
336            k = n - (kb + 1);
337            kp1 = k + 1;
338            b[k] += ddot(n - (kp1), a[k], kp1, 1, b, kp1, 1);
339            l = ipvt[k];
340            if (l != k) {
341              t = b[l];
342              b[l] = b[k];
343              b[k] = t;
344            }
345          }
346        }
347      }
348    }
349
350    private void daxpy(int n, double da, double[] dx, int dx_off, int incx, double[] dy, int dy_off, int incy) {
351      int i, ix, iy;
352
353      if ((n > 0) && (da != 0)) {
354        if (incx != 1 || incy != 1) {
355
356          // code for unequal increments or equal increments not equal to 1
357
358          ix = 0;
359          iy = 0;
360          if (incx < 0) ix = (-n + 1) * incx;
361          if (incy < 0) iy = (-n + 1) * incy;
362          for (i = 0; i < n; i++) {
363            dy[iy + dy_off] += da * dx[ix + dx_off];
364            ix += incx;
365            iy += incy;
366          }
367          return;
368        } else {
369          // code for both increments equal to 1
370
371          for (i = 0; i < n; i++)
372            dy[i + dy_off] += da * dx[i + dx_off];
373        }
374      }
375    }
376
377    private double ddot(int n, double[] dx, int dx_off, int incx, double[] dy, int dy_off, int incy) {
378      double dtemp = 0;
379      int i, ix, iy;
380
381      if (n > 0) {
382        if (incx != 1 || incy != 1) {
383          // code for unequal increments or equal increments not equal to 1
384
385          ix = 0;
386          iy = 0;
387          if (incx < 0) ix = (-n + 1) * incx;
388          if (incy < 0) iy = (-n + 1) * incy;
389          for (i = 0; i < n; i++) {
390            dtemp += dx[ix + dx_off] * dy[iy + dy_off];
391            ix += incx;
392            iy += incy;
393          }
394        } else {
395          // code for both increments equal to 1
396
397          for (i = 0; i < n; i++)
398            dtemp += dx[i + dx_off] * dy[i + dy_off];
399        }
400      }
401      return (dtemp);
402    }
403
404    private void dscal(int n, double da, double[] dx, int dx_off, int incx) {
405      int i, nincx;
406
407      if (n > 0) {
408        if (incx != 1) {
409          // code for increment not equal to 1
410
411          nincx = n * incx;
412          for (i = 0; i < nincx; i += incx)
413            dx[i + dx_off] *= da;
414        } else {
415          // code for increment equal to 1
416
417          for (i = 0; i < n; i++)
418            dx[i + dx_off] *= da;
419        }
420      }
421    }
422
423    private int idamax(int n, double[] dx, int dx_off, int incx) {
424      double dmax, dtemp;
425      int i, ix, itemp = 0;
426
427      if (n < 1) {
428        itemp = -1;
429      } else if (n == 1) {
430        itemp = 0;
431      } else if (incx != 1) {
432        // code for increment not equal to 1
433
434        dmax = (dx[dx_off] < 0.0) ? -dx[dx_off] : dx[dx_off];
435        ix = 1 + incx;
436        for (i = 0; i < n; i++) {
437          dtemp = (dx[ix + dx_off] < 0.0) ? -dx[ix + dx_off] : dx[ix + dx_off];
438          if (dtemp > dmax) {
439            itemp = i;
440            dmax = dtemp;
441          }
442          ix += incx;
443        }
444      } else {
445        // code for increment equal to 1
446
447        itemp = 0;
448        dmax = (dx[dx_off] < 0.0) ? -dx[dx_off] : dx[dx_off];
449        for (i = 0; i < n; i++) {
450          dtemp = (dx[i + dx_off] < 0.0) ? -dx[i + dx_off] : dx[i + dx_off];
451          if (dtemp > dmax) {
452            itemp = i;
453            dmax = dtemp;
454          }
455        }
456      }
457      return (itemp);
458    }
459
460    private double epslon(double x) {
461      double a, b, c, eps;
462
463      a = 4.0e0 / 3.0e0;
464      eps = 0;
465      while (eps == 0) {
466        b = a - 1.0;
467        c = b + b + b;
468        eps = abs(c - 1.0);
469      }
470      return (eps * abs(x));
471    }
472
473    private void dmxpy(int n1, double[] y, int n2, int ldm, double[] x, double[][] m) {
474      int j, i;
475
476      // cleanup odd vector
477      for (j = 0; j < n2; j++) {
478        for (i = 0; i < n1; i++) {
479          y[i] += x[j] * m[j][i];
480        }
481      }
482    }
483
484    #endregion
485  }
486}
Note: See TracBrowser for help on using the repository browser.