1
2 // g++-4.4 bench_gemm.cpp -I .. -O2 -DNDEBUG -lrt -fopenmp && OMP_NUM_THREADS=2 ./a.out
3 // icpc bench_gemm.cpp -I .. -O3 -DNDEBUG -lrt -openmp && OMP_NUM_THREADS=2 ./a.out
4
5 #include <iostream>
6 #include <Eigen/Core>
7 #include <bench/BenchTimer.h>
8
9 using namespace std;
10 using namespace Eigen;
11
12 #ifndef SCALAR
13 // #define SCALAR std::complex<float>
14 #define SCALAR float
15 #endif
16
17 typedef SCALAR Scalar;
18 typedef NumTraits<Scalar>::Real RealScalar;
19 typedef Matrix<RealScalar,Dynamic,Dynamic> A;
20 typedef Matrix</*Real*/Scalar,Dynamic,Dynamic> B;
21 typedef Matrix<Scalar,Dynamic,Dynamic> C;
22 typedef Matrix<RealScalar,Dynamic,Dynamic> M;
23
24 #ifdef HAVE_BLAS
25
26 extern "C" {
27 #include <Eigen/src/misc/blas.h>
28 }
29
30 static float fone = 1;
31 static float fzero = 0;
32 static double done = 1;
33 static double szero = 0;
34 static std::complex<float> cfone = 1;
35 static std::complex<float> cfzero = 0;
36 static std::complex<double> cdone = 1;
37 static std::complex<double> cdzero = 0;
38 static char notrans = 'N';
39 static char trans = 'T';
40 static char nonunit = 'N';
41 static char lower = 'L';
42 static char right = 'R';
43 static int intone = 1;
44
blas_gemm(const MatrixXf & a,const MatrixXf & b,MatrixXf & c)45 void blas_gemm(const MatrixXf& a, const MatrixXf& b, MatrixXf& c)
46 {
47 int M = c.rows(); int N = c.cols(); int K = a.cols();
48 int lda = a.rows(); int ldb = b.rows(); int ldc = c.rows();
49
50 sgemm_(¬rans,¬rans,&M,&N,&K,&fone,
51 const_cast<float*>(a.data()),&lda,
52 const_cast<float*>(b.data()),&ldb,&fone,
53 c.data(),&ldc);
54 }
55
blas_gemm(const MatrixXd & a,const MatrixXd & b,MatrixXd & c)56 EIGEN_DONT_INLINE void blas_gemm(const MatrixXd& a, const MatrixXd& b, MatrixXd& c)
57 {
58 int M = c.rows(); int N = c.cols(); int K = a.cols();
59 int lda = a.rows(); int ldb = b.rows(); int ldc = c.rows();
60
61 dgemm_(¬rans,¬rans,&M,&N,&K,&done,
62 const_cast<double*>(a.data()),&lda,
63 const_cast<double*>(b.data()),&ldb,&done,
64 c.data(),&ldc);
65 }
66
blas_gemm(const MatrixXcf & a,const MatrixXcf & b,MatrixXcf & c)67 void blas_gemm(const MatrixXcf& a, const MatrixXcf& b, MatrixXcf& c)
68 {
69 int M = c.rows(); int N = c.cols(); int K = a.cols();
70 int lda = a.rows(); int ldb = b.rows(); int ldc = c.rows();
71
72 cgemm_(¬rans,¬rans,&M,&N,&K,(float*)&cfone,
73 const_cast<float*>((const float*)a.data()),&lda,
74 const_cast<float*>((const float*)b.data()),&ldb,(float*)&cfone,
75 (float*)c.data(),&ldc);
76 }
77
blas_gemm(const MatrixXcd & a,const MatrixXcd & b,MatrixXcd & c)78 void blas_gemm(const MatrixXcd& a, const MatrixXcd& b, MatrixXcd& c)
79 {
80 int M = c.rows(); int N = c.cols(); int K = a.cols();
81 int lda = a.rows(); int ldb = b.rows(); int ldc = c.rows();
82
83 zgemm_(¬rans,¬rans,&M,&N,&K,(double*)&cdone,
84 const_cast<double*>((const double*)a.data()),&lda,
85 const_cast<double*>((const double*)b.data()),&ldb,(double*)&cdone,
86 (double*)c.data(),&ldc);
87 }
88
89
90
91 #endif
92
matlab_cplx_cplx(const M & ar,const M & ai,const M & br,const M & bi,M & cr,M & ci)93 void matlab_cplx_cplx(const M& ar, const M& ai, const M& br, const M& bi, M& cr, M& ci)
94 {
95 cr.noalias() += ar * br;
96 cr.noalias() -= ai * bi;
97 ci.noalias() += ar * bi;
98 ci.noalias() += ai * br;
99 }
100
matlab_real_cplx(const M & a,const M & br,const M & bi,M & cr,M & ci)101 void matlab_real_cplx(const M& a, const M& br, const M& bi, M& cr, M& ci)
102 {
103 cr.noalias() += a * br;
104 ci.noalias() += a * bi;
105 }
106
matlab_cplx_real(const M & ar,const M & ai,const M & b,M & cr,M & ci)107 void matlab_cplx_real(const M& ar, const M& ai, const M& b, M& cr, M& ci)
108 {
109 cr.noalias() += ar * b;
110 ci.noalias() += ai * b;
111 }
112
113 template<typename A, typename B, typename C>
gemm(const A & a,const B & b,C & c)114 EIGEN_DONT_INLINE void gemm(const A& a, const B& b, C& c)
115 {
116 c.noalias() += a * b;
117 }
118
main(int argc,char ** argv)119 int main(int argc, char ** argv)
120 {
121 std::ptrdiff_t l1 = internal::queryL1CacheSize();
122 std::ptrdiff_t l2 = internal::queryTopLevelCacheSize();
123 std::cout << "L1 cache size = " << (l1>0 ? l1/1024 : -1) << " KB\n";
124 std::cout << "L2/L3 cache size = " << (l2>0 ? l2/1024 : -1) << " KB\n";
125 typedef internal::gebp_traits<Scalar,Scalar> Traits;
126 std::cout << "Register blocking = " << Traits::mr << " x " << Traits::nr << "\n";
127
128 int rep = 1; // number of repetitions per try
129 int tries = 2; // number of tries, we keep the best
130
131 int s = 2048;
132 int cache_size = -1;
133
134 bool need_help = false;
135 for (int i=1; i<argc; ++i)
136 {
137 if(argv[i][0]=='s')
138 s = atoi(argv[i]+1);
139 else if(argv[i][0]=='c')
140 cache_size = atoi(argv[i]+1);
141 else if(argv[i][0]=='t')
142 tries = atoi(argv[i]+1);
143 else if(argv[i][0]=='p')
144 rep = atoi(argv[i]+1);
145 else
146 need_help = true;
147 }
148
149 if(need_help)
150 {
151 std::cout << argv[0] << " s<matrix size> c<cache size> t<nb tries> p<nb repeats>\n";
152 return 1;
153 }
154
155 if(cache_size>0)
156 setCpuCacheSizes(cache_size,96*cache_size);
157
158 int m = s;
159 int n = s;
160 int p = s;
161 A a(m,p); a.setRandom();
162 B b(p,n); b.setRandom();
163 C c(m,n); c.setOnes();
164 C rc = c;
165
166 std::cout << "Matrix sizes = " << m << "x" << p << " * " << p << "x" << n << "\n";
167 std::ptrdiff_t mc(m), nc(n), kc(p);
168 internal::computeProductBlockingSizes<Scalar,Scalar>(kc, mc, nc);
169 std::cout << "blocking size (mc x kc) = " << mc << " x " << kc << "\n";
170
171 C r = c;
172
173 // check the parallel product is correct
174 #if defined EIGEN_HAS_OPENMP
175 int procs = omp_get_max_threads();
176 if(procs>1)
177 {
178 #ifdef HAVE_BLAS
179 blas_gemm(a,b,r);
180 #else
181 omp_set_num_threads(1);
182 r.noalias() += a * b;
183 omp_set_num_threads(procs);
184 #endif
185 c.noalias() += a * b;
186 if(!r.isApprox(c)) std::cerr << "Warning, your parallel product is crap!\n\n";
187 }
188 #elif defined HAVE_BLAS
189 blas_gemm(a,b,r);
190 c.noalias() += a * b;
191 if(!r.isApprox(c)) std::cerr << "Warning, your product is crap!\n\n";
192 #else
193 gemm(a,b,c);
194 r.noalias() += a.cast<Scalar>() * b.cast<Scalar>();
195 if(!r.isApprox(c)) std::cerr << "Warning, your product is crap!\n\n";
196 #endif
197
198 #ifdef HAVE_BLAS
199 BenchTimer tblas;
200 c = rc;
201 BENCH(tblas, tries, rep, blas_gemm(a,b,c));
202 std::cout << "blas cpu " << tblas.best(CPU_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/tblas.best(CPU_TIMER))*1e-9 << " GFLOPS \t(" << tblas.total(CPU_TIMER) << "s)\n";
203 std::cout << "blas real " << tblas.best(REAL_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/tblas.best(REAL_TIMER))*1e-9 << " GFLOPS \t(" << tblas.total(REAL_TIMER) << "s)\n";
204 #endif
205
206 BenchTimer tmt;
207 c = rc;
208 BENCH(tmt, tries, rep, gemm(a,b,c));
209 std::cout << "eigen cpu " << tmt.best(CPU_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/tmt.best(CPU_TIMER))*1e-9 << " GFLOPS \t(" << tmt.total(CPU_TIMER) << "s)\n";
210 std::cout << "eigen real " << tmt.best(REAL_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/tmt.best(REAL_TIMER))*1e-9 << " GFLOPS \t(" << tmt.total(REAL_TIMER) << "s)\n";
211
212 #ifdef EIGEN_HAS_OPENMP
213 if(procs>1)
214 {
215 BenchTimer tmono;
216 omp_set_num_threads(1);
217 Eigen::internal::setNbThreads(1);
218 c = rc;
219 BENCH(tmono, tries, rep, gemm(a,b,c));
220 std::cout << "eigen mono cpu " << tmono.best(CPU_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/tmono.best(CPU_TIMER))*1e-9 << " GFLOPS \t(" << tmono.total(CPU_TIMER) << "s)\n";
221 std::cout << "eigen mono real " << tmono.best(REAL_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/tmono.best(REAL_TIMER))*1e-9 << " GFLOPS \t(" << tmono.total(REAL_TIMER) << "s)\n";
222 std::cout << "mt speed up x" << tmono.best(CPU_TIMER) / tmt.best(REAL_TIMER) << " => " << (100.0*tmono.best(CPU_TIMER) / tmt.best(REAL_TIMER))/procs << "%\n";
223 }
224 #endif
225
226 #ifdef DECOUPLED
227 if((NumTraits<A::Scalar>::IsComplex) && (NumTraits<B::Scalar>::IsComplex))
228 {
229 M ar(m,p); ar.setRandom();
230 M ai(m,p); ai.setRandom();
231 M br(p,n); br.setRandom();
232 M bi(p,n); bi.setRandom();
233 M cr(m,n); cr.setRandom();
234 M ci(m,n); ci.setRandom();
235
236 BenchTimer t;
237 BENCH(t, tries, rep, matlab_cplx_cplx(ar,ai,br,bi,cr,ci));
238 std::cout << "\"matlab\" cpu " << t.best(CPU_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/t.best(CPU_TIMER))*1e-9 << " GFLOPS \t(" << t.total(CPU_TIMER) << "s)\n";
239 std::cout << "\"matlab\" real " << t.best(REAL_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/t.best(REAL_TIMER))*1e-9 << " GFLOPS \t(" << t.total(REAL_TIMER) << "s)\n";
240 }
241 if((!NumTraits<A::Scalar>::IsComplex) && (NumTraits<B::Scalar>::IsComplex))
242 {
243 M a(m,p); a.setRandom();
244 M br(p,n); br.setRandom();
245 M bi(p,n); bi.setRandom();
246 M cr(m,n); cr.setRandom();
247 M ci(m,n); ci.setRandom();
248
249 BenchTimer t;
250 BENCH(t, tries, rep, matlab_real_cplx(a,br,bi,cr,ci));
251 std::cout << "\"matlab\" cpu " << t.best(CPU_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/t.best(CPU_TIMER))*1e-9 << " GFLOPS \t(" << t.total(CPU_TIMER) << "s)\n";
252 std::cout << "\"matlab\" real " << t.best(REAL_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/t.best(REAL_TIMER))*1e-9 << " GFLOPS \t(" << t.total(REAL_TIMER) << "s)\n";
253 }
254 if((NumTraits<A::Scalar>::IsComplex) && (!NumTraits<B::Scalar>::IsComplex))
255 {
256 M ar(m,p); ar.setRandom();
257 M ai(m,p); ai.setRandom();
258 M b(p,n); b.setRandom();
259 M cr(m,n); cr.setRandom();
260 M ci(m,n); ci.setRandom();
261
262 BenchTimer t;
263 BENCH(t, tries, rep, matlab_cplx_real(ar,ai,b,cr,ci));
264 std::cout << "\"matlab\" cpu " << t.best(CPU_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/t.best(CPU_TIMER))*1e-9 << " GFLOPS \t(" << t.total(CPU_TIMER) << "s)\n";
265 std::cout << "\"matlab\" real " << t.best(REAL_TIMER)/rep << "s \t" << (double(m)*n*p*rep*2/t.best(REAL_TIMER))*1e-9 << " GFLOPS \t(" << t.total(REAL_TIMER) << "s)\n";
266 }
267 #endif
268
269 return 0;
270 }
271
272