Я не собираюсь нести ответственность за этот код (или платить алименты на его содержание), поэтому: нет, код не мой. Это всё, что я могу ответить на заданный вопрос.
По хорошему, тогда нужно на asm писать под конкретный процессор.
Для RISC-V я бы может и заюзал «мнемоники машинного кода», писать на asm для современного AMD64-совместимого CPU, который по факту CISC-on-VLIW и работает по принципу «фиг пойми», мне неохота.
Что касается NVIDIA, то я почти на 100% уверен, что Thrust «под капотом» написан на asm'е инженерами с Ph.D и 7-значной зарплатой: Huang'у нужны «красивые» цифры в бенчмарках, чтобы продавать.
Денис, описание бенчмарка на гитхаб не соответствует действительности: вместо 'GPU vs CPU' должно быть 'Thrust vs STL'. На задаче сортировки это сравнение выглядит особенно «нелепо»: однопоточная «логарифмическая» сортировка сравнением, запущенная на 1-ом ядре с TDP 10Вт и ценой ~$100 против «линейной» поразрядной сортировки подсчётом, утилизирующей «видяху» с TDP 250Вт и ценой $5000.
#ifndef ___MSORT_H
#define ___MSORT_H
#include <pthread.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define MCMP(func_name) _Bool func_name(const void *a, const void *b)
typedef _Bool (* ___cmp_t) (const void *, const void *);
#define MSORT_DEC(func_name) \
void* func_name(void *Arr, size_t Arr_Len, void *Ext_Buff, ___cmp_t Is_Unordered);
#define MSORT_MT_DEC(func_name) \
void* func_name(size_t T, void *Arr, size_t Arr_Len, void *Ext_Buff);
#define ___I *x++ = *i++
#define ___J *x++ = *j++
#define MAIN_VOID int main(void)
#define MAIN_ARGS int main(const int argc, const char *const *const argv)
#define TIME_PAIR struct timespec ___t1, ___t2
#define TIME_GET_START clock_gettime(CLOCK_REALTIME, &___t1)
#define TIME_GET_STOP clock_gettime(CLOCK_REALTIME, &___t2)
#define TIME_DIFF_NS ((___t2.tv_nsec + ___t2.tv_sec * (int64_t)1e9) - \
(___t1.tv_nsec + ___t1.tv_sec * (int64_t)1e9))
#define TIME_DIFF_MS (TIME_DIFF_NS / (int64_t)1e6)
#define TIME_DIFF_MCS (TIME_DIFF_NS / (int64_t)1e3)
#define TIME_DIFF_EXEC(EXEC, DIFF_SUFFIX, DIFF_LVAL) \
TIME_GET_START; EXEC; TIME_GET_STOP; DIFF_LVAL = TIME_DIFF_##DIFF_SUFFIX
#define MSORT(func_name, TYPE) /* TYPE: size_t for pointers */ \
void* func_name(void *const A, const size_t L, void *const E, ___cmp_t cmp) { \
if (L < 2) return A; size_t z = 2, y = 1, l = L / 2, t = 0; \
TYPE tmp, *i = A, *j = i + 1, *I, *J, *a = A, *e = E, *x = a + L, *p; \
for (; l--; i += 2, j += 2) { if (cmp(i, j)) { tmp = *j; *j = *i; *i = tmp; } } goto t; \
y2: for (; l--; I += 4, J += 4, i += 2, j += 2) { \
if (cmp(i, j)) ___J; else ___I; \
if (cmp(i, j)) { ___J; if (j == J) { ___I; ___I; continue; } } \
else { ___I; if (i == I) { ___J; ___J; continue; } } \
if (cmp(i, j)) { ___J; ___I; } else { ___I; ___J; } } goto t; \
y4: for (; l--; I += 8, J += 8, i += 4, j += 4) { \
for (size_t k = 4; --k;) { if (cmp(i, j)) ___J; else ___I; } \
for (;;) { if (cmp(i, j)) { ___J; if (j == J) { for (; i < I;) ___I; break; } } \
else { ___I; if (i == I) { for (; j < J;) ___J; break; } } } } goto t; \
y_: for (; l--; I += z, J += z, i += y, j += y) { \
for (size_t a, b, k = y;;) { for (; --k;) { if (cmp(i, j)) ___J; else ___I; } \
a = I - i, b = J - j, k = a > b ? b : a; if (k < 4) break; } \
for (;;) { if (cmp(i, j)) { ___J; if (j == J) { for (; i < I;) ___I; break; } } \
else { ___I; if (i == I) { for (; j < J;) ___J; break; } } } } \
t: if (t == 0) { if (L % z == 0) goto n; else { \
if (z == 2) { i = (I = j = (J = a + L) - 1) - 2; x = e + (L - 3); t = 3; } \
else { i = (I = x) - z; x = (j = (J = a + L) - y) - z; t = z + y; } } } \
else { J = (j = p) + t; if (L % z < t) { t += y; } \
else { i = (I = x) - z; x = a + (i - e); t += z; } } p = x; \
for (size_t a, b, k = y;;) { for (; --k;) { if (cmp(i, j)) ___J; else ___I; } \
a = I - i, b = J - j, k = a > b ? b : a; if (k < 4) break; } \
for (;;) { if (cmp(i, j)) { ___J; if (j == J) { for (; i < I;) ___I; break; } } \
else { ___I; \
if (i == I) { if (x != j) { for (; j < J;) ___J; } else x = J; break; } } } \
n: if (L / z == 1) return (x == a + L) ? a : e; if (z > 2) a = e, e = (e == E) ? A : E; \
z = 2 * (y *= 2), l = L / z, J = (j = (I = (i = a) + y)) + y, x = e; if (L % z < t) l--; \
if (y > 4) goto y_; else if (y == 2) goto y2; else goto y4; }
typedef struct ___mttfd_t { void *A, *E, *V; size_t L, X; } ___mttfd_t;
#define MSORT_MT(before, func_name, TYPE, msort_func_name, cmp) /* TYPE: size_t for pointers */ \
static void* func_name##_T1(void *data) { ___mttfd_t *d = data; \
d->A = msort_func_name(d->A, d->L, d->E, cmp); return NULL; } \
static void* func_name##_T2(void *data) { ___mttfd_t *d = data; \
TYPE *i = d->A, *I = i + d->L, \
*j = (d->V == NULL) ? I : d->V, *J = j + (d->L + d->X), *x = d->E; \
for (size_t a, b, k = d->L;;) { for (; --k;) { if (cmp(i, j)) ___J; else ___I; } \
a = I - i, b = J - j, k = a > b ? b : a; if (k < 4) break; } \
for (;;) { if (cmp(i, j)) { ___J; if (j == J) { for (; i < I;) ___I; break; } } \
else { ___I; if (i == I) { if (x != j) { for (; j < J;) ___J; } break; } } } \
return NULL; } \
before void* func_name(size_t T, void *const A, const size_t L, void *const E) { \
if (T < 2 || T > 256 || (T & (T - 1))) return NULL; \
if (L < 4096) return msort_func_name(A, L, E, cmp); \
for (; T > 2; T /= 2) if (L / T > 2047) break; ___mttfd_t d[T]; _Bool R = 0; pthread_t th[T]; \
size_t t, P = T / 2, y = L / T, z = y * 2, X = L % T; TYPE *a = A, *e = E, *i = a, *j = e, *V; \
for (t = 0; t < T; t++) { \
d[t].A = i; i += y; d[t].E = j; j += y; d[t].L = (t < T - 1) ? y : y + X; \
if (pthread_create(th + t, NULL, func_name##_T1, d + t) != 0) { R = 1; T = t; break; } } \
for (t = 0; t < T; t++) if (pthread_join(th[t], NULL) != 0) R = 1; if (R) return NULL; \
if (d[0].A == d[0].E) { a = E, e = A; V = (d[T - 1].A != d[T - 1].E) ? d[T - 1].A : NULL; } \
else V = (d[T - 1].A == d[T - 1].E) ? d[T - 1].A : NULL; \
for(i = a, j = e, T /= 2; T; T /= 2, y *= 2, z *= 2, a = e, e = (e == E) ? A : E, i = a, j = e) { \
for (t = 0; t < T; t++) { \
d[t].A = i; i += z; d[t].E = j; j += z; d[t].L = y; d[t].X = (t < T - 1) ? 0 : X; \
d[t].V = (T == P && t == T - 1) ? V : NULL; \
if (pthread_create(th + t, NULL, func_name##_T2, d + t) != 0) { R = 1; T = t; break; } } \
for (t = 0; t < T; t++) if (pthread_join(th[t], NULL) != 0) R = 1; if (R) return NULL; } \
return a; }
#endif
m.c:
#include <stdio.h>
#include "msort.h"
#define N (size_t)1e8
#define TYPE int32_t
static inline MCMP(cmp) { return *(const TYPE *)a > *(const TYPE *)b; }
static inline MSORT(sort, TYPE)
MSORT_MT(static inline, sort_mt, TYPE, sort, cmp)
TYPE R[N], A[N], B[N], *v, *e; size_t ms, t = 2, T, i;
MAIN_ARGS {
if (argc != 2 || sscanf(argv[1], "%lu" , &T) == 0
|| T < 2 || T > 256 || (T & (T - 1))) goto ae;
srand(time(NULL)); for (; i < N; i++) R[i] = rand(); TIME_PAIR;
for (; t <= T; t *= 2) {
printf("%31ld threads:", t); memcpy(A, R, sizeof(TYPE) * N);
TIME_DIFF_EXEC(e = sort_mt(t, A, N, B), MS, ms); if (e == NULL) goto te;
for (v = e; v < e + (N - 1); v++) if (cmp(v, v + 1)) goto ve;
printf("%40ld\n", ms);
}
return 0;
ae: printf("Usage: './m <max_thread_count>', where\n"
" <max_thread_count> - power of 2 in inclusive range 2..256\n");
return 1;
te: printf(" ...THREAD_ERROR\n");
return 1;
ve: printf(" ...VERIFICATION_ERROR, #%ld\n", v - e);
return 1;
}
Big_O:
1e8 elem:
- MSORT: 27.00 N
- MSORT_MT:
- 2 threads: 14.00 N
- 4 threads: 7.75 N
- 8 threads: 4.75 N
- 16 threads: 3.30 N
- 32 threads: 2.60 N
- 64 threads: 2.30 N
- 128 threads: 2.15 N
- 256 threads: 2.05 N
def big_O log_N
log_N = Float log_N
for i in 1..8
return if i > log_N
x = (log_N - i) / 2 ** i
a = 1.0; i.times { x += a; a /= 2 }
puts "\t- %3d: %13.10f" % [2 ** i, x]
end
end
https://github.com/CppCon/CppCon2019/blob/master/Presentations/speed_is_found_in_the_minds_of_people/speed_is_found_in_the_minds_of_people__andrei_alexandrescu__cppcon_2019.pdf
http://www.win-vector.com/blog/2008/04/sorting-in-anger/
- for (i = 0; i < T; i++) *R += (TYPE)(r[i]); return 0; \
+ for (i = 0; i < T; i++) *R += *((TYPE *)(r + i)); return 0; \
Файл назван 'dram.c', потому что на моей машине всё упирается в 40ГБ/с пропускной способности памяти: наращивание числа птоков для 8-байтовых типов данных теряет смысл при 2-х тредах для 'transform_*' и при 4-х для 'sum_*'. На «железе» с числом каналов памяти большем, чем у 1700X производительность должна быть выше.
#include <pthread.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#define MAIN_VOID int main(void)
#define MAIN_ARGS int main(const int argc, const char *const *const argv)
#define TIME_PAIR struct timespec ___t1, ___t2
#define TIME_GET_START clock_gettime(CLOCK_REALTIME, &___t1)
#define TIME_GET_STOP clock_gettime(CLOCK_REALTIME, &___t2)
#define TIME_DIFF_NS ((___t2.tv_nsec + ___t2.tv_sec * (int64_t)1e9) - \
(___t1.tv_nsec + ___t1.tv_sec * (int64_t)1e9))
#define TIME_DIFF_MS (TIME_DIFF_NS / (int64_t)1e6)
#define TIME_DIFF_MCS (TIME_DIFF_NS / (int64_t)1e3)
#define TIME_DIFF_EXEC(EXEC, DIFF_SUFFIX, DIFF_LVAL) \
TIME_GET_START; EXEC; TIME_GET_STOP; DIFF_LVAL = TIME_DIFF_##DIFF_SUFFIX
enum { N = 100800000, T = 4 };
#define SUM(func_name, TYPE) \
void* func_name##_T(void *_) { \
TYPE *p = *(TYPE **)_; TYPE r = 0; \
for (int i = 0; i < N / T; i++) r += p[i]; \
*(TYPE *)_ = r; \
return NULL; \
} \
_Bool func_name(TYPE *p, TYPE *R) { \
pthread_t th[T]; size_t r[T]; int t, i, L = N / T; _Bool e = 0; \
for (t = 0; t < T; t++) { \
r[t] = (size_t)(p + L * t); \
if (pthread_create(th + t, NULL, func_name##_T, r + t) != 0) { e = 1; break; } \
} \
for (i = 0; i < t; i++) if (pthread_join(th[i], NULL) != 0) e = 1; \
if (e) return e; *R = 0; \
for (i = 0; i < T; i++) *R += (TYPE)(r[i]); return 0; \
}
#define TRANSFORM(func_name, TYPE) \
void* func_name##_T(void *_) { \
TYPE *p = _; \
for (int i = 0; i < N / T; i++) p[i] = p[i] * p[i] + p[i] / 3; \
return NULL; \
} \
_Bool func_name(TYPE *p) { \
pthread_t th[T]; int t, i, L = N / T; _Bool e = 0; \
for (t = 0; t < T; t++) { \
if (pthread_create(th + t, NULL, func_name##_T, p + L * t) != 0) { e = 1; break; } \
} \
for (i = 0; i < t; i++) if (pthread_join(th[i], NULL) != 0) e = 1; \
return e; \
}
SUM(sum_d, double)
SUM(sum_l, int64_t)
SUM(sum_i, int32_t)
TRANSFORM(transform_d, double)
TRANSFORM(transform_l, int64_t)
TRANSFORM(transform_i, int32_t)
MAIN_VOID {
if (N % T) { printf("line 64\n"); return 1; } _Bool e;
double *d = malloc(N * 8), rd; if (d == NULL) { printf("line 66\n"); return 1; }
int64_t *l = (int64_t *)d, rl;
int32_t *i = (int32_t *)d, ri;
srand(time(NULL)); TIME_PAIR; size_t sd, td, sl, tl, si, ti;
for (int i = 0; i < N; i++) d[i] = rand();
TIME_DIFF_EXEC(e = sum_d(d, &rd), MS, sd); if (e) { printf("line 72\n"); return 1; }
printf("sum d = %30.1lf, %7zu ms\n", rd, sd);
for (int i = 0; i < N; i++) d[i] = rand();
TIME_DIFF_EXEC(e = transform_d(d), MS, td); if (e) { printf("line 75\n"); return 1; }
printf("d[rand] = %30.1lf, %7zu ms\n", d[rand() % N], td);
for (int i = 0; i < N; i++) l[i] = rand();
TIME_DIFF_EXEC(e = sum_l(l, &rl), MS, sl); if (e) { printf("line 79\n"); return 1; }
printf("sum l = %30ld, %7zu ms\n", rl, sl);
for (int i = 0; i < N; i++) l[i] = rand();
TIME_DIFF_EXEC(e = transform_l(l), MS, tl); if (e) { printf("line 82\n"); return 1; }
printf("l[rand] = %30ld, %7zu ms\n", l[rand() % N], tl);
for (int j = 0; j < N; j++) i[j] = rand();
TIME_DIFF_EXEC(e = sum_i(i, &ri), MS, si); if (e) { printf("line 86\n"); return 1; }
printf("sum i = %30d, %7zu ms\n", ri, si);
for (int j = 0; j < N; j++) i[j] = rand();
TIME_DIFF_EXEC(e = transform_i(i), MS, ti); if (e) { printf("line 89\n"); return 1; }
printf("i[rand] = %30d, %7zu ms\n", i[rand() % N], ti);
return 0;
}
$ gcc -flto -O3 dram.c -lpthread --static && ./a.out
sum d = 19410551262083469312.0, 28 ms
d[rand] = 2111606310059096576.0, 51 ms
sum l = 108243668190261742, 19 ms
l[rand] = 4035494864057467728, 53 ms
sum i = -529253480, 9 ms
i[rand] = 841548210, 29 ms
$ cpuid | grep 'brand =' | uniq
brand = "AMD Ryzen 7 1700X Eight-Core Processor "
dram_again.c:
Что касается NVIDIA, то я почти на 100% уверен, что Thrust «под капотом» написан на asm'е инженерами с Ph.D и 7-значной зарплатой: Huang'у нужны «красивые» цифры в бенчмарках, чтобы продавать.
m.c:
Big_O:
Файл назван 'dram.c', потому что на моей машине всё упирается в 40ГБ/с пропускной способности памяти: наращивание числа птоков для 8-байтовых типов данных теряет смысл при 2-х тредах для 'transform_*' и при 4-х для 'sum_*'. На «железе» с числом каналов памяти большем, чем у 1700X производительность должна быть выше.