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 


39  private TimeSpan timeLimit;


40 


41  private bool stopBenchmark;


42 


43  private CancellationToken cancellationToken;


44 


45  #region Benchmark Fields


46 


47  private const int DEFAULT_PSIZE = 1500;


48 


49  private double eps_result = 0.0;


50  private double mflops_result = 0.0;


51  private double residn_result = 0.0;


52  private double time_result = 0.0;


53  private double total = 0.0;


54 


55  private Stopwatch sw = new Stopwatch();


56 


57  #endregion


58 


59  #region Properties


60 


61  public byte[][] ChunkData {


62  get { return chunk; }


63  set { chunk = value; }


64  }


65 


66  public TimeSpan TimeLimit {


67  get { return timeLimit; }


68  set { timeLimit = value; }


69  }


70 


71  public string ItemName {


72  get { return ItemAttribute.GetName(this.GetType()); }


73  }


74 


75  public string ItemDescription {


76  get { return ItemAttribute.GetDescription(this.GetType()); }


77  }


78 


79  public Version ItemVersion {


80  get { return ItemAttribute.GetVersion(this.GetType()); }


81  }


82 


83  public Image ItemImage {


84  get { return HeuristicLab.Common.Resources.VSImageLibrary.Event; }


85  }


86 


87  #endregion


88 


89  #region Costructors


90 


91  public LinpackBenchmark() {


92  }


93 


94  public LinpackBenchmark(LinpackBenchmark original, Cloner cloner) {


95  cloner.RegisterClonedObject(original, this);


96  }


97 


98  #endregion


99 


100  #region Linpack Benchmark


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


102 


103  public void Run(CancellationToken token, ResultCollection results) {


104  cancellationToken = token;


105  stopBenchmark = false;


106  TimeSpan executionTime = new TimeSpan();


107  bool resultAchieved = false;


108  do {


109  int n = DEFAULT_PSIZE;


110  int ldaa = DEFAULT_PSIZE;


111  int lda = DEFAULT_PSIZE + 1;


112 


113  double[][] a = new double[ldaa][];


114  double[] b = new double[ldaa];


115  double[] x = new double[ldaa];


116 


117  double ops;


118  double norma;


119  double normx;


120  double resid;


121  int i;


122  int info;


123  int[] ipvt = new int[ldaa];


124 


125  for (i = 0; i < ldaa; i++) {


126  a[i] = new double[lda];


127  }


128 


129  ops = (2.0e0 * (((double)n) * n * n)) / 3.0 + 2.0 * (n * n);


130 


131  norma = mathGen(a, lda, n, b);


132 


133  if (cancellationToken.IsCancellationRequested) {


134  throw new OperationCanceledException(cancellationToken);


135  }


136 


137  sw.Reset();


138  sw.Start();


139 


140  info = dgefa(a, lda, n, ipvt);


141 


142  if (cancellationToken.IsCancellationRequested) {


143  throw new OperationCanceledException(cancellationToken);


144  }


145 


146  dgesl(a, lda, n, ipvt, b, 0);


147 


148  sw.Stop();


149  total = sw.Elapsed.TotalMilliseconds / 1000;


150 


151  if (cancellationToken.IsCancellationRequested) {


152  throw new OperationCanceledException(cancellationToken);


153  }


154 


155  for (i = 0; i < n; i++) {


156  x[i] = b[i];


157  }


158 


159  norma = mathGen(a, lda, n, b);


160 


161  for (i = 0; i < n; i++) {


162  b[i] = b[i];


163  }


164 


165  dmxpy(n, b, n, lda, x, a);


166 


167  resid = 0.0;


168  normx = 0.0;


169 


170  for (i = 0; i < n; i++) {


171  resid = (resid > abs(b[i])) ? resid : abs(b[i]);


172  normx = (normx > abs(x[i])) ? normx : abs(x[i]);


173  }


174 


175  eps_result = epslon((double)1.0);


176 


177  residn_result = resid / (n * norma * normx * eps_result);


178  residn_result += 0.005; // for rounding


179  residn_result = (int)(residn_result * 100);


180  residn_result /= 100;


181 


182  time_result = total;


183  time_result += 0.005; // for rounding


184  time_result = (int)(time_result * 100);


185  time_result /= 100;


186 


187  mflops_result = ops / (1.0e6 * total);


188  mflops_result += 0.0005; // for rounding


189  mflops_result = (int)(mflops_result * 1000);


190  mflops_result /= 1000;


191 


192  if (!resultAchieved) {


193  results.Add(new Result("Mflops/s", new DoubleValue(mflops_result)));


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


195  resultAchieved = true;


196  }


197 


198  executionTime += sw.Elapsed;


199  if (timeLimit == null)


200  stopBenchmark = true;


201  else if (executionTime > timeLimit)


202  stopBenchmark = true;


203  } while (!stopBenchmark);


204  }


