#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/times.h>
#include <unistd.h>
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
initGeMatrix(size_t m, size_t n,
double *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
A[i*incRowA + j*incColA] = i*n + j +1;
}
}
}
void
dgeaxpy(size_t m, size_t n,
double alpha,
const double *A,
ptrdiff_t incRowA, ptrdiff_t incColA,
double *B,
ptrdiff_t incRowB, ptrdiff_t incColB)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
B[i*incRowB + j*incColB] += alpha*A[i*incRowA + j*incColA];
}
}
}
#define MAX_M 2*16*16*64
#define MAX_N 2*16*16*64
int
main()
{
double *A = malloc(MAX_M*MAX_N*sizeof(double));
double *B = malloc(MAX_M*MAX_N*sizeof(double));
printf("#%5s %5s %10s %10s\n", "m", "n", "t (colm)", "t (rowm)");
for (size_t m=128, n=128; m<=MAX_M && n<=MAX_N; m+=16, n+=16) {
printf(" %5ld %5ld ", (long)m, (long)n);
{
ptrdiff_t incRowA = MAX_N;
ptrdiff_t incColA = 1;
ptrdiff_t incRowB = MAX_N;
ptrdiff_t incColB = 1;
double alpha = 3;
initGeMatrix(m, n, A, incRowA, incColA);
double t = walltime();
dgeaxpy(m, n, alpha, A, incRowA, incColA, B, incRowB, incColB);
t = walltime() - t;
printf("%10.2lf ", t);
}
{
ptrdiff_t incRowA = 1;
ptrdiff_t incColA = MAX_M;
ptrdiff_t incRowB = 1;
ptrdiff_t incColB = MAX_M;
double alpha = 3;
initGeMatrix(m, n, A, incRowA, incColA);
double t = walltime();
dgeaxpy(m, n, alpha, A, incRowA, incColA, B, incRowB, incColB);
t = walltime() - t;
printf("%10.2lf ", t);
}
printf("\n");
}
free(B);
free(A);
}