/***************************************************************************/
/* GAUSS Benchmark program version 3.0                                     */
/* Author : Stefan Steinhaus                                               */
/* EMAIL  : stefan@steinhaus-net.de                                        */
/* This program is public domain. Feel free to copy it freely.             */
/*                                                                         */
/* Note: This program will not run under GAUSS 3.2 for DOS but you can     */
/*       run it on any UNIX version if you want. However for running it    */
/*       with GAUSS for UNIX it might be necessary to adjust the           */
/*       max_workspace variable .. for GAUSS 3.5 this is not necessary.    */
/***************************************************************************/

proc BenchFunc1(runs);
  local result,zeit,timing,a,b,i,j;
  result=0;
  for i (1,runs,1);
    timing=0; zeit=hsec;
    for j (1,20,1);
      b=abs(rndn(1000,1000)/10); a=b'; b=reshape(a,500,20000); a=b';
    endfor;
    timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc2(runs);
  local result,zeit,timing,a,b;
  result=0;
  for i (1,runs,1);
    timing=0; a=rndn(3000000,1);
    zeit=hsec; b=sortc(a,1); a=rev(b); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc3(runs);
  local result,zeit,timing,a,b;
  result=0;
  for i (1,runs,1);
    timing=0; b=abs(rndn(600,600)/10);
    zeit=hsec; a=polychar(b); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc4(runs);
  local result,zeit,timing,a,b;
  result=0;
  for i (1,runs,1);
    timing=0; a=rndn(2^22,1);
    zeit=hsec; b=fft(a); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc5(runs);
  local result,zeit,timing,q,r,s,a,b,c,d,a1,a2,a3,a4,a5,a6;
  result=0;
  for i (1,runs,1);
    timing=0;
    s=seqa(1,1,600); q=rndn(600,600);
    r=seqa(1,1,600); a=rndn(600,600);
    b=seqa(1,1,600); c=a'; d=seqa(1,2,600);
    zeit=hsec; {a1,a2,a3,a4,a5,a6}=qprog(s,q,r,a,b,c,d,0); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc6(runs);
  local result,zeit,timing,a,b,j;
  result=0;
  for i (1,runs,1);
    timing=0;
    for j (1,12,1);
      a=rndn(600,600);
      zeit=hsec; b=det(a); timing=timing+(hsec-zeit)/100;
    endfor;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc8(runs);
  local result,zeit,timing,a,b;
  result=0;
  for i (1,runs,1);
    timing=0; a=rndu(1000,1000);
    zeit=hsec; b=inv(a); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc9(runs);
  local result,zeit,timing,a,b,c;
  result=0;
  for i (1,runs,1);
    timing=0; a=rndn(650,650);
    zeit=hsec; c=eig(a); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc10(runs);
  local result,zeit,timing,a,b;
  result=0;
  for i (1,runs,1);
    timing=0;
    for j (1,20,1);
      a=moment(rndn(700,700),0); zeit=hsec;
      b=chol(a);
      timing=timing+(hsec-zeit)/100;
    endfor;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc11(runs);
  local result,zeit,timing,a,b;
  result=0;
  for i (1,runs,1);
    timing=0; a=rndn(1000,1000);
    zeit=hsec; b=moment(a,0); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc13(runs);
  local result,zeit,timing,a,b,j;
  result=0;
  for i (1,runs,1);
    timing=0; zeit=hsec;
    for j (1,20,1);
      a=rndn(1000,1000)^2;
      b=gamma(a);
    endfor;
    timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc14(runs);
  local result,zeit,timing,a,b,j;
  result=0;
  for i (1,runs,1);
    timing=0; zeit=hsec;
      for j (1,10,1);
        a=rndn(1000,1000)^2;
        b=erf(a); a=erfc(b);
    endfor;
    timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

proc BenchFunc15(runs);
  local result,zeit,timing,a,b;
  result=0;
  for i (1,runs,1);
    timing=0; a=rndn(1000,1000); b=seqa(1,1,1000);
    zeit=hsec; b=olsqr(b,a); timing=(hsec-zeit)/100;
    result=result+timing;
  endfor;
  retp(result/runs);
endp;

/**********************************/
/* Main program                   */
/**********************************/