205 


206  private double abs(double d) {


207  return (d >= 0) ? d : d;


208  }


209 


210  private double mathGen(double[][] a, int lda, int n, double[] b) {


211  Random gen;


212  double norma;


213  int init, i, j;


214 


215  init = 1325;


216  norma = 0.0;


217 


218  gen = new Random(init);


219 


220  if (cancellationToken.IsCancellationRequested) {


221  throw new OperationCanceledException(cancellationToken);


222  }


223 


224  // Next two for() statements switched. Solver wants


225  // matrix in column order. dmd 3/3/97


226 


227  for (i = 0; i < n; i++) {


228  for (j = 0; j < n; j++) {


229  a[j][i] = gen.NextDouble()  .5;


230  norma = (a[j][i] > norma) ? a[j][i] : norma;


231  }


232  }


233 


234  for (i = 0; i < n; i++) {


235  b[i] = 0.0;


236  }


237 


238  for (j = 0; j < n; j++) {


239  for (i = 0; i < n; i++) {


240  b[i] += a[j][i];


241  }


242  }


243 


244  return norma;


245  }


246 


247  private int dgefa(double[][] a, int lda, int n, int[] ipvt) {


248  double[] col_k, col_j;


249  double t;


250  int j, k, kp1, l, nm1;


251  int info;


252 


253  if (cancellationToken.IsCancellationRequested) {


254  throw new OperationCanceledException(cancellationToken);


255  }


256 


257  // gaussian elimination with partial pivoting


258 


259  info = 0;


260  nm1 = n  1;


261  if (nm1 >= 0) {


262  for (k = 0; k < nm1; k++) {


263  col_k = a[k];


264  kp1 = k + 1;


265 


266  // find l = pivot index


267 


268  l = idamax(n  k, col_k, k, 1) + k;


269  ipvt[k] = l;


270 


271  // zero pivot implies this column already triangularized


272 


273  if (col_k[l] != 0) {


274  // interchange if necessary


275 


276  if (l != k) {


277  t = col_k[l];


278  col_k[l] = col_k[k];


279  col_k[k] = t;


280  }


281 


282  if (cancellationToken.IsCancellationRequested) {


283  throw new OperationCanceledException(cancellationToken);


284  }


285 


286  // compute multipliers


287 


288  t = 1.0 / col_k[k];


289  dscal(n  (kp1), t, col_k, kp1, 1);


290 


291  if (cancellationToken.IsCancellationRequested) {


292  throw new OperationCanceledException(cancellationToken);


293  }


294 


295  // row elimination with column indexing


296 


297  for (j = kp1; j < n; j++) {


298  col_j = a[j];


299  t = col_j[l];


300  if (l != k) {


301  col_j[l] = col_j[k];


302  col_j[k] = t;


303  }


304  daxpy(n  (kp1), t, col_k, kp1, 1,


305  col_j, kp1, 1);


306  }


307  } else {


308  info = k;


309  }


310  }


311  }


312 


313  ipvt[n  1] = n  1;


314  if (a[(n  1)][(n  1)] == 0) info = n  1;


315 


316  return info;


317  }


