SingleSource/Benchmarks/CoyoteBench/lpbench.c - third_party/llvm-test-suite - Git at Google

 /*
     lpbench
     Standard C version
     17 April 2003

     Written by Scott Robert Ladd (scott@coyotegulch.com)
     No rights reserved. This is public domain software, for use by anyone.

     A number-crunching benchmark that can be used as a fitness test for
     evolving optimal compiler options via genetic algorithm.

     This program is a modernized version of the classic Linpack
     benchmark. My implementation is based on Bonnie Toy's translation
     of the original Fortran source code, including modifications by Jack
     Dongarra and Roy Longbottom.

     I've updated the code as follows:

       1) ANSI C code
       2) A timing mechanism based on the POSIX clock_gettime() function
       3) Replaces C's rather weak random number generator
       4) Compute results for a 2500x2500 (6.25 million element) matrix

     Note that the code herein is design for the purpose of testing
     computational performance; error handling and other such "niceties"
     is virtually non-existent.

     Actual benchmark results can be found at:
             http://www.coyotegulch.com

     Please do not use this information or algorithm in any way that might
     upset the balance of the universe or otherwise cause a disturbance in
     the space-time continuum.
 */

 #include <time.h>
 #include <string.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <stdbool.h>

 // embedded random number generator; ala Park and Miller
 static       long seed = 1325;
 static const long IA   = 16807;
 static const long IM   = 2147483647;
 static const double AM = 4.65661287525E-10;
 static const long IQ   = 127773;
 static const long IR   = 2836;
 static const long MASK = 123459876;

 static double random_double()
 {
     long k;
     double result;

     seed ^= MASK;
     k = seed / IQ;
     seed = IA * (seed - k * IQ) - IR * k;

     if (seed < 0)
         seed += IM;

     result = AM * seed;
     seed ^= MASK;

     return result;
 }

 #if defined(ACOVEA)
 #if defined(arch_pentium4)
 static const int N   = 1000;
 static const int NM1 =  999; // N - 1
 static const int NP1 = 1001; // N + 1
 #else
 static const int N   =  600;
 static const int NM1 =  599; // N - 1
 static const int NP1 =  601; // N + 1
 #endif
 #else
 #ifdef SMALL_PROBLEM_SIZE
 static const int N   = 400;
 static const int NM1 = 399; // N - 1
 static const int NP1 = 401; // N + 1
 #else
 static const int N   = 2000;
 static const int NM1 = 1999; // N - 1
 static const int NP1 = 2001; // N + 1
 #endif
 #endif

 // benchmark code
 void matgen (double ** a, double * b)
 {
     // fill arrays
     int i, j;

     for (i = 0; i < N; i++)
         for (j = 0; j < N; j++)
             a[j][i] = random_double();

     for (i = 0; i < N; i++)
         b[i] = 0.0;

     for (j = 0; j < N; j++)
         for (i = 0; i < N; i++)
             b[i] += a[j][i];
 }

 // finds the index of element having max. absolute value.
 int idamax(int n, double * dx, int dx_off, int incx)
 {
     double dmax, dtemp;
     int i, ix, itemp = 0;

     if (n < 1)
         itemp = -1;
     else
     {
         if (n ==1)
             itemp = 0;
         else
         {
             if (incx != 1)
             {
                 // code for increment not equal to 1
                 dmax = fabs(dx[dx_off]);
                 ix = 1 + incx;

                 for (i = 1; i < n; i++)
                 {
 	                dtemp = fabs(dx[ix + dx_off]);

                     if (dtemp > dmax)
                     {
 	                    itemp = i;
 	                    dmax = dtemp;
                     }

                     ix += incx;
                 }
             }
             else
             {
                 // code for increment equal to 1
                 itemp = 0;
                 dmax = fabs(dx[dx_off]);

                 for (i = 1; i < n; i++)
                 {
 	                dtemp = fabs(dx[i + dx_off]);

 	                if (dtemp > dmax)
                     {
                         itemp = i;
                         dmax = dtemp;
                     }
                 }
             }
         }
     }

     return itemp;
 }

 // forms the dot product of two vectors.
 static double ddot(int n, double * dx, int dx_off, int incx, double * dy, int dy_off, int incy)
 {
     int i;
     double dtemp = 0.0;

     if (n > 0)
     {
         if (incx != 1 || incy != 1)
         {
             // code for unequal increments or equal increments not equal to 1
             int ix = 0;
 	        int iy = 0;

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

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

             for (i = 0;i < n; i++)
             {
                 dtemp += dx[ix + dx_off] * dy[iy + dy_off];
                 ix    += incx;
                 iy    += incy;
             }
         }
         else
         {
             // code for both increments equal to 1
             for (i=0; i < n; i++)
                 dtemp += dx[i + dx_off] * dy[i + dy_off];
         }
     }

     return(dtemp);
 }

 // scales a vector by a constant.
 void dscal(int n, double da, double * dx, int dx_off, int incx)
 {
     int i;

     if (n > 0)
     {
         if (incx != 1)
         {
             // code for increment not equal to 1
 	        int nincx = n * incx;

             for (i = 0; i < nincx; i += incx)
                 dx[i + dx_off] *= da;
         }
         else
         {
             // code for increment equal to 1
             for (i = 0; i < n; i++)
                 dx[i + dx_off] *= da;
         }
     }
 }

 //  constant times a vector plus a vector.
 void daxpy(int n, double da, double * dx, int dx_off, int incx, double * dy, int dy_off, int incy)
 {
     int i;

     if ((n > 0) && (da != 0))
     {
         if (incx != 1 || incy != 1)
         {
             // code for unequal increments or equal increments not equal to 1
             int ix = 0;
             int iy = 0;

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

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

             for (i = 0; i < n; i++)
             {
                 dy[iy + dy_off] += da * dx[ix + dx_off];
                 ix += incx;
                 iy += incy;
 	        }

 	        return;
         }
         else
         {
             // code for both increments equal to 1
             for (i = 0; i < n; i++)
                 dy[i + dy_off] += da * dx[i + dx_off];
         }
     }
 }

 // Factors a double precision matrix by gaussian elimination.
 void dgefa(double ** a, int * ipvt)
 {
     double temp;
     int k, j;

     for (k = 0; k < NM1; k++)
     {
         double * col_k = a[k];
         int kp1 = k + 1;

         // find l = pivot index
         int l = idamax(N - k, col_k, k, 1) + k;
 	    ipvt[k] = l;

     	// zero pivot implies this column already triangularized
         if (col_k[l] != 0)
         {
 	        // interchange if necessary
     	    if (l != k)
             {
     	        temp     = col_k[l];
 	            col_k[l] = col_k[k];
 	            col_k[k] = temp;
             }

             // compute multipliers
             temp = -1.0 / col_k[k];
             dscal(N - kp1, temp, col_k, kp1, 1);

     	    // row elimination with column indexing
             for (j = kp1; j < N; j++)
             {
 	            double * col_j = a[j];
 	            temp = col_j[l];

     	        if (l != k)
                 {
                     col_j[l] = col_j[k];
                     col_j[k] = temp;
                 }

                 daxpy(N - kp1, temp, col_k, kp1, 1, col_j, kp1, 1);
             }
         }
     }

     ipvt[N - 1] = N - 1;
 }

 //  Solves the double precision system a * x = b  or trans(a) * x = b
 //  using the factors computed by dgeco or dgefa.
 void dgesl(double ** a, int * ipvt, double * b)
 {
     double t;
     int k, kb;

     // solve  a * x = b.  first solve  l*y = b
     for (k = 0; k < NM1; k++)
     {
         int l = ipvt[k];
         t = b[l];

         if (l != k)
         {
             b[l] = b[k];
             b[k] = t;
         }

         int kp1 = k + 1;
         daxpy(N - kp1, t, a[k], kp1, 1, b, kp1, 1);
     }

     // now solve  u*x = y
     for (kb = 0; kb < N; kb++)
     {
         k     = N - (kb + 1);
         b[k] /= a[k][k];
         t     = -b[k];
         daxpy(k, t, a[k], 0, 1, b, 0, 1);
     }
 }

 int main(int argc, char ** argv)
 {
     int i;

     // do we have verbose output?
     bool ga_testing = false;

     if (argc > 1)
     {
         for (i = 1; i < argc; ++i)
         {
             if (!strcmp(argv[1],"-ga"))
             {
                 ga_testing = true;
                 break;
             }
         }
     }

     double ** a = (double **)malloc(sizeof(double) * N);

     for (i = 0; i < N; ++i)
         a[i] = (double *)malloc(sizeof(double) * NP1);

     double * b = (double *)malloc(sizeof(double) * N);
     double * x = (double *)malloc(sizeof(double) * N);
     int * ipvt = (int *)malloc(sizeof(int)    * N);

     // calculate operations per timeInSeconds
     double ops = (2.0 * ((double)N * N * N)) / 3.0 + 2.0 * ((double)N * N);

     // generate matrix
     matgen(a,b);

     // get starting time
     //struct timespec start, stop;
     //clock_gettime(CLOCK_REALTIME,&start);

     // what we're timing
     dgefa(a,ipvt);
     dgesl(a,ipvt,b);

     // calculate run time
     //clock_gettime(CLOCK_REALTIME,&stop);
     double run_time = 0;//(stop.tv_sec - start.tv_sec) + (double)(stop.tv_nsec - start.tv_nsec) / 1000000000.0;

     // clean up
     free(ipvt);
     free(x);
     free(b);

     for (i = 0; i < N; ++i)
         free(a[i]);

     free(a);

     // report runtime
     if (ga_testing)
         fprintf(stdout,"%f",run_time);
     else
         fprintf(stdout,"\nlpbench (Std. C) run time: %f\n\n",run_time);

     fflush(stdout);

     // done
     return 0;
 }
	/*
	lpbench
	Standard C version
	17 April 2003

	Written by Scott Robert Ladd (scott@coyotegulch.com)
	No rights reserved. This is public domain software, for use by anyone.

	A number-crunching benchmark that can be used as a fitness test for
	evolving optimal compiler options via genetic algorithm.

	This program is a modernized version of the classic Linpack
	benchmark. My implementation is based on Bonnie Toy's translation
	of the original Fortran source code, including modifications by Jack
	Dongarra and Roy Longbottom.

	I've updated the code as follows:

	1) ANSI C code
	2) A timing mechanism based on the POSIX clock_gettime() function
	3) Replaces C's rather weak random number generator
	4) Compute results for a 2500x2500 (6.25 million element) matrix

	Note that the code herein is design for the purpose of testing
	computational performance; error handling and other such "niceties"
	is virtually non-existent.

	Actual benchmark results can be found at:
	http://www.coyotegulch.com

	Please do not use this information or algorithm in any way that might
	upset the balance of the universe or otherwise cause a disturbance in
	the space-time continuum.
	*/

	#include <time.h>
	#include <string.h>
	#include <stdio.h>
	#include <stdlib.h>
	#include <math.h>
	#include <stdbool.h>

	// embedded random number generator; ala Park and Miller
	static long seed = 1325;
	static const long IA = 16807;
	static const long IM = 2147483647;
	static const double AM = 4.65661287525E-10;
	static const long IQ = 127773;
	static const long IR = 2836;
	static const long MASK = 123459876;

	static double random_double()
	{
	long k;
	double result;

	seed ^= MASK;
	k = seed / IQ;
	seed = IA * (seed - k * IQ) - IR * k;

	if (seed < 0)
	seed += IM;

	result = AM * seed;
	seed ^= MASK;

	return result;
	}

	#if defined(ACOVEA)
	#if defined(arch_pentium4)
	static const int N = 1000;
	static const int NM1 = 999; // N - 1
	static const int NP1 = 1001; // N + 1
	#else
	static const int N = 600;
	static const int NM1 = 599; // N - 1
	static const int NP1 = 601; // N + 1
	#endif
	#else
	#ifdef SMALL_PROBLEM_SIZE
	static const int N = 400;
	static const int NM1 = 399; // N - 1
	static const int NP1 = 401; // N + 1
	#else
	static const int N = 2000;
	static const int NM1 = 1999; // N - 1
	static const int NP1 = 2001; // N + 1
	#endif
	#endif

	// benchmark code
	void matgen (double ** a, double * b)
	{
	// fill arrays
	int i, j;

	for (i = 0; i < N; i++)
	for (j = 0; j < N; j++)
	a[j][i] = random_double();

	for (i = 0; i < N; i++)
	b[i] = 0.0;

	for (j = 0; j < N; j++)
	for (i = 0; i < N; i++)
	b[i] += a[j][i];
	}

	// finds the index of element having max. absolute value.
	int idamax(int n, double * dx, int dx_off, int incx)
	{
	double dmax, dtemp;
	int i, ix, itemp = 0;

	if (n < 1)
	itemp = -1;
	else
	{
	if (n ==1)
	itemp = 0;
	else
	{
	if (incx != 1)
	{
	// code for increment not equal to 1
	dmax = fabs(dx[dx_off]);
	ix = 1 + incx;

	for (i = 1; i < n; i++)
	{
	dtemp = fabs(dx[ix + dx_off]);

	if (dtemp > dmax)
	{
	itemp = i;
	dmax = dtemp;
	}

	ix += incx;
	}
	}
	else
	{
	// code for increment equal to 1
	itemp = 0;
	dmax = fabs(dx[dx_off]);

	for (i = 1; i < n; i++)
	{
	dtemp = fabs(dx[i + dx_off]);

	if (dtemp > dmax)
	{
	itemp = i;
	dmax = dtemp;
	}
	}
	}
	}
	}

	return itemp;
	}

	// forms the dot product of two vectors.
	static double ddot(int n, double * dx, int dx_off, int incx, double * dy, int dy_off, int incy)
	{
	int i;
	double dtemp = 0.0;

	if (n > 0)
	{
	if (incx != 1 \|\| incy != 1)
	{
	// code for unequal increments or equal increments not equal to 1
	int ix = 0;
	int iy = 0;

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

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

	for (i = 0;i < n; i++)
	{
	dtemp += dx[ix + dx_off] * dy[iy + dy_off];
	ix += incx;
	iy += incy;
	}
	}
	else
	{
	// code for both increments equal to 1
	for (i=0; i < n; i++)
	dtemp += dx[i + dx_off] * dy[i + dy_off];
	}
	}

	return(dtemp);
	}

	// scales a vector by a constant.
	void dscal(int n, double da, double * dx, int dx_off, int incx)
	{
	int i;

	if (n > 0)
	{
	if (incx != 1)
	{
	// code for increment not equal to 1
	int nincx = n * incx;

	for (i = 0; i < nincx; i += incx)
	dx[i + dx_off] *= da;
	}
	else
	{
	// code for increment equal to 1
	for (i = 0; i < n; i++)
	dx[i + dx_off] *= da;
	}
	}
	}

	// constant times a vector plus a vector.
	void daxpy(int n, double da, double * dx, int dx_off, int incx, double * dy, int dy_off, int incy)
	{
	int i;

	if ((n > 0) && (da != 0))
	{
	if (incx != 1 \|\| incy != 1)
	{
	// code for unequal increments or equal increments not equal to 1
	int ix = 0;
	int iy = 0;

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

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

	for (i = 0; i < n; i++)
	{
	dy[iy + dy_off] += da * dx[ix + dx_off];
	ix += incx;
	iy += incy;
	}

	return;
	}
	else
	{
	// code for both increments equal to 1
	for (i = 0; i < n; i++)
	dy[i + dy_off] += da * dx[i + dx_off];
	}
	}
	}

	// Factors a double precision matrix by gaussian elimination.
	void dgefa(double ** a, int * ipvt)
	{
	double temp;
	int k, j;

	for (k = 0; k < NM1; k++)
	{
	double * col_k = a[k];
	int kp1 = k + 1;

	// find l = pivot index
	int l = idamax(N - k, col_k, k, 1) + k;
	ipvt[k] = l;

	// zero pivot implies this column already triangularized
	if (col_k[l] != 0)
	{
	// interchange if necessary
	if (l != k)
	{
	temp = col_k[l];
	col_k[l] = col_k[k];
	col_k[k] = temp;
	}

	// compute multipliers
	temp = -1.0 / col_k[k];
	dscal(N - kp1, temp, col_k, kp1, 1);

	// row elimination with column indexing
	for (j = kp1; j < N; j++)
	{
	double * col_j = a[j];
	temp = col_j[l];

	if (l != k)
	{
	col_j[l] = col_j[k];
	col_j[k] = temp;
	}

	daxpy(N - kp1, temp, col_k, kp1, 1, col_j, kp1, 1);
	}
	}
	}

	ipvt[N - 1] = N - 1;
	}

	// Solves the double precision system a * x = b or trans(a) * x = b
	// using the factors computed by dgeco or dgefa.
	void dgesl(double ** a, int * ipvt, double * b)
	{
	double t;
	int k, kb;

	// solve a * x = b. first solve l*y = b
	for (k = 0; k < NM1; k++)
	{
	int l = ipvt[k];
	t = b[l];

	if (l != k)
	{
	b[l] = b[k];
	b[k] = t;
	}

	int kp1 = k + 1;
	daxpy(N - kp1, t, a[k], kp1, 1, b, kp1, 1);
	}

	// now solve u*x = y
	for (kb = 0; kb < N; kb++)
	{
	k = N - (kb + 1);
	b[k] /= a[k][k];
	t = -b[k];
	daxpy(k, t, a[k], 0, 1, b, 0, 1);
	}
	}

	int main(int argc, char ** argv)
	{
	int i;

	// do we have verbose output?
	bool ga_testing = false;

	if (argc > 1)
	{
	for (i = 1; i < argc; ++i)
	{
	if (!strcmp(argv[1],"-ga"))
	{
	ga_testing = true;
	break;
	}
	}
	}

	double a = (double )malloc(sizeof(double) * N);

	for (i = 0; i < N; ++i)
	a[i] = (double )malloc(sizeof(double) NP1);

	double * b = (double )malloc(sizeof(double) N);
	double * x = (double )malloc(sizeof(double) N);
	int * ipvt = (int )malloc(sizeof(int) N);

	// calculate operations per timeInSeconds
	double ops = (2.0 * ((double)N * N * N)) / 3.0 + 2.0 * ((double)N * N);

	// generate matrix
	matgen(a,b);

	// get starting time
	//struct timespec start, stop;
	//clock_gettime(CLOCK_REALTIME,&start);

	// what we're timing
	dgefa(a,ipvt);
	dgesl(a,ipvt,b);

	// calculate run time
	//clock_gettime(CLOCK_REALTIME,&stop);
	double run_time = 0;//(stop.tv_sec - start.tv_sec) + (double)(stop.tv_nsec - start.tv_nsec) / 1000000000.0;

	// clean up
	free(ipvt);
	free(x);
	free(b);

	for (i = 0; i < N; ++i)
	free(a[i]);

	free(a);

	// report runtime
	if (ga_testing)
	fprintf(stdout,"%f",run_time);
	else
	fprintf(stdout,"\nlpbench (Std. C) run time: %f\n\n",run_time);

	fflush(stdout);

	// done
	return 0;
	}