runs=2;              /* Amount of runs for avg. timing -> suggestion at least 2 */
ver=sysstate(1,0);
timing = {1 48.009 1,1 13.710 1,1 35.631 1,
          1 36.793 1,1 45.365 1,1 40.650 1,
          1 45.075 1,1 48.049 1,1 26.858 1,
          1 25.917 1,1 40.458 1,1 36.362 1,1 32.036 1};
format /RD 8,3;

output file = bench.out reset;
print; print; print;
print "                  !!! GAUSS - Benchmark program !!!";
print "                  ================================="; print; print;
print "Platform: GAUSS " ftos(ver[1],"%*lf",2,0) "." ftos(ver[2],"%*lf",1,0) "." ftos(ver[3],"%*lf",2,0);
print "Displaying the avg. timings over " runs " runs.";
print; print;
print "Function / timing [sec.]                                   result / reference / ratio";
print "----------------------------------------------------------------------------------------";

timing[1,1]=BenchFunc1(runs);
timing[1,3]=timing[1,2]/timing[1,1];
print /flush "20xCreation, trans. & reshaping of a 1000x1000 matrix : " timing[1,1] " " timing[1,2] " " timing[1,3];
timing[2,1]=BenchFunc2(runs);
timing[2,3]=timing[2,2]/timing[2,1];
print /flush "3000000 values sorted ascending and descending_______ : " timing[2,1] " " timing[2,2] " " timing[2,3];
timing[3,1]=BenchFunc3(runs);
timing[3,3]=timing[3,2]/timing[3,1];
print /flush "Characteristic polynom of 600x600 random matrix______ : " timing[3,1] " " timing[3,2] " " timing[3,3];
timing[4,1]=BenchFunc4(runs);
timing[4,3]=timing[4,2]/timing[4,1];
print /flush "FFT over 2^22 values_________________________________ : " timing[4,1] " " timing[4,2] " " timing[4,3];
timing[5,1]=BenchFunc5(runs);
timing[5,3]=timing[5,2]/timing[5,1];
print /flush "QP problem with constraints 600x600 model matrix_____ : " timing[5,1] " " timing[5,2] " " timing[5,3];
timing[6,1]=BenchFunc6(runs);
timing[6,3]=timing[6,2]/timing[6,1];
print /flush "Det. of 12 600x600 random matrices___________________ : " timing[6,1] " " timing[6,2] " " timing[6,3];
timing[7,1]=BenchFunc8(runs);
timing[7,3]=timing[7,2]/timing[7,1];
print /flush "Inverse of a 1000x1000 random matrix_________________ : " timing[7,1] " " timing[7,2] " " timing[7,3];
timing[8,1]=BenchFunc9(runs);
timing[8,3]=timing[8,2]/timing[8,1];
print /flush "Eigenvalues of a 650x650 randommatrix________________ : " timing[8,1] " " timing[8,2] " " timing[8,3];
timing[9,1]=BenchFunc10(runs);
timing[9,3]=timing[9,2]/timing[9,1];
print /flush "20xCholesky decomposition of a 700x700-matrix________ : " timing[9,1] " " timing[9,2] " " timing[9,3];
timing[10,1]=BenchFunc11(runs);
timing[10,3]=timing[10,2]/timing[10,1];
print /flush "1000x1000 cross-product matrix_______________________ : " timing[10,1] " " timing[10,2] " " timing[10,3];
timing[11,1]=BenchFunc13(runs);
timing[11,3]=timing[11,2]/timing[11,1];
print /flush "20xGamma function over a 1000x1000 matrix____________ : " timing[11,1] " " timing[11,2] " " timing[11,3];
timing[12,1]=BenchFunc14(runs);
timing[12,3]=timing[12,2]/timing[12,1];
print /flush "10xGaussian err & inv. err func. of a 1000x1000 matrix: " timing[12,1] " " timing[12,2] " " timing[12,3];
timing[13,1]=BenchFunc15(runs);
timing[13,3]=timing[13,2]/timing[13,1];
print /flush "Linear regression over a 1000x1000 matrix____________ : " timing[13,1] " " timing[13,2] " " timing[13,3];

print "----------------------------------------------------------------------------------------";
print "Overall timing_______________________________________ : " sumc(timing[1:13,1]) " " sumc(timing[1:13,2]) " " sumc(timing[1:13,2])/sumc(timing[1:13,1]);
print;
print "(reference PC: Pentium III, 550 MHz, 384MB, Win. NT 4, GAUSS 3.5.10)";
output off;
end;