source: trunk/sources/HeuristicLab.Algorithms.Benchmarks/3.3/LinpackBenchmark.cs @ 7007

Last change on this file since 7007 was 7007, checked in by ascheibe, 9 years ago

#1659 fixed small typo

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