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

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

#1659:

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