小さい行列積をたくさん計算するサンプルを追加

[cuda.git] / mult_matrices / mult_matrices.cpp

diff --git a/mult_matrices/mult_matrices.cpp b/mult_matrices/mult_matrices.cpp
new file mode 100644 (file)
index 0000000..bddb96e
--- /dev/null
+++ b/mult_matrices/mult_matrices.cpp
@@ -0,0 +1,164 @@
+/*
+       Copyright (C) 2012, 2013  fmaj7b5.info
+
+       This program is free software: you can redistribute it and/or modify
+       it under the terms of the GNU General Public License as published by
+       the Free Software Foundation, either version 2 of the License, or
+       (at your option) any later version.
+
+       This program is distributed in the hope that it will be useful,
+       but WITHOUT ANY WARRANTY; without even the implied warranty of
+       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+       GNU General Public License for more details.
+
+       You should have received a copy of the GNU General Public License
+       along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+// mult_matrices.cpp : \83R\83\93\83\\81[\83\8b \83A\83v\83\8a\83P\81[\83V\83\87\83\93\82Ì\83G\83\93\83g\83\8a \83|\83C\83\93\83g\82ð\92è\8b`\82µ\82Ü\82·\81B
+//
+
+#include "stdafx.h"
+
+#include "mult_matrices.h"
+
+using namespace FM7b5;
+
+const size_t g_num(1024 * 1024);
+const size_t g_num_run(10);
+
+std::default_random_engine g_re;
+
+static void random_matrices(float* mat44, const size_t num);
+
+static void disp_matrix(const float* m44);
+static inline void mult_matrix(float* C, const float* A, const float* B);
+static void mult_matrices(float* C, const float* A, const float* B, const size_t num);
+
+static double flops(const size_t num, const double elapsed_ms);
+
+#define FM7b5_USE_CPU
+//#define FM7b5_TRANSPOSED
+
+int _tmain(int argc, _TCHAR* argv[])
+{
+       /* Column-major 4x4 matrices (i.e. M(i, j) = M[i + 4*j] */
+       std::vector<float> A(g_num * 4*4), B(g_num * 4*4), C(g_num * 4*4), C_gpu(g_num * 4*4);
+
+       random_matrices(&A[0], g_num);
+       random_matrices(&B[0], g_num);
+       
+       ULONGLONG start, finish;
+       double elapsed_ms;
+
+#ifdef FM7b5_USE_CPU
+       /* CPU */
+       start = GetTickCount64();
+       for (size_t nrun = 0; nrun < g_num_run; ++nrun) {
+               mult_matrices(&C[0], &A[0], &B[0], g_num);
+       }
+       finish = GetTickCount64();
+       elapsed_ms = static_cast<double>(finish - start);
+       printf("1 CPU: %.1f [ms] (%.1f GFLOPS)\n\n", elapsed_ms, static_cast<double>(g_num_run) * flops(g_num, elapsed_ms) / 1e9);
+#endif
+
+       /* GPU */
+       mult_matrices_init_gpu();
+       start = GetTickCount64();
+       for (size_t nrun = 0; nrun < g_num_run; ++nrun) {
+               mult_matrices_gpu(&C_gpu[0], &A[0], &B[0], g_num);
+       }
+       finish = GetTickCount64();
+       elapsed_ms = static_cast<double>(finish - start);
+       printf("GPU: %.1f [ms] (%.1f GFLOPS)\n\n", elapsed_ms, static_cast<double>(g_num_run) * flops(g_num, elapsed_ms) / 1e9);
+
+#if 0
+       disp_matrix(&A[g_num-1]);
+       disp_matrix(&B[g_num-1]);
+
+# ifdef FM7b5_USE_CPU
+       disp_matrix(&C[g_num-1]);
+# endif
+       disp_matrix(&C_gpu[g_num-1]);
+#endif
+
+       return 0;
+}
+
+void
+random_matrices(float* mat44, const size_t num)
+{
+       if (mat44 == NULL || num < 1) {
+               return;
+       }
+
+       std::uniform_real_distribution<float> rand_dist;
+
+       for (size_t i = 0; i < num; ++i) {
+               float* p(mat44 + 16*i);
+
+               for (size_t j = 0; j < 16; ++j)
+               {
+                       p[j] = rand_dist(g_re);
+               }
+       }
+}
+
+void
+disp_matrix(const float* m44)
+{
+       for (size_t r = 0; r < 4; ++r) {
+               for (size_t c = 0; c < 4; ++c) {
+                       printf("% 7.3f ", m44[r + 4*c]);
+               }
+               printf("\n");
+       }
+       printf("\n");
+}
+
+inline void
+mult_matrix(float* __restrict C, const float* __restrict A, const float* __restrict B)
+{
+       for (size_t i = 0; i < 16; ++i) {
+               C[i] = 0.0;
+       }
+
+#ifdef FM7b5_TRANSPOSED
+       float Bt[16];
+       for (size_t c = 0; c < 4; ++c) {
+               for (size_t r = 0; r < 4; ++r) {
+                       Bt[r + 4*c] = B[c + 4*r];
+               }
+       }
+#endif
+
+       for (size_t k = 0; k < 4; ++k) {
+               for (size_t c = 0; c < 4; ++c) {
+                       for (size_t r = 0; r < 4; ++r) {
+#ifndef FM7b5_TRANSPOSED
+                               C[r + 4*c] += A[r + 4*k] * B[k + 4*c];
+#else
+                               C[r + 4*c] += A[r + 4*k] * Bt[c + 4*k];
+#endif
+                       }
+               }
+       }
+}
+
+void
+mult_matrices(float* C, const float* A, const float* B, const size_t num)
+{
+#pragma omp parallel for
+       for (int i = 0; i < static_cast<int>(num); ++i) {
+               mult_matrix(C + 16*i, A + 16*i, B + 16*i);
+       }
+}
+
+double
+flops(const size_t num, const double elapsed_ms)
+{
+       /* num of multiplications and additions in a single matrix-matrix multiplication */
+       const int ops(4 * 4 * (4 + 3));
+
+       return static_cast<double>(num * ops) * 1.0e3 / elapsed_ms;
+}