Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Algorithms.Benchmarks/3.3/Linpack.cs @ 11123

Last change on this file since 11123 was 9456, checked in by swagner, 12 years ago

Updated copyright year and added some missing license headers (#1889)

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