Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Benchmarking/sources/HeuristicLab.Algorithms.Benchmarks/3.3/LinpackBenchmark.cs @ 6948

Last change on this file since 6948 was 6948, checked in by spimming, 11 years ago

#1659:

  • restructuring of the benchmarking algorithms
  • common interface for benchmarks
  • 'benchmark' class discovers benchmarking algorithms
  • added 'TimeLimit' and 'ChunkSize' parameters to benchmark
File size: 15.0 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2011 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Diagnostics;
24using System.Drawing;
25using System.Threading;
26using HeuristicLab.Common;
27using HeuristicLab.Core;
28using HeuristicLab.Data;
29using HeuristicLab.Optimization;
30using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
31
32namespace HeuristicLab.Algorithms.Benchmarks {
33  [Item("Linpack Algorithm", "Linpack benchmarking algorithm.")]
34  [StorableClass]
35  public class LinpackBenchmark : IBenchmark {
36    [Storable]
37    private byte[][] chunk;
38
39    private TimeSpan timeLimit;
40
41    private bool stopBenchmark;
42
43    private CancellationToken cancellationToken;
44
45    #region Benchmark Fields
46
47    private const int DEFAULT_PSIZE = 1500;
48
49    private double eps_result = 0.0;
50    private double mflops_result = 0.0;
51    private double residn_result = 0.0;
52    private double time_result = 0.0;
53    private double total = 0.0;
54
55    private Stopwatch sw = new Stopwatch();
56
57    #endregion
58
59    #region Properties
60
61    public byte[][] ChunkData {
62      get { return chunk; }
63      set { chunk = value; }
64    }
65
66    public TimeSpan TimeLimit {
67      get { return timeLimit; }
68      set { timeLimit = value; }
69    }
70
71    public string ItemName {
72      get { return ItemAttribute.GetName(this.GetType()); }
73    }
74
75    public string ItemDescription {
76      get { return ItemAttribute.GetDescription(this.GetType()); }
77    }
78
79    public Version ItemVersion {
80      get { return ItemAttribute.GetVersion(this.GetType()); }
81    }
82
83    public Image ItemImage {
84      get { return HeuristicLab.Common.Resources.VSImageLibrary.Event; }
85    }
86
87    #endregion
88
89    #region Costructors
90
91    public LinpackBenchmark() {
92    }
93
94    public LinpackBenchmark(LinpackBenchmark original, Cloner cloner) {
95      cloner.RegisterClonedObject(original, this);
96    }
97
98    #endregion
99
100    #region Linpack Benchmark
101    // implementation based on Java version: http://www.netlib.org/benchmark/linpackjava/
102
103    public void Run(CancellationToken token, ResultCollection results) {
104      cancellationToken = token;
105      stopBenchmark = false;
106      TimeSpan executionTime = new TimeSpan();
107      bool resultAchieved = false;
108      do {
109        int n = DEFAULT_PSIZE;
110        int ldaa = DEFAULT_PSIZE;
111        int lda = DEFAULT_PSIZE + 1;
112
113        double[][] a = new double[ldaa][];
114        double[] b = new double[ldaa];
115        double[] x = new double[ldaa];
116
117        double ops;
118        double norma;
119        double normx;
120        double resid;
121        int i;
122        int info;
123        int[] ipvt = new int[ldaa];
124
125        for (i = 0; i < ldaa; i++) {
126          a[i] = new double[lda];
127        }
128
129        ops = (2.0e0 * (((double)n) * n * n)) / 3.0 + 2.0 * (n * n);
130
131        norma = mathGen(a, lda, n, b);
132
133        if (cancellationToken.IsCancellationRequested) {
134          throw new OperationCanceledException(cancellationToken);
135        }
136
137        sw.Reset();
138        sw.Start();
139
140        info = dgefa(a, lda, n, ipvt);
141
142        if (cancellationToken.IsCancellationRequested) {
143          throw new OperationCanceledException(cancellationToken);
144        }
145
146        dgesl(a, lda, n, ipvt, b, 0);
147
148        sw.Stop();
149        total = sw.Elapsed.TotalMilliseconds / 1000;
150
151        if (cancellationToken.IsCancellationRequested) {
152          throw new OperationCanceledException(cancellationToken);
153        }
154
155        for (i = 0; i < n; i++) {
156          x[i] = b[i];
157        }
158
159        norma = mathGen(a, lda, n, b);
160
161        for (i = 0; i < n; i++) {
162          b[i] = -b[i];
163        }
164
165        dmxpy(n, b, n, lda, x, a);
166
167        resid = 0.0;
168        normx = 0.0;
169
170        for (i = 0; i < n; i++) {
171          resid = (resid > abs(b[i])) ? resid : abs(b[i]);
172          normx = (normx > abs(x[i])) ? normx : abs(x[i]);
173        }
174
175        eps_result = epslon((double)1.0);
176
177        residn_result = resid / (n * norma * normx * eps_result);
178        residn_result += 0.005; // for rounding
179        residn_result = (int)(residn_result * 100);
180        residn_result /= 100;
181
182        time_result = total;
183        time_result += 0.005; // for rounding
184        time_result = (int)(time_result * 100);
185        time_result /= 100;
186
187        mflops_result = ops / (1.0e6 * total);
188        mflops_result += 0.0005; // for rounding
189        mflops_result = (int)(mflops_result * 1000);
190        mflops_result /= 1000;
191
192        if (!resultAchieved) {
193          results.Add(new Result("Mflops/s", new DoubleValue(mflops_result)));
194          results.Add(new Result("total Mflops/s", new DoubleValue(mflops_result * Environment.ProcessorCount)));
195          resultAchieved = true;
196        }
197
198        executionTime += sw.Elapsed;
199        if (timeLimit == null)
200          stopBenchmark = true;
201        else if (executionTime > timeLimit)
202          stopBenchmark = true;
203      } while (!stopBenchmark);
204    }
205
206    private double abs(double d) {
207      return (d >= 0) ? d : -d;
208    }
209
210    private double mathGen(double[][] a, int lda, int n, double[] b) {
211      Random gen;
212      double norma;
213      int init, i, j;
214
215      init = 1325;
216      norma = 0.0;
217
218      gen = new Random(init);
219
220      if (cancellationToken.IsCancellationRequested) {
221        throw new OperationCanceledException(cancellationToken);
222      }
223
224      // Next two for() statements switched.  Solver wants
225      // matrix in column order. --dmd 3/3/97
226
227      for (i = 0; i < n; i++) {
228        for (j = 0; j < n; j++) {
229          a[j][i] = gen.NextDouble() - .5;
230          norma = (a[j][i] > norma) ? a[j][i] : norma;
231        }
232      }
233
234      for (i = 0; i < n; i++) {
235        b[i] = 0.0;
236      }
237
238      for (j = 0; j < n; j++) {
239        for (i = 0; i < n; i++) {
240          b[i] += a[j][i];
241        }
242      }
243
244      return norma;
245    }
246
247    private int dgefa(double[][] a, int lda, int n, int[] ipvt) {
248      double[] col_k, col_j;
249      double t;
250      int j, k, kp1, l, nm1;
251      int info;
252
253      if (cancellationToken.IsCancellationRequested) {
254        throw new OperationCanceledException(cancellationToken);
255      }
256
257      // gaussian elimination with partial pivoting
258
259      info = 0;
260      nm1 = n - 1;
261      if (nm1 >= 0) {
262        for (k = 0; k < nm1; k++) {
263          col_k = a[k];
264          kp1 = k + 1;
265
266          // find l = pivot index
267
268          l = idamax(n - k, col_k, k, 1) + k;
269          ipvt[k] = l;
270
271          // zero pivot implies this column already triangularized
272
273          if (col_k[l] != 0) {
274            // interchange if necessary
275
276            if (l != k) {
277              t = col_k[l];
278              col_k[l] = col_k[k];
279              col_k[k] = t;
280            }
281
282            if (cancellationToken.IsCancellationRequested) {
283              throw new OperationCanceledException(cancellationToken);
284            }
285
286            // compute multipliers
287
288            t = -1.0 / col_k[k];
289            dscal(n - (kp1), t, col_k, kp1, 1);
290
291            if (cancellationToken.IsCancellationRequested) {
292              throw new OperationCanceledException(cancellationToken);
293            }
294
295            // row elimination with column indexing
296
297            for (j = kp1; j < n; j++) {
298              col_j = a[j];
299              t = col_j[l];
300              if (l != k) {
301                col_j[l] = col_j[k];
302                col_j[k] = t;
303              }
304              daxpy(n - (kp1), t, col_k, kp1, 1,
305                col_j, kp1, 1);
306            }
307          } else {
308            info = k;
309          }
310        }
311      }
312
313      ipvt[n - 1] = n - 1;
314      if (a[(n - 1)][(n - 1)] == 0) info = n - 1;
315
316      return info;
317    }
318
319    private void dgesl(double[][] a, int lda, int n, int[] ipvt, double[] b, int job) {
320      double t;
321      int k, kb, l, nm1, kp1;
322
323      if (cancellationToken.IsCancellationRequested) {
324        throw new OperationCanceledException(cancellationToken);
325      }
326
327      nm1 = n - 1;
328      if (job == 0) {
329        // job = 0 , solve  a * x = b.  first solve  l*y = b
330
331        if (nm1 >= 1) {
332          for (k = 0; k < nm1; k++) {
333            l = ipvt[k];
334            t = b[l];
335            if (l != k) {
336              b[l] = b[k];
337              b[k] = t;
338            }
339            kp1 = k + 1;
340            daxpy(n - (kp1), t, a[k], kp1, 1, b, kp1, 1);
341          }
342        }
343
344        if (cancellationToken.IsCancellationRequested) {
345          throw new OperationCanceledException(cancellationToken);
346        }
347
348        // now solve  u*x = y
349
350        for (kb = 0; kb < n; kb++) {
351          k = n - (kb + 1);
352          b[k] /= a[k][k];
353          t = -b[k];
354          daxpy(k, t, a[k], 0, 1, b, 0, 1);
355        }
356      } else {
357        // job = nonzero, solve  trans(a) * x = b.  first solve  trans(u)*y = b
358
359        for (k = 0; k < n; k++) {
360          t = ddot(k, a[k], 0, 1, b, 0, 1);
361          b[k] = (b[k] - t) / a[k][k];
362        }
363
364        if (cancellationToken.IsCancellationRequested) {
365          throw new OperationCanceledException(cancellationToken);
366        }
367
368        // now solve trans(l)*x = y
369
370        if (nm1 >= 1) {
371          //for (kb = 1; kb < nm1; kb++) {
372          for (kb = 0; kb < nm1; kb++) {
373            k = n - (kb + 1);
374            kp1 = k + 1;
375            b[k] += ddot(n - (kp1), a[k], kp1, 1, b, kp1, 1);
376            l = ipvt[k];
377            if (l != k) {
378              t = b[l];
379              b[l] = b[k];
380              b[k] = t;
381            }
382          }
383        }
384      }
385    }
386
387    private void daxpy(int n, double da, double[] dx, int dx_off, int incx, double[] dy, int dy_off, int incy) {
388      int i, ix, iy;
389
390      if (cancellationToken.IsCancellationRequested) {
391        throw new OperationCanceledException(cancellationToken);
392      }
393
394      if ((n > 0) && (da != 0)) {
395        if (incx != 1 || incy != 1) {
396
397          // code for unequal increments or equal increments not equal to 1
398
399          ix = 0;
400          iy = 0;
401          if (incx < 0) ix = (-n + 1) * incx;
402          if (incy < 0) iy = (-n + 1) * incy;
403          for (i = 0; i < n; i++) {
404            dy[iy + dy_off] += da * dx[ix + dx_off];
405            ix += incx;
406            iy += incy;
407          }
408          return;
409        } else {
410          // code for both increments equal to 1
411
412          for (i = 0; i < n; i++)
413            dy[i + dy_off] += da * dx[i + dx_off];
414        }
415      }
416    }
417
418    private double ddot(int n, double[] dx, int dx_off, int incx, double[] dy, int dy_off, int incy) {
419      double dtemp = 0;
420      int i, ix, iy;
421
422      if (cancellationToken.IsCancellationRequested) {
423        throw new OperationCanceledException(cancellationToken);
424      }
425
426      if (n > 0) {
427        if (incx != 1 || incy != 1) {
428          // code for unequal increments or equal increments not equal to 1
429
430          ix = 0;
431          iy = 0;
432          if (incx < 0) ix = (-n + 1) * incx;
433          if (incy < 0) iy = (-n + 1) * incy;
434          for (i = 0; i < n; i++) {
435            dtemp += dx[ix + dx_off] * dy[iy + dy_off];
436            ix += incx;
437            iy += incy;
438          }
439        } else {
440          // code for both increments equal to 1
441
442          for (i = 0; i < n; i++)
443            dtemp += dx[i + dx_off] * dy[i + dy_off];
444        }
445      }
446      return (dtemp);
447    }
448
449    private void dscal(int n, double da, double[] dx, int dx_off, int incx) {
450      int i, nincx;
451
452      if (cancellationToken.IsCancellationRequested) {
453        throw new OperationCanceledException(cancellationToken);
454      }
455
456      if (n > 0) {
457        if (incx != 1) {
458          // code for increment not equal to 1
459
460          nincx = n * incx;
461          for (i = 0; i < nincx; i += incx)
462            dx[i + dx_off] *= da;
463        } else {
464          // code for increment equal to 1
465
466          for (i = 0; i < n; i++)
467            dx[i + dx_off] *= da;
468        }
469      }
470    }
471
472    private int idamax(int n, double[] dx, int dx_off, int incx) {
473      double dmax, dtemp;
474      int i, ix, itemp = 0;
475
476      if (cancellationToken.IsCancellationRequested) {
477        throw new OperationCanceledException(cancellationToken);
478      }
479
480      if (n < 1) {
481        itemp = -1;
482      } else if (n == 1) {
483        itemp = 0;
484      } else if (incx != 1) {
485        // code for increment not equal to 1
486
487        dmax = (dx[dx_off] < 0.0) ? -dx[dx_off] : dx[dx_off];
488        ix = 1 + incx;
489        for (i = 0; i < n; i++) {
490          dtemp = (dx[ix + dx_off] < 0.0) ? -dx[ix + dx_off] : dx[ix + dx_off];
491          if (dtemp > dmax) {
492            itemp = i;
493            dmax = dtemp;
494          }
495          ix += incx;
496        }
497      } else {
498        // code for increment equal to 1
499
500        itemp = 0;
501        dmax = (dx[dx_off] < 0.0) ? -dx[dx_off] : dx[dx_off];
502        for (i = 0; i < n; i++) {
503          dtemp = (dx[i + dx_off] < 0.0) ? -dx[i + dx_off] : dx[i + dx_off];
504          if (dtemp > dmax) {
505            itemp = i;
506            dmax = dtemp;
507          }
508        }
509      }
510      return (itemp);
511    }
512
513    private double epslon(double x) {
514      double a, b, c, eps;
515
516      a = 4.0e0 / 3.0e0;
517      eps = 0;
518      while (eps == 0) {
519        b = a - 1.0;
520        c = b + b + b;
521        eps = abs(c - 1.0);
522      }
523      return (eps * abs(x));
524    }
525
526    private void dmxpy(int n1, double[] y, int n2, int ldm, double[] x, double[][] m) {
527      int j, i;
528
529      // cleanup odd vector
530      for (j = 0; j < n2; j++) {
531        for (i = 0; i < n1; i++) {
532          y[i] += x[j] * m[j][i];
533        }
534      }
535    }
536
537    #endregion
538
539    #region Clone
540
541    public IDeepCloneable Clone(Cloner cloner) {
542      return new LinpackBenchmark(this, cloner);
543    }
544
545    public object Clone() {
546      return Clone(new Cloner());
547    }
548
549    #endregion
550
551    #region Events
552
553    public event EventHandler ItemImageChanged;
554    protected virtual void OnItemImageChanged() {
555      EventHandler handler = ItemImageChanged;
556      if (handler != null) handler(this, EventArgs.Empty);
557    }
558
559    public event EventHandler ToStringChanged;
560    protected virtual void OnToStringChanged() {
561      EventHandler handler = ToStringChanged;
562      if (handler != null) handler(this, EventArgs.Empty);
563    }
564
565    #endregion
566  }
567}
Note: See TracBrowser for help on using the repository browser.