#include <stdio.h>
#include <stdlib.h>
#include <sys/times.h>
#include <unistd.h>
#include <math.h>
#ifndef MINDIM_M
#define MINDIM_M 1000
#endif
#ifndef MINDIM_N
#define MINDIM_N 1000
#endif
#ifndef MINDIM_K
#define MINDIM_K 1000
#endif
#ifndef MAXDIM_M
#define MAXDIM_M 5000
#endif
#ifndef MAXDIM_N
#define MAXDIM_N 5000
#endif
#ifndef MAXDIM_K
#define MAXDIM_K 5000
#endif
#ifndef INC_M
#define INC_M 100
#endif
#ifndef INC_N
#define INC_N 100
#endif
#ifndef INC_K
#define INC_K 100
#endif
#ifndef MIN_T
#define MIN_T 1
#endif
double buffer[MAXDIM_M*MAXDIM_N];
double buffer1[MAXDIM_M*MAXDIM_N];
double buffer2[MAXDIM_M*MAXDIM_N];
double
walltime()
{
struct tms ts;
static double ClockTick=0.0;
if (ClockTick==0.0) {
ClockTick = 1.0 / ((double) sysconf(_SC_CLK_TCK));
}
return ((double) times(&ts)) * ClockTick;
}
void
gecopy(long m, long n, const double *A, long incRowA, long incColA,
double *B, long incRowB, long incColB)
{
long i, j;
for (i=0; i<m; ++i) {
for (j=0; j<n; ++j) {
B[i*incRowB+j*incColB] = A[i*incRowA+j*incColA];
}
}
}
double
asumDiff(long m, long n, const double *A, long incRowA, long incColA,
const double *B, long incRowB, long incColB)
{
double diff = 0;
long i, j;
for (i=0; i<m; ++i) {
for (j=0; j<n; ++j) {
diff += fabs(B[i*incRowB+j*incColB] - A[i*incRowA+j*incColA]);
}
}
return diff;
}
//------------------------------------------------------------------------------
void
initMatrix(long m, long n, double *A, long incRowA, long incColA)
{
long i, j;
for (j=0; j<n; ++j) {
for (i=0; i<m; ++i) {
A[i*incRowA+j*incColA] = j*n+i+1;
}
}
}
void
printMatrix(long m, long n, const double *A, long incRowA, long incColA)
{
long i, j;
for (i=0; i<m; ++i) {
printf(" ");
for (j=0; j<n; ++j) {
printf("%4.1lf ", A[i*incRowA+j*incColA]);
}
printf("\n");
}
printf("\n");
}
int
main()
{
long m, n, runs;
double t0, t1, t2, diff;
printf("# M N t1 (colMajor) t2 (rowMajor) t2/t1 diff\n");
printf("#================================================================\n");
for (m=MINDIM_M, n=MINDIM_N; m<=MAXDIM_M && n<MAXDIM_N; m+=INC_M, n+=INC_N) {
printf("%4ld %4ld ", m, n);
runs = 0;
t1 = 0;
do {
t0 = walltime();
initMatrix(m, n, buffer, 1, m);
t1 += walltime() - t0;
++runs;
} while (t1<MIN_T);
t1 /= runs;
gecopy(m, n, buffer, 1, m, buffer1, 1, m);
printf(" %12.2lf", t1);
runs = 0;
t2 = 0;
do {
t0 = walltime();
initMatrix(m, n, buffer, n, 1);
t2 += walltime() - t0;
++runs;
} while (t2<MIN_T);
t2 /= runs;
gecopy(m, n, buffer, n, 1, buffer2, n, 1);
diff = asumDiff(m, n, buffer1, 1, m, buffer2, n, 1);
printf(" %12.2lf %16.2lf", t2, t2/t1);
printf(" %11.2le", diff);
printf("\n");
}
return 0;
}