1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022018 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 


22  using System;


23  using System.Diagnostics;


24  using System.Threading;


25  using HeuristicLab.Common;


26  using HeuristicLab.Core;


27  using HeuristicLab.Data;


28  using HeuristicLab.Optimization;


29  using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;


30 


31  namespace HeuristicLab.Algorithms.Benchmarks {


32  [Item("Linpack", "Linpack performance benchmark.")]


33  [StorableClass]


34  public sealed class Linpack : Benchmark {


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 


43  private CancellationToken cancellationToken;


44  private Stopwatch sw = new Stopwatch();


45 


46  [StorableConstructor]


47  private Linpack(bool deserializing) : base(deserializing) { }


48  private Linpack(Linpack original, Cloner cloner) : base(original, cloner) { }


49  public Linpack() { }


50 


51  public override IDeepCloneable Clone(Cloner cloner) {


52  return new Linpack(this, cloner);


53  }


54 


55  // implementation based on Java version: http://www.netlib.org/benchmark/linpackjava/


56  public override void Run(CancellationToken token, ResultCollection results) {


57  cancellationToken = token;


58  bool stopBenchmark = false;


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)));


147  results.Add(new Result("Total Mflops/s", new DoubleValue(mflops_result * Environment.ProcessorCount)));


148  resultAchieved = true;


149  }


150 


151  executionTime += sw.Elapsed;


152  if ((TimeLimit == null)  (TimeLimit.TotalMilliseconds == 0))


153  stopBenchmark = true;


154  else if (executionTime > TimeLimit)


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  }

