1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022011 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.Drawing;


25  using System.Threading;


26  using HeuristicLab.Common;


27  using HeuristicLab.Core;


28  using HeuristicLab.Data;


29  using HeuristicLab.Optimization;


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


31 


32  namespace 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  }

