Последно редактирано на 27.09.2020 от |, видяно: 1985 пъти.
Не разбрах защо сравнявате нещо дето ще се чете от харда с нещо дето ще се чете от паметта?
И аз нямам никаква идея, но специалистите по програмиране твърдят, че производителността на релационните бази данни, и особено sqlite е значително по-голяма от специално написан код решаващ определения проблем.
От друга страна, с новите подобрения кода на GPU сравнява 200К стринга от колекция Б със 100К стринга от колекция А за 1335 секунди, или 0.0068 секунди на стринг.
Или, ако използваме мерките и теглилки препоръчани от гениалния програмист на Асемблер, един стринг от Б се обработва за .00000001103670634920 седмици.
Последно редактирано на 27.09.2020 от Golden Gega, видяно: 1959 пъти.
Не разбрах защо сравнявате нещо дето ще се чете от харда с нещо дето ще се чете от паметта?
И аз нямам никаква идея, но специалистите по програмиране твърдят, че производителността на релационните бази данни, и особено sqlite е значително по-голяма от специално написан код решаващ определения проблем.
От друга страна, с новите подобрения кода на GPU сравнява 200К стринга от колекция Б със 100К стринга от колекция А за 1335 секунди, или 0.0068 секунди на стринг.
Или, ако използваме мерките и теглилки препоръчани от гениалния програмист на Асемблер, един стринг от Б се обработва за .00000001103670634920 седмици.
Хайде, хайде, има и програмисти дето просто им е интересно.
Аз подходих леко отпускарски към проблема, в SQL Server Management Studio една голяма част от времето се яде от извеждането на резултатите, затова се позадълбах в проблема и отчетох само чистата сметка с левенщайна, без подреждането по резултати (елиминирах показването на резултати и пълненето на таблиците с данни). Това го направих от чиста прагматичност, защото може да се окаже полезно (евентуално) в един бъдещ проект. Резултата го давам за сравнение, явно е бавно в сравнение с чисти сметки с език от по-ниско ниво.
И така, два варианта:
1) Една таблица с 1 000 000 реда (същата като резултатния сет от пример 2), старт/стоп/време в секунди
За по-големи сметки честно казано не ми се чака, а и резултата явно става само за информация.
P.S. Ако някой много държи мога да постна кода с теста, стъпките с настройки и импорта на външната функция (SQL Server много държи импортираните функции да са подписани), както и опциите по Management Studio за провеждане на теста. Ако бях гъзар щях да ползвам конзола ама нали съм мързеливо архитектче трябва да пазя имидж.
Реших да го нахвърлям набързо да видим за колко ще се изпълни върху GPU-то (1080ti). 100 хиляди сравнения, низовете са случайни с дължина до 200 байта, kernel-а е тотално неоптимизиран, кажи-речи копи-пейст на примерна имплементация. Обаче естествено има измама, в паметта, низовете са оформени като "паскалски" такива, с първите 4 байта дължината на низа и нататък, това спестява strlen-а имплементиран в kernel-а, което не че е сложно, ама като нагази в глобалната памет да го смята, ще стане лееееееекинко по-тегаво.
main.c:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <CL/cl.h>
#define MAX_SOURCE_SIZE (0x100000)
#define SET_SIZE 100032
#define STRING_MAX 200
#define MY_STRING "boiko borisov pederast pederast! boiko borisov pederast pederast!"
int main(void)
{
// Create the inputs
int i,j,sz;
char *A = (char*)malloc(SET_SIZE*(STRING_MAX + sizeof(int)));
char *B = (char*)malloc(STRING_MAX + sizeof(int));
char *results = (char*)malloc(SET_SIZE * sizeof(int));
char temp[STRING_MAX + sizeof(int)];
double exec_time = 0.0;
memset(A,SET_SIZE*(STRING_MAX + sizeof(int)), 0);
memset(B,STRING_MAX + sizeof(int), 0);
// Create the set A, SET_SIZE pascal-type strings up to STRING_MAX length
for(i = 0; i < SET_SIZE; i++)
{
sz = 1 + rand() % (STRING_MAX-1);
for (j = 0; j < sz; j++)
{
A[i * (STRING_MAX + sizeof(int)) + sizeof(int) + j] = (char)(rand() % 57 + 65);
}
memcpy((char*)(A + i * (STRING_MAX + sizeof(int))), &sz, sizeof(int));
}
// Create the set B, single pascal-type string
memcpy((char*)(B + sizeof(int)), MY_STRING, strlen(MY_STRING));
sz = strlen(MY_STRING);
memcpy(B, &sz, 4);
// Load the kernel
FILE *fp;
char *source_str;
size_t source_size;
fp = fopen("kernel.cl", "r");
if (!fp) {
fprintf(stderr, "Failed to load kernel.\n");
exit(1);
}
source_str = (char*)malloc(MAX_SOURCE_SIZE);
source_size = fread( source_str, 1, MAX_SOURCE_SIZE, fp);
fclose( fp );
// OpenCL stuff initialization
cl_platform_id platform_id = NULL;
cl_device_id device_id = NULL;
cl_uint ret_num_devices;
cl_uint ret_num_platforms;
cl_int ret = clGetPlatformIDs(1, &platform_id, &ret_num_platforms);
ret = clGetDeviceIDs( platform_id, CL_DEVICE_TYPE_DEFAULT, 1,
&device_id, &ret_num_devices);
cl_context context = clCreateContext( NULL, 1, &device_id, NULL, NULL, &ret);
cl_command_queue command_queue = clCreateCommandQueue(context, device_id, 0, &ret);
// Create the buffers and copy to device
cl_mem set_a = clCreateBuffer(context, CL_MEM_READ_ONLY, SET_SIZE*(STRING_MAX + sizeof(int)), NULL, &ret);
cl_mem set_b = clCreateBuffer(context, CL_MEM_READ_ONLY, STRING_MAX + sizeof(int), NULL, &ret);
cl_mem results_set = clCreateBuffer(context, CL_MEM_WRITE_ONLY, SET_SIZE * sizeof(int), NULL, &ret);
// Copy the lists A and B to their respective memory buffers
ret = clEnqueueWriteBuffer(command_queue, set_a, CL_TRUE, 0, SET_SIZE*(STRING_MAX + sizeof(int)), A, 0, NULL, NULL);
ret = clEnqueueWriteBuffer(command_queue, set_b, CL_TRUE, 0, STRING_MAX + sizeof(int), B, 0, NULL, NULL);
// Compile
cl_program program = clCreateProgramWithSource(context, 1, (const char **)&source_str, (const size_t *)&source_size, &ret);
ret = clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);
if (ret != 0)
{
size_t log_size;
clGetProgramBuildInfo(program, device_id, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
char *log = (char *) malloc(log_size);
clGetProgramBuildInfo(program, device_id, CL_PROGRAM_BUILD_LOG, log_size, log, NULL);
printf("%s\n", log);
}
cl_kernel kernel = clCreateKernel(program, "levenshtein", &ret);
ret = clSetKernelArg(kernel, 0, sizeof(cl_mem), (void *)&set_a);
ret = clSetKernelArg(kernel, 1, sizeof(cl_mem), (void *)&set_b);
ret = clSetKernelArg(kernel, 2, sizeof(cl_mem), (void *)&results_set);
// Execute
size_t global_item_size = SET_SIZE;
size_t local_item_size = 64;
clock_t start = clock();
ret = clEnqueueNDRangeKernel(command_queue, kernel, 1, NULL, &global_item_size, &local_item_size, 0, NULL, NULL);
printf("execution retstatus=%d\n",ret);
clock_t end = clock();
exec_time += (double)(end - start) / CLOCKS_PER_SEC;
printf("execution time = %f seconds\n", exec_time);
unsigned int *C = (int*)malloc(sizeof(int)*SET_SIZE);
ret = clEnqueueReadBuffer(command_queue, results_set, CL_TRUE, 0, SET_SIZE * sizeof(int), C, 0, NULL, NULL);
// Print distances
//for(i = 0; i < SET_SIZE; i++)
// printf("%d = %u\n", i, C[i]);
ret = clFlush(command_queue);
ret = clFinish(command_queue);
ret = clReleaseKernel(kernel);
ret = clReleaseProgram(program);
ret = clReleaseMemObject(set_a);
ret = clReleaseMemObject(set_b);
ret = clReleaseMemObject(results_set);
ret = clReleaseCommandQueue(command_queue);
ret = clReleaseContext(context);
free(A);
free(B);
free(C);
return 0;
}
Не разбрах защо сравнявате нещо дето ще се чете от харда с нещо дето ще се чете от паметта?
И аз нямам никаква идея, но специалистите по програмиране твърдят, че производителността на релационните бази данни, и особено sqlite е значително по-голяма от специално написан код решаващ определения проблем.
От друга страна, с новите подобрения кода на GPU сравнява 200К стринга от колекция Б със 100К стринга от колекция А за 1335 секунди, или 0.0068 секунди на стринг.
Или, ако използваме мерките и теглилки препоръчани от гениалния програмист на Асемблер, един стринг от Б се обработва за .00000001103670634920 седмици.
Споменах за външните индекси де ги писахме за MSSQL - дигнати tries в паметта на сървър и реално правеха възможно бързото търсене в релационната база данни. Ta вземал съм си пуканките и чакам да ти бият решението :)
Дъртия хари говореше за монго. Като що годе уж "експерт" по монго ще кажа, че въпреки дигането на двете колекции в паметта, което става автоматик при достатъчно наличие на памет, то почти всякакви индекси ще се четат от харда. Отделно понятието join не съществува в МонгоДъйБъ. Има агрегативен оператор $lookup, който действа като left outer join, изключително бавен и в pipeline-a на агрегацията ти ще трябва да имаш нещо от типа на {$match: {joinedRecords: {$ne: []}}}, което допълнително ще закопае нещата.
А, обаче излъгах, nvidia-та са педали и изпълнението на кернела немое се мери така, това е времето за което драйвера schedule-ва kernel-а върху GPU-то. Затва реших да го измеря от изпълнение до след трансфера на резултатите обратно на хоста и сега резултатът е дооооста различен:
execution retstatus=0
execution time = 0.030778 seconds
Реално обаче това включва и host-device трансфера. Но пък то е интересно като цяло какво трябва да се мери - самият левенщайн вероятно минава за по-малко от това, сетъпването на буфери, трансферите, естествено и "подготвянето" на входните данни във вид на паскалски низове със сигурност отнема доста повече време.
А, обаче излъгах, nvidia-та са педали и изпълнението на кернела немое се мери така, това е времето за което драйвера schedule-ва kernel-а върху GPU-то. Затва реших да го измеря от изпълнение до след трансфера на резултатите обратно на хоста и сега резултатът е дооооста различен:
execution retstatus=0
execution time = 0.030778 seconds
Реално обаче това включва и host-device трансфера. Но пък то е интересно като цяло какво трябва да се мери - самият левенщайн вероятно минава за по-малко от това, сетъпването на буфери, трансферите, естествено и "подготвянето" на входните данни във вид на паскалски низове със сигурност отнема доста повече време.
Аз затова меря милиони сравнения. Ако кода работи половин час, разни подробности като сетъпване на буфери и четене на данните от диска стават съвсе маловажни.
Моят кърнъл е по-сложен, но той сравнява всички стрингове от А със всички стрингове от Б и използва shared memory:
typedef struct Oligo {
short len;
char* seq;
double qubu;
} Oligo;
typedef struct Pool {
int olnum; // number of oligos in the pool
Oligo* ols;
int olmaxlen; // maximum length of an oligo
} Pool;
__global__ void cudnaMatch(Pool *dataset, Pool *reads, Match *matches) {
extern __shared__ char shbuf[];
char *dsol, *rol;
unsigned char *rbuf;
int i, n, mdist;
int dsolen, rolen, dsolnum;
Oligo *ol, *mol;
// setup the shared pointers
dsol = shbuf;
rol = shbuf + dataset->olmaxlen + 1 + (reads->olmaxlen * 2 + 1) * threadIdx.x;
rbuf = (unsigned char *) (rol + reads->olmaxlen);
// copy the sequence the local thread works on the shared memory
n = blockIdx.x * blockDim.x + threadIdx.x;
if (n < reads->olnum) {
ol = &reads->ols[n];
rolen = ol->len;
for(i = 0; i < rolen; i++) {
rol[i] = ol->seq[i];
}
} else {
rolen = 0;
}
dsolnum = dataset->olnum;
mdist = 10000;
mol = NULL;
__syncthreads();
for(i = 0; i < dsolnum; i++) {
ol = &dataset->ols[i];
dsolen = ol->len;
// all threads copy some of the oligo
for(n = threadIdx.x; n < dsolen; n += blockDim.x) {
dsol[n] = ol->seq[n];
}
__syncthreads();
// at this point everything should be in the shared memory, calculate the distance
if (rolen > 0) {
n = ldist(dsolen, dsol, rolen, rol, rbuf);
if (n < mdist) {
mdist = n;
mol = ol;
}
}
__syncthreads();
}
if (rolen > 0) {
n = blockIdx.x * blockDim.x + threadIdx.x;
Match *m = &matches[n];
m->dist = mdist;
m->o1 = &reads->ols[blockIdx.x * blockDim.x + threadIdx.x];
m->o2 = mol;
}
}
Левещейн функцията ми е почти като твоята, но използвам два umin вместо MIN3. Забързва производителността малко, защото има по-малко бранчове.
Споменах за външните индекси де ги писахме за MSSQL - дигнати tries в паметта на сървър и реално правеха възможно бързото търсене в релационната база данни. Ta вземал съм си пуканките и чакам да ти бият решението :)
Дъртия хари говореше за монго. Като що годе уж "експерт" по монго ще кажа, че въпреки дигането на двете колекции в паметта, което става автоматик при достатъчно наличие на памет, то почти всякакви индекси ще се четат от харда. Отделно понятието join не съществува в МонгоДъйБъ. Има агрегативен оператор $lookup, който действа като left outer join, изключително бавен и в pipeline-a на агрегацията ти ще трябва да имаш нещо от типа на {$match: {joinedRecords: {$ne: []}}}, което допълнително ще закопае нещата.
Според мен всякакви бази данни са напълно безмислени в случая. Първо, данните не са чак толкова много (дори 23 милиона стринга са под 4GB) и могат да се прочетат в паметта в началото. И второ, На алгоритъма дори не му трябват всичките данни в паметта. Трябва му по-малката колекция (в случая А), а стринговете от колекцията Б може да ги стриймва. 100К стринга от колекция А са под 20MB, L3 кеша на процесора е 36MB.
Така и трябва. Иначе, ако ще викаш C библиотека, трябва да го копираш другаде с достатъчно място да сложиш нулата.
Това му е идеята, предназначен е за interop, обаче ставаше весело когато със strcpy или морален еквивалент се изпляска по-къс стринг върху BSTR променлива и после се хвърля като параметър на ком функция която чете дължината от префикса понеже по-често не крашваше :)
Споменах за външните индекси де ги писахме за MSSQL - дигнати tries в паметта на сървър и реално правеха възможно бързото търсене в релационната база данни. Ta вземал съм си пуканките и чакам да ти бият решението :)
Дъртия хари говореше за монго. Като що годе уж "експерт" по монго ще кажа, че въпреки дигането на двете колекции в паметта, което става автоматик при достатъчно наличие на памет, то почти всякакви индекси ще се четат от харда. Отделно понятието join не съществува в МонгоДъйБъ. Има агрегативен оператор $lookup, който действа като left outer join, изключително бавен и в pipeline-a на агрегацията ти ще трябва да имаш нещо от типа на {$match: {joinedRecords: {$ne: []}}}, което допълнително ще закопае нещата.
Според мен всякакви бази данни са напълно безмислени в случая. Първо, данните не са чак толкова много (дори 23 милиона стринга са под 4GB) и могат да се прочетат в паметта в началото. И второ, На алгоритъма дори не му трябват всичките данни в паметта. Трябва му по-малката колекция (в случая А), а стринговете от колекцията Б може да ги стриймва. 100К стринга от колекция А са под 20MB, L3 кеша на процесора е 36MB.
Един случай трябва да се гледа в контекста на цяла система, ако говорим за реалния свят. В академичен план например случая може да е изсипваме наведнъж 1м стрингове в едната и 1м в другата колекция и го решаваме наведнъж, реално ако тия колекции се попълват в разтегнат вариант от време в някаква база решението директно в базата може да се окаже оптимално - ако скоростта на пресмятане е съизмерима със скоростта на постъпване на данните, например.
За мен важното е че вместо за козите на Рамбо поразцъкахме малко код, а това че се и понапцувахме на тъпанар и рабиноид си е бонус.
Споменах за външните индекси де ги писахме за MSSQL - дигнати tries в паметта на сървър и реално правеха възможно бързото търсене в релационната база данни. Ta вземал съм си пуканките и чакам да ти бият решението :)
Дъртия хари говореше за монго. Като що годе уж "експерт" по монго ще кажа, че въпреки дигането на двете колекции в паметта, което става автоматик при достатъчно наличие на памет, то почти всякакви индекси ще се четат от харда. Отделно понятието join не съществува в МонгоДъйБъ. Има агрегативен оператор $lookup, който действа като left outer join, изключително бавен и в pipeline-a на агрегацията ти ще трябва да имаш нещо от типа на {$match: {joinedRecords: {$ne: []}}}, което допълнително ще закопае нещата.
Според мен всякакви бази данни са напълно безмислени в случая. Първо, данните не са чак толкова много (дори 23 милиона стринга са под 4GB) и могат да се прочетат в паметта в началото. И второ, На алгоритъма дори не му трябват всичките данни в паметта. Трябва му по-малката колекция (в случая А), а стринговете от колекцията Б може да ги стриймва. 100К стринга от колекция А са под 20MB, L3 кеша на процесора е 36MB.
Един случай трябва да се гледа в контекста на цяла система, ако говорим за реалния свят. В академичен план например случая може да е изсипваме наведнъж 1м стрингове в едната и 1м в другата колекция и го решаваме наведнъж, реално ако тия колекции се попълват в разтегнат вариант от време в някаква база решението директно в базата може да се окаже оптимално - ако скоростта на пресмятане е съизмерима със скоростта на постъпване на данните, например.
За мен важното е че вместо за козите на Рамбо поразцъкахме малко код, а това че се и понапцувахме на тъпанар и рабиноид си е бонус.
Скоростта на постъпване на данните е много висока. Един run на най-често използваните машини на Илумина 25 милиона стринга (най-голямата 1 милиард стринга). Наведнъж. После няколко месеца нищо. :)
Но този проблем който описвам тук не е толкова важен, с данни от биологично ДНК проблемите са съвсем други и биоинформатиците имат някакви решения за тях.