#include <stdio.h>
#include <sys/time.h>
#include "nrutil.h"
#include <math.h>
double mult(int n,double **M1,double **M2,double **M3,
double *v1,double *v2,double *v3)
{
int i;
double x=1.;
for (i=0; i<n*n*n; i++)
x=x*1.000001;
return x;
}
double timer(char *name,double flop,
double (*func)(int,double **,double **,double **,
double *,double *,double *),
int n)
{
/* Timer variables */
struct timeval time1,time2;
double ctime,mflop;
double **M1,*v1;
double **M2,*v2;
double **M3,*v3;
int i,j;
M1=dmatrix(0,n,0,n);
M2=dmatrix(0,n,0,n);
M3=dmatrix(0,n,0,n);
v1=dvector(0,n);
v2=dvector(0,n);
v3=dvector(0,n);
/* Fill in values for input vectors */
for (i=0; i<n; i++) {
v1[i]=3.14159;
v2[i]=0.314159;
}
/* Fill in values for input matrices */
for (i=0; i<n; i++)
for (j=0; j<n; j++) {
M1[i][j]=1./(1+i+j);
M2[i][j]=1./(2+i+j);
}
/* Record time of function call */
gettimeofday(&time1,NULL);
/* Call function */
func(n,M1,M2,M3,v1,v2,v3);
/* Record time of function return */
gettimeofday(&time2,NULL);
/* Coupute time in seconds */
ctime=
(time2.tv_sec-time1.tv_sec) /* Integer seconds */
+(time2.tv_usec-time1.tv_usec)/1.e6; /* Integer giving microseconds part */
/* Convert from FLOP to MFLOP */
mflop=flop/1e6;
/* Output result */
fprintf(stderr,"%20s: %4.0lf MFLOPS = %.1lf MFLOP / %.3f sec\n",
name,mflop/ctime,mflop,ctime);
/* Free memory before return */
free_dvector(v1,0,n);
free_dvector(v2,0,n);
free_dvector(v3,0,n);
free_dmatrix(M1,0,n,0,n);
free_dmatrix(M2,0,n,0,n);
free_dmatrix(M3,0,n,0,n);
return ctime;
}
main()
{
/* Place to store run times of various routines */
double tm,td,ts,ti,tal,tac1,tac2;
/* Variables for verification of matmat3 */
int i,j;
double **M1,**M2,**M3,**M4;
double err=0.;
/* Problem size (n x n) */
int n;
printf("Problem size: ");
scanf("%d",&n);
/* Calls to timers which don't access much memory */
tm=timer("mult",(double) n*n*n,mult,n);
/* td=timer("divide",(double) n*n*n,divide,n); */
/* ts=timer("sub ",(double) n*n*n,multsub,n); */
/* ti=timer("if ",(double) n*n*n,multif,n); */
/* tal=timer("alloc",(double) n*n*n,multal,n); */
/* timer("multadd",(double) 2*n*n*n,multadd,n); */
/* Time strides through memory */
/* tac1=timer("cmat1",(double) n*n*n,cmat1,n); */
/* tac2=timer("cmat2",(double) n*n*n,cmat2,n); */
/* Print out time ratios */
/* printf("Notes:\n"); */
/* printf("Division costs %4.1f mults\n",td/tm); */
/* printf("Subroutine costs %4.1f mults\n",ts/tm-1); */
/* printf("If costs %4.1f mults\n",ti/tm-1); */
/* printf("malloc costs %4.1f mults\n",tal/tm-1); */
/* printf("Memory access costs %4.1f mults\n",tac1/tm-1); */
/* printf("Cache miss costs %4.1f accesses\n",tac2/tac1); */
/* printf("\n"); */
/* Matrix-vector operations */
/* timer("vecmat1",(double) 2*n*n*n,vecmat1,n); */
/* timer("vecmat2",(double) 2*n*n*n,vecmat2,n); */
/* timer("matmat1",(double) 2*n*n*n,matmat1,n); */
/* timer("matmat2",(double) 2*n*n*n,matmat2,n); */
/* Verify matmat2 */
/* M1=dmatrix(0,n,0,n); */
/* M2=dmatrix(0,n,0,n); */
/* M3=dmatrix(0,n,0,n); */
/* M4=dmatrix(0,n,0,n); */
/* Fill in values for input matrices */
/* for (i=0; i<n; i++) */
/* for (j=0; j<n; j++) { */
/* M1[i][j]=1./(1+i+j); */
/* M2[i][j]=1./(2+i+j); */
/* } */
/* Note, the matmat routines don't use their last three arguments anyway */
/* matmat1(n,M1,M2,M3,M1[0],M2[0],M3[0]); */
/* matmat2(n,M1,M2,M4,M1[0],M2[0],M3[0]); */
/* for (i=0; i<n; i++) */
/* for (j=0; j<n; j++) */
/* err+=pow(M4[i][j]-M3[i][j],2); */
/* printf("\nIMPORTANT: Test of matmat2; should be near zero... %e\n\n",err); */
/* Free up space */
/* free_dmatrix(M1,0,n,0,n); */
/* free_dmatrix(M2,0,n,0,n); */
/* free_dmatrix(M3,0,n,0,n); */
/* free_dmatrix(M4,0,n,0,n); */
}