318 


319  private void dgesl(double[][] a, int lda, int n, int[] ipvt, double[] b, int job) {


320  double t;


321  int k, kb, l, nm1, kp1;


322 


323  if (cancellationToken.IsCancellationRequested) {


324  throw new OperationCanceledException(cancellationToken);


325  }


326 


327  nm1 = n  1;


328  if (job == 0) {


329  // job = 0 , solve a * x = b. first solve l*y = b


330 


331  if (nm1 >= 1) {


332  for (k = 0; k < nm1; k++) {


333  l = ipvt[k];


334  t = b[l];


335  if (l != k) {


336  b[l] = b[k];


337  b[k] = t;


338  }


339  kp1 = k + 1;


340  daxpy(n  (kp1), t, a[k], kp1, 1, b, kp1, 1);


341  }


342  }


343 


344  if (cancellationToken.IsCancellationRequested) {


345  throw new OperationCanceledException(cancellationToken);


346  }


347 


348  // now solve u*x = y


349 


350  for (kb = 0; kb < n; kb++) {


351  k = n  (kb + 1);


352  b[k] /= a[k][k];


353  t = b[k];


354  daxpy(k, t, a[k], 0, 1, b, 0, 1);


355  }


356  } else {


357  // job = nonzero, solve trans(a) * x = b. first solve trans(u)*y = b


358 


359  for (k = 0; k < n; k++) {


360  t = ddot(k, a[k], 0, 1, b, 0, 1);


361  b[k] = (b[k]  t) / a[k][k];


362  }


363 


364  if (cancellationToken.IsCancellationRequested) {


365  throw new OperationCanceledException(cancellationToken);


366  }


367 


368  // now solve trans(l)*x = y


369 


370  if (nm1 >= 1) {


371  //for (kb = 1; kb < nm1; kb++) {


372  for (kb = 0; kb < nm1; kb++) {


373  k = n  (kb + 1);


374  kp1 = k + 1;


375  b[k] += ddot(n  (kp1), a[k], kp1, 1, b, kp1, 1);


376  l = ipvt[k];


377  if (l != k) {


378  t = b[l];


379  b[l] = b[k];


380  b[k] = t;


381  }


382  }


383  }


384  }


385  }


386 


387  private void daxpy(int n, double da, double[] dx, int dx_off, int incx, double[] dy, int dy_off, int incy) {


388  int i, ix, iy;


389 


390  if (cancellationToken.IsCancellationRequested) {


391  throw new OperationCanceledException(cancellationToken);


392  }


393 


394  if ((n > 0) && (da != 0)) {


395  if (incx != 1  incy != 1) {


396 


397  // code for unequal increments or equal increments not equal to 1


398 


399  ix = 0;


400  iy = 0;


401  if (incx < 0) ix = (n + 1) * incx;


402  if (incy < 0) iy = (n + 1) * incy;


403  for (i = 0; i < n; i++) {


404  dy[iy + dy_off] += da * dx[ix + dx_off];


405  ix += incx;


406  iy += incy;


407  }


408  return;


409  } else {


410  // code for both increments equal to 1


411 


412  for (i = 0; i < n; i++)


413  dy[i + dy_off] += da * dx[i + dx_off];


414  }


415  }


416  }


417 


418  private double ddot(int n, double[] dx, int dx_off, int incx, double[] dy, int dy_off, int incy) {


419  double dtemp = 0;


420  int i, ix, iy;


421 


422  if (cancellationToken.IsCancellationRequested) {


423  throw new OperationCanceledException(cancellationToken);


424  }


425 


426  if (n > 0) {


427  if (incx != 1  incy != 1) {


428  // code for unequal increments or equal increments not equal to 1


429 


430  ix = 0;


431  iy = 0;


432  if (incx < 0) ix = (n + 1) * incx;


433  if (incy < 0) iy = (n + 1) * incy;


434  for (i = 0; i < n; i++) {


435  dtemp += dx[ix + dx_off] * dy[iy + dy_off];


436  ix += incx;


437  iy += incy;


438  }


439  } else {


440  // code for both increments equal to 1


441 


442  for (i = 0; i < n; i++)


443  dtemp += dx[i + dx_off] * dy[i + dy_off];


444  }


445  }


446  return (dtemp);


447  }


448 


449  private void dscal(int n, double da, double[] dx, int dx_off, int incx) {


450  int i, nincx;


451 


452  if (cancellationToken.IsCancellationRequested) {


453  throw new OperationCanceledException(cancellationToken);


454  }


455 


456  if (n > 0) {


457  if (incx != 1) {


458  // code for increment not equal to 1


459 


460  nincx = n * incx;


461  for (i = 0; i < nincx; i += incx)


462  dx[i + dx_off] *= da;


463  } else {


464  // code for increment equal to 1


465 


466  for (i = 0; i < n; i++)


467  dx[i + dx_off] *= da;


468  }


469  }


470  }


471 


472  private int idamax(int n, double[] dx, int dx_off, int incx) {


473  double dmax, dtemp;


474  int i, ix, itemp = 0;


475 


476  if (cancellationToken.IsCancellationRequested) {


477  throw new OperationCanceledException(cancellationToken);


478  }


479 


480  if (n < 1) {


481  itemp = 1;


482  } else if (n == 1) {


483  itemp = 0;


484  } else if (incx != 1) {


485  // code for increment not equal to 1


486 


487  dmax = (dx[dx_off] < 0.0) ? dx[dx_off] : dx[dx_off];


488  ix = 1 + incx;


489  for (i = 0; i < n; i++) {


490  dtemp = (dx[ix + dx_off] < 0.0) ? dx[ix + dx_off] : dx[ix + dx_off];


491  if (dtemp > dmax) {


492  itemp = i;


493  dmax = dtemp;


494  }


495  ix += incx;


496  }


497  } else {


498  // code for increment equal to 1


499 


500  itemp = 0;


501  dmax = (dx[dx_off] < 0.0) ? dx[dx_off] : dx[dx_off];


502  for (i = 0; i < n; i++) {


503  dtemp = (dx[i + dx_off] < 0.0) ? dx[i + dx_off] : dx[i + dx_off];


504  if (dtemp > dmax) {


505  itemp = i;


506  dmax = dtemp;


507  }


508  }


509  }


510  return (itemp);


511  }


512 


513  private double epslon(double x) {


514  double a, b, c, eps;


515 


516  a = 4.0e0 / 3.0e0;


517  eps = 0;


518  while (eps == 0) {


519  b = a  1.0;


520  c = b + b + b;


521  eps = abs(c  1.0);


522  }


523  return (eps * abs(x));


524  }


525 


526  private void dmxpy(int n1, double[] y, int n2, int ldm, double[] x, double[][] m) {


527  int j, i;


528 


529  // cleanup odd vector


530  for (j = 0; j < n2; j++) {


531  for (i = 0; i < n1; i++) {


532  y[i] += x[j] * m[j][i];


533  }


534  }


535  }


536 


537  #endregion


538 


539  #region Clone


540 


541  public IDeepCloneable Clone(Cloner cloner) {


542  return new LinpackBenchmark(this, cloner);


543  }


544 


545  public object Clone() {


546  return Clone(new Cloner());


547  }


548 


549  #endregion


550 


551  #region Events


552 


553  public event EventHandler ItemImageChanged;


554  protected virtual void OnItemImageChanged() {


555  EventHandler handler = ItemImageChanged;


556  if (handler != null) handler(this, EventArgs.Empty);


557  }


558 


559  public event EventHandler ToStringChanged;


560  protected virtual void OnToStringChanged() {


561  EventHandler handler = ToStringChanged;


562  if (handler != null) handler(this, EventArgs.Empty);


563  }


564 


565  #endregion


566  }


567  }

