0 1 2 3 4 ...6 7 8 9 10 ...14 15 16 17 18 ...24 25 26 27 28 ...32 33 34 35 36
Хванахте се стринговете, като безсмислената работа е цялото дърво, тая функция може да си сложиш и в sqlite.
Ти можеш да навреш каквото си искаш, където си искаш. Но въпросът е защо? :)
Да не си като оня, дето имал чук...
Аз това имах предвид, че твоят е бавен, не стринговете, а това че зареждаш бинарно дърво на всеки старт, точно затова sqlite е по-бърз. Функцията на Левенщайн е най-неинтересното в случая, дали ще я правя на openCL или не - си е твое решение.
Абе, ти от коя Вселена пишеш? :)
Зареждането на данните в паметта е под секунда. Сравняването на 100К стринга със 200К стринга е 1300 секунди на GPU.
Хайде опитай се да си поне малко по-адекватен?
Както почнеш да работиш с данни над милион пак пиши. Дори и 100к пак не е под секунда.
И какво, като данните са над милион, сравнението между двете колекции няма ли да се увеличи също?
Казах ти, опитай се да си поне малко адекватен.
Опитай ти да си поне малко умен, всяко пускане на твоята програма е бавене и нищо друго. Данните седят база, за да не се зареждат излишно, това което правиш ти.
УжасТ, моята програма ще се бави и вместо една седмица ще отнеме седмица и 5 секунди! Край, съгласих се, пускам sqlite.
P.S. Идиот. :)
P.P.S. Я разкажи КОИ данни се зареждат излишно, да ти се посмеем. :)
Абе нещо пак не е така. Задачата е за всеки от 100-те хиляди стринга да намериш най-близкия от другата колекция. С други думи, решението ти трябва да произведе 100_000 реда на изхода.
ПП: задачата много прилича на spellchecker обаче там алгоритмите се издъниха здраво :Д Както казва мастър Линус - когато теорията и практиката се сблъскат, практиката винаги печели.
Да, само един срещу 100 хиляди, макар че произволна бройка срещу произволна бройка няма да е кой знае колко различно, няма да нарасне баш линейно времето, защото достъпите до паметта са лайняна работа, но от друга страна, няма да е и нещо много по-зле.
На AMD APP SDK-то имаше много приличен profiler, дето за нвидията не му знам еквивалента (сигурно има), та няма как да го твърдя със сигурност, но това мога да се обзаложа че е относително memory-интензивен kernel и повечето време не отива за ALU операции, а за достъпване на памет. Вероятно компилатора тук-там успява да мине тънко с cmov-подобни изпълнения, вместо да бранчва наистина, но като цяло повечето усилия за оптимизиране бих ги хвърлил за по-умно използване на паметта (и утилизация на локалната доколкото може).
Все пак обаче подозирам че каквото и както да било върху CPU-то ще върви по-бавно не заради друго, а заради паралелизма, задачата си е почти embarassingly parallel, не е чак толкова memory-интензивен kernel-а, нема достатъчно бранчване че да се обезсмисли, не вярвам да ползва особено много регистри, поради което винаги ще си намазва върху GPU-та.
Броят на резултатите трябва да е равен на броят на стринговете в колекция Б. За всеки стринг от нея трябва да намерим най-близкия стринг от А.
Иначе биоинформатиците обикновено цепят тези стрингове на подстрингове (kmers) и правят разни неща с тях. Всеки стринг се цепи на (n-k) подстринга с дължина k. Ако избереш правилното k (стрингове с дължина до 32 приятно се събират в 64-битово число), можеш да правиш разни интересни неща, de Bruijn graphs и т.н.
Но не знам дали алгоритмите им ще работят ако има прекалено много грешки.
Като цяло, специалистите по бази данни изместиха проблема в области, които са много безинтересни. Но, when you have a hammer ...
Всъщност кефещото поне за мен в темата е точно противопоставянето на database и 'true' програмиране подходите
Засега администратора е водещ за наградата идиотщина на темата с изцепката (перифразирам) "ти обработваш данни, sql е създаден за обработка на данни, следователно е логично да използваш sql". :)
Всъщност the promised land of SQL е точно идеята че ще бъдем свободни от ненужните досадни детайли и само ще изразяваме чистите си идеи ползвайки синтаксис близък до set theory а могъщият енджин отдолу ще транслира това това до макс ефективен код обаче както се вижда каубойското пуцане с пойнтери на някои нива бие с нокаут :)
Което пък е забавно понеже ти се ебаваше с Джон на тема hand tuned assembly vs compilers :)
Е няма как да сравняваш "дърво" като релационните данни, и особено идиотщина като SQL с език за програмиране, от колкото и да е високо ниво. Това наистина е все едно да използваш чук за да завиеш болт.
SQL е примордиален DSL с ниско качество като абстракция и съответно куп грозотии отгоре и особености при имплементациите ама това е положението, Еволюцията бачка така, очевидно няма време за ПРАВИЛЕН grand design щом джаваскрипт и C царуват :)
Ти мани това, ама SQL-а, понякога става по-засукан отколкото да си го напишеш на чисто С, акто си въртиш цикли и си избираш каквото ти трябва.
Те SQL сенсеите дишащи stored procedures и user-defined functions са кунгфу мастъри, с case/coalesce фатки само сигурно мое да напишат ос кърнъл :)
IMNSHO SQL може да се конкурира с Бейсик за най-дървен популярен език. Дори Фортран (игнорирайки дали е популярен или не) е по-добре от него.
Реших да изпълня пълната задача, да се намери стринга с най-близка дистанция бла-бла, 1024 срещу 100 хиляди.
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 SET2_SIZE (1024)
#define STRING_MAX 200
int main(void)
{
// Create the inputs
int i,j,sz;
char *A = (char*)malloc(SET_SIZE*(STRING_MAX + sizeof(int)));
char *B = (char*)malloc(SET2_SIZE*(STRING_MAX + sizeof(int)));
char *results = (char*)malloc(SET_SIZE * sizeof(int));
char temp[STRING_MAX + sizeof(int)];
double exec_time = 0.0;
srand(1234);
memset(A,SET_SIZE*(STRING_MAX + sizeof(int)), 0);
memset(B,SET2_SIZE*(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
for(i = 0; i < SET2_SIZE; i++)
{
sz = 1 + rand() % (STRING_MAX-1);
for (j = 0; j < sz; j++)
{
B[i * (STRING_MAX + sizeof(int)) + sizeof(int) + j] = (char)(rand() % 57 + 65);
}
memcpy((char*)(B + i * (STRING_MAX + sizeof(int))), &sz, sizeof(int));
}
// 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, SET2_SIZE*(STRING_MAX + sizeof(int)), NULL, &ret);
cl_mem temp_set = clCreateBuffer(context, CL_MEM_READ_WRITE, SET_SIZE * sizeof(int), NULL, &ret);
cl_mem results_set = clCreateBuffer(context, CL_MEM_WRITE_ONLY, SET2_SIZE * sizeof(int) * 2, 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, SET2_SIZE*(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_l = clCreateKernel(program, "levenshtein", &ret);
cl_kernel kernel_r = clCreateKernel(program, "reduce", &ret);
ret = clSetKernelArg(kernel_l, 0, sizeof(cl_mem), (void *)&set_a);
ret = clSetKernelArg(kernel_l, 1, sizeof(cl_mem), (void *)&set_b);
ret = clSetKernelArg(kernel_l, 2, sizeof(cl_mem), (void *)&temp_set);
ret = clSetKernelArg(kernel_r, 0, sizeof(cl_mem), (void *)&temp_set);
ret = clSetKernelArg(kernel_r, 1, sizeof(cl_mem), (void *)&results_set);
i = SET_SIZE;
ret = clSetKernelArg(kernel_r, 3, sizeof(cl_uint), (void *)&i);
size_t global_item_size = SET_SIZE;
size_t global_item2_size = 64;
size_t local_item_size = 64;
unsigned int *C = (int*)malloc(sizeof(int)*SET2_SIZE*2);
clock_t start = clock();
for (i=0;i < SET2_SIZE; i++)
{
// Execute
ret = clSetKernelArg(kernel_l, 3, sizeof(cl_uint), (void *)&i);
ret = clEnqueueNDRangeKernel(command_queue, kernel_l, 1, NULL, &global_item_size, &local_item_size, 0, NULL, NULL);
ret = clSetKernelArg(kernel_r, 2, sizeof(cl_uint), (void *)&i);
ret = clEnqueueNDRangeKernel(command_queue, kernel_r, 1, NULL, &global_item2_size, &local_item_size, 0, NULL, NULL);
ret = clFinish(command_queue);
}
ret = clEnqueueReadBuffer(command_queue, results_set, CL_TRUE, 0, SET2_SIZE * sizeof(int) * 2, C, 0, NULL, NULL);
clock_t end = clock();
exec_time += (double)(end - start) / CLOCKS_PER_SEC;
printf("execution time = %f seconds\n", exec_time);
// Print best
//for (i = 0; i < SET2_SIZE; i++)
// printf("%d = %u (index %u)\n", i, C[i*2], C[i*2+1]);
ret = clFlush(command_queue);
ret = clFinish(command_queue);
ret = clReleaseKernel(kernel_l);
ret = clReleaseKernel(kernel_r);
ret = clReleaseProgram(program);
ret = clReleaseMemObject(set_a);
ret = clReleaseMemObject(set_b);
ret = clReleaseMemObject(results_set);
ret = clReleaseMemObject(temp_set);
ret = clReleaseCommandQueue(command_queue);
ret = clReleaseContext(context);
free(A);
free(B);
free(C);
return 0;
}
kernel.cl:
#define MIN3(a, b, c) ((a) < (b) ? ((a) < (c) ? (a) : (c)) : ((b) < (c) ? (b) : (c)))
#define MAX_STRING 200
__kernel void levenshtein(__global const uchar *A, __global const uchar *B, __global uint *C, uint index)
{
uint s1len, s2len, x, y, lastdiag, olddiag;
uint temp, temp2;
__local uint column[MAX_STRING];
int i = get_global_id(0);
temp = i * (MAX_STRING+sizeof(uint));
s1len = (A[temp+3] << 24) | (A[temp+2] << 16) | (A[temp+1] << 8) | A[temp];
temp = index * (MAX_STRING+sizeof(uint));
s2len = (B[temp+3] << 24) | (B[temp+2] << 16) | (B[temp+1] << 8) | B[temp];
for (y = 1; y <= s1len; y++)
column[y] = y;
for (x = 1; x <= s2len; x++)
{
column[0] = x;
for (y = 1, lastdiag = x-1; y <= s1len; y++)
{
olddiag = column[y];
column[y] = MIN3(column[y] + 1, column[y-1] + 1, lastdiag + (A[temp + sizeof(int) + y-1] == B[temp + sizeof(int) + x - 1] ? 0 : 1));
lastdiag = olddiag;
}
}
C[i] = column[s1len];
}
// THE WORST IMPLEMENTATION EVER!!!
__kernel void reduce(__global uint *B, __global uint *C, uint index, uint set_size)
{
uint i;
uint ind=0,dist=10000;
__local uint res[64];
__local uint off[64];
for (i=(set_size/64) * get_local_id(0); i < ((set_size/64)*(get_local_id(0)+1)); i++)
{
if (B[i] < dist)
{
dist = B[i];
ind = i;
}
}
res[get_local_id(0)] = dist;
off[get_local_id(0)] = ind;
barrier(CLK_LOCAL_MEM_FENCE);
if (get_local_id(0) == 0)
{
ind = 0;
dist = 1000;
for (i=0; i<64; i++)
{
if (res[i] < dist)
{
dist = res[i];
ind = off[i];
}
}
C[index*2] = dist;
C[index*2+1] = ind;
}
}
резултат:
execution time = 35.197846 seconds
Като за силно неоптимизиран код, толкова.
P.S забавно 1 срещу 100 хиляди за 0.03 секунди, 1000 срещу 100 хиляди за 35 секунди, при това с тоя потресаващо тъп reduce кернел, очевидно доста съм го подценил, явно нвидията има некакви скрити "costs" когато изпълнява за първи път кернела, при АМД смътни такива спомени имах от същото, но нямам идея на какво точно се дължи.
0 1 2 3 4 ...6 7 8 9 10 ...14 15 16 17 18 ...24 25 26 27 28 ...32 33 34 35 36