Обчислюємо «магічні квадрати» за допомогою GPU

Привіт habr.

Тема «магічних квадратів» досить цікава, оскільки з одного боку, вони відомі ще з давнини, з іншого боку, обчислення «магічного квадрата» навіть сьогодні являє собою досить непросту обчислювальну задачу. Нагадаємо, щоб побудувати «магічний квадрат» NxN, потрібно вписати числа 1..N*N так, щоб сума його горизонталей, вертикалей та діагоналей була дорівнює одному й тому ж числу. Якщо просто перебрати число всіх варіантів розстановки цифр для квадрата 4х4, то отримаємо 16! = 20 922 789 888 000 варіантів.

Подумаємо, як це зробити більш ефективно.

Для початку, повторимо умову задачі. Потрібно розставити числа в квадраті так, щоб вони не повторювалися, і сума горизонталей, вертикалей та діагоналей була дорівнює одному й тому ж числу.

Легко довести, що ця сума завжди однакова, і обчислюється за формулою для будь-якого n:

Ми будемо розглядати квадрати 4х4, так що сума = 34.
Позначимо всі змінні черех X, наш квадрат буде мати такий вигляд:

Перше і очевидне, властивість: т. к. сума квадрата відома, крайні стоблцы можна виразити через інші 3:
X14 = S - X11 - X12 - X13
X24 = S - X21 - X22 - X23
...
X41 = S - X11 - X21 - X31

Таким чином, квадрат 4х4 фактично перетворюється в квадрат 3х3, що зменшує кількість варіантів перебору з 16! до 9!, тобто в 57млн раз. Знаючи це, приступаємо до написання коду, подивимося наскільки складний такий перебір для сучасних комп’ютерів.

С++ — однопотоковий варіант

Принцип програми досить простий. Беремо безліч чисел 1..16 і цикл for з цього безлічі, це буде х11. Потім беремо друге безліч, що складається з першого за винятком числа x11, і так далі.

Приблизний вигляд програми виглядає так:

int squares = 0;
int digits[] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 };
Set mset(digits, digits + N*N);
for (int x11 = 1; x11 <= MAX; x11++) {
 Set set12(mset); set12.erase(x11);
 for (SetIterator it12 = set12.begin(); it12 != set12.end(); it12++) {
 int x12 = *it12;
 Set set13(set12); set13.erase(x12);
 for (SetIterator it13 = set13.begin(); it13 != set13.end(); it13++) {
 int x13 = *it13;
 int x14 = S - x11 - x12 - x13;
 if (x14 < 1 || x14 > MAX) continue;
 if (x14 == x11 || x14 == x12 || x14 == x13) continue;
...
 int sh1 = x11 + x12 + x13 + x14, sh2 = x21 + x22 + x23 + x24, sh3 = x31 + x32 + x33 + x34, sh4 = x41 + x42 + x43 + x44;
 int sv1 = x11 + x21 + x31 + x41, sv2 = x12 + x22 + x32 + x42, sv3 = x13 + x23 + x33 + x43, sv4 = x14 + x24 + x34 + x44;
 int sd1 = x11 + x22 + x33 + x44, sd2 = x14 + x23 + x32 + x41;
 if (sh1 != S || sh2 != S || sh3 != S || sh4 != S || sv1 != S || sv2 != S || sv3 != S || sv4 != S || sd1 != S || sd2 != S)
continue;
 // Якщо числа пройшли всі перевірки на перетину, то квадрат знайдено
 printf("%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %dn", x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34, x41, x42, x43, x44);
squares++;
}
}
printf("CNT: %dn", squares);

Повний текст програми можна знайти під спойлером.
Вихідний текст цілком

#include <set>
#include <stdio.h>
#include <ctime>
#include "stdafx.h"

typedef std::set<int> Set;
typedef Set::iterator SetIterator;

#define N 4
#define MAX (N*N)
#define S 34

int main(int argc, char *argv[])
{
 // x11 x12 x13 x14 
 // x21 x22 x23 x24 
 // x31 x32 x33 x34 
 // x41 x42 x43 x44 

 const clock_t begin_time = clock();

 int squares = 0;
 int digits[] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 };
 Set mset(digits, digits + N*N);
 for (int x11 = 1; x11 <= MAX; x11++) {
 Set set12(mset); set12.erase(x11);
 for (SetIterator it12 = set12.begin(); it12 != set12.end(); it12++) {
 int x12 = *it12;
 Set set13(set12); set13.erase(x12);
 for (SetIterator it13 = set13.begin(); it13 != set13.end(); it13++) {
 int x13 = *it13;
 int x14 = S - x11 - x12 - x13;
 if (x14 < 1 || x14 > MAX) continue;
 if (x14 == x11 || x14 == x12 || x14 == x13) continue;

 Set set21(set13); set21.erase(x13); set21.erase(x14);
 for (SetIterator it21 = set21.begin(); it21 != set21.end(); it21++) {
 int x21 = *it21;
 Set set22(set21); set22.erase(x21);
 for (SetIterator it22 = set22.begin(); it22 != set22.end(); it22++) {
 int x22 = *it22;
 Set set23(set22); set23.erase(x22);
 for (SetIterator it23 = set23.begin(); it23 != set23.end(); it23++) {
 int x23 = *it23, x24 = S - x21 - x22 - x23;
 if (x24 < 1 || x24 > MAX) continue;
 if (x24 == x11 || x24 == x12 || x24 == x13 || x24 == x14 || x24 == x21 || x24 == x22 || x24 == x23) continue;

 Set set31(set23);
 set31.erase(x23); set31.erase(x24);
 for (SetIterator it31 = set31.begin(); it31 != set31.end(); it31++) {
 int x31 = *it31;
 Set set32(set31); set32.erase(x31);
 for (SetIterator it32 = set32.begin(); it32 != set32.end(); it32++) {
 int x32 = *it32;
 Set set33(set32); set33.erase(x32);
 for (SetIterator it33 = set33.begin(); it33 != set33.end(); it33++) {
 int x33 = *it33, x34 = S - x31 - x32 - x33;
 if (x34 < 1 || x34 > MAX) continue;
 if (x34 == x11 || x34 == x12 || x34 == x13 || x34 == x14 || x34 == x21 || x34 == x22 || x34 == x23 || x34 == x24 || x34 == x31 || x34 == x32 || x34 == x33) continue;

 int x41 = S - x11 - x21 - x31, x42 = S - x12 - x22 - x32, x43 = S - x13 - x23 - x33, x44 = S - x14 - x24 - x34;
 if (x41 < 1 || x41 > MAX || x42 < 1 || x42 > MAX || x43 < 1 || x43 > MAX || x44 < 1 || x41 > MAX) continue;

 if (x41 == x11 || x41 == x12 || x41 == x13 || x41 == x14 || x41 == x21 || x41 == x22 || x41 == x23 || x41 == x24 ||
 x41 == x31 || x41 == x32 || x41 == x33 || x41 == x34)
continue;
 if (x42 == x11 || x42 == x12 || x42 == x13 || x42 == x14 || x42 == x21 || x42 == x22 || x42 == x23 || x42 == x24 ||
 x42 == x31 || x42 == x32 || x42 == x33 || x42 == x34 || x42 == x41)
continue;
 if (x43 == x11 || x43 == x12 || x43 == x13 || x43 == x14 || x43 == x21 || x43 == x22 || x43 == x23 || x43 == x24 ||
 x43 == x31 || x43 == x32 || x43 == x33 || x43 == x34 || x43 == x41 || x43 == x42)
continue;
 if (x44 == x11 || x44 == x12 || x44 == x13 || x44 == x14 || x44 == x21 || x44 == x22 || x44 == x23 || x44 == x24 ||
 x44 == x31 || x44 == x32 || x44 == x33 || x44 == x34 || x44 == x41 || x44 == x42 || x44 == x43)
continue;

 int sh1 = x11 + x12 + x13 + x14, sh2 = x21 + x22 + x23 + x24, sh3 = x31 + x32 + x33 + x34, sh4 = x41 + x42 + x43 + x44;
 int sv1 = x11 + x21 + x31 + x41, sv2 = x12 + x22 + x32 + x42, sv3 = x13 + x23 + x33 + x43, sv4 = x14 + x24 + x34 + x44;
 int sd1 = x11 + x22 + x33 + x44, sd2 = x14 + x23 + x32 + x41;
 if (sh1 != S || sh2 != S || sh3 != S || sh4 != S || sv1 != S || sv2 != S || sv3 != S || sv4 != S || sd1 != S || sd2 != S)
continue;

 printf("%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %dn", x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34, x41, x42, x43, x44);
squares++;
}
}
}
}
}
}
}
}
}

 printf("CNT: %dn", squares);

 float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
 printf("T = %.2fsn", diff_t);

 return 0;
}

Результат: всього було знайдено 7040 варіантів «магічних квадратів» 4х4, а час пошуку склало 102с.

До речі, цікаво перевірити, чи є в списку квадратів той самий, який зображений на гравюрі Дюрера. Зрозуміло є, т. к. програма виводить всі квадрати розмірності 4х4:

Читайте також  Мене оточують ідіоти або як працювати в команді

Не можна не відзначити, що Дюрер вставив квадрат зображення не просто так, цифри 1514 також позначають рік створення гравюри.

Як можна бачити, програма працює (відзначаємо завдання як verified at 1514 by Albrecht Дюрера;), проте час виконання не так і мало для комп’ютера з процесором Core i7. Очевидно, що програма виконується в один потік, і доцільно задіяти всі інші ядра.

С++ — багатопотоковий варіант

Переписати програму з використанням потоків в принципі нескладно, хоча і трохи громіздко. На щастя, є майже забутий сьогодні варіант — використання підтримки OpenMP (Open Multi-Processing). Ця технологія існує аж з 1998р, і дозволяє директивами процесора вказати компілятору, які частини програми виконувати паралельно. Підтримка OpenMP є і в Visual Studio, так що для перетворення програми в багатопотокову, досить додати код лише один рядок:

int squares = 0;
#pragma omp parallel for reduction(+: squares)
for (int x11 = 1; x11 <= MAX; x11++) {
...
}
printf("CNT: %dn", squares);

Директива #pragma omp parallel for вказує, що наступний цикл for можна виконувати паралельно, а додатковий параметр squares задає ім’я змінної, яка буде спільною для паралельних потоків (без цього інкремент працює некоректно).

Результат очевидний: час виконання скоротилося з 102с до 18с.

Вихідний текст цілком

#include <set>
#include <stdio.h>
#include <ctime>
#include "stdafx.h"

typedef std::set<int> Set;
typedef Set::iterator SetIterator;

#define N 4
#define MAX (N*N)
#define S 34

int main(int argc, char *argv[])
{
 // x11 x12 x13 x14 
 // x21 x22 x23 x24 
 // x31 x32 x33 x34 
 // x41 x42 x43 x44 

 const clock_t begin_time = clock();

 int squares = 0;
 #pragma omp parallel for reduction(+: squares)
 for (int x11 = 1; x11 <= MAX; x11++) {
 int digits[] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 };
 Set mset(digits, digits + N*N);
 Set set12(mset); set12.erase(x11);
 for (SetIterator it12 = set12.begin(); it12 != set12.end(); it12++) {
 int x12 = *it12;
 Set set13(set12); set13.erase(x12);
 for (SetIterator it13 = set13.begin(); it13 != set13.end(); it13++) {
 int x13 = *it13;
 int x14 = S - x11 - x12 - x13;
 if (x14 < 1 || x14 > MAX) continue;
 if (x14 == x11 || x14 == x12 || x14 == x13) continue;

 Set set21(set13); set21.erase(x13); set21.erase(x14);
 for (SetIterator it21 = set21.begin(); it21 != set21.end(); it21++) {
 int x21 = *it21;
 Set set22(set21); set22.erase(x21);
 for (SetIterator it22 = set22.begin(); it22 != set22.end(); it22++) {
 int x22 = *it22;
 Set set23(set22); set23.erase(x22);
 for (SetIterator it23 = set23.begin(); it23 != set23.end(); it23++) {
 int x23 = *it23, x24 = S - x21 - x22 - x23;
 if (x24 < 1 || x24 > MAX) continue;
 if (x24 == x11 || x24 == x12 || x24 == x13 || x24 == x14 || x24 == x21 || x24 == x22 || x24 == x23) continue;

 Set set31(set23);
 set31.erase(x23); set31.erase(x24);
 for (SetIterator it31 = set31.begin(); it31 != set31.end(); it31++) {
 int x31 = *it31;
 Set set32(set31); set32.erase(x31);
 for (SetIterator it32 = set32.begin(); it32 != set32.end(); it32++) {
 int x32 = *it32;
 Set set33(set32); set33.erase(x32);
 for (SetIterator it33 = set33.begin(); it33 != set33.end(); it33++) {
 int x33 = *it33, x34 = S - x31 - x32 - x33;
 if (x34 < 1 || x34 > MAX) continue;
 if (x34 == x11 || x34 == x12 || x34 == x13 || x34 == x14 || x34 == x21 || x34 == x22 || x34 == x23 || x34 == x24 || x34 == x31 || x34 == x32 || x34 == x33) continue;

 int x41 = S - x11 - x21 - x31, x42 = S - x12 - x22 - x32, x43 = S - x13 - x23 - x33, x44 = S - x14 - x24 - x34;
 if (x41 < 1 || x41 > MAX || x42 < 1 || x42 > MAX || x43 < 1 || x43 > MAX || x44 < 1 || x41 > MAX) continue;

 if (x41 == x11 || x41 == x12 || x41 == x13 || x41 == x14 || x41 == x21 || x41 == x22 || x41 == x23 || x41 == x24 ||
 x41 == x31 || x41 == x32 || x41 == x33 || x41 == x34)
continue;
 if (x42 == x11 || x42 == x12 || x42 == x13 || x42 == x14 || x42 == x21 || x42 == x22 || x42 == x23 || x42 == x24 ||
 x42 == x31 || x42 == x32 || x42 == x33 || x42 == x34 || x42 == x41)
continue;
 if (x43 == x11 || x43 == x12 || x43 == x13 || x43 == x14 || x43 == x21 || x43 == x22 || x43 == x23 || x43 == x24 ||
 x43 == x31 || x43 == x32 || x43 == x33 || x43 == x34 || x43 == x41 || x43 == x42)
continue;
 if (x44 == x11 || x44 == x12 || x44 == x13 || x44 == x14 || x44 == x21 || x44 == x22 || x44 == x23 || x44 == x24 ||
 x44 == x31 || x44 == x32 || x44 == x33 || x44 == x34 || x44 == x41 || x44 == x42 || x44 == x43)
continue;

 int sh1 = x11 + x12 + x13 + x14, sh2 = x21 + x22 + x23 + x24, sh3 = x31 + x32 + x33 + x34, sh4 = x41 + x42 + x43 + x44;
 int sv1 = x11 + x21 + x31 + x41, sv2 = x12 + x22 + x32 + x42, sv3 = x13 + x23 + x33 + x43, sv4 = x14 + x24 + x34 + x44;
 int sd1 = x11 + x22 + x33 + x44, sd2 = x14 + x23 + x32 + x41;
 if (sh1 != S || sh2 != S || sh3 != S || sh4 != S || sv1 != S || sv2 != S || sv3 != S || sv4 != S || sd1 != S || sd2 != S)
continue;

 printf("%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %dn", x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34, x41, x42, x43, x44);
squares++;
}
}
}
}
}
}
}
}
}

 printf("CNT: %dn", squares);

 float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
 printf("T = %.2fsn", diff_t);

 return 0;
}

Це набагато краще — т. к. завдання практично ідеально распараллеливается (розрахунки у кожній гілці не залежать один від одного), менше приблизно в число разів, рівне кількості ядер процесора. Але, на жаль, принципово бпрольшего з цього коду не отримати, хоча якимись оптимизациями може і можна виграти кілька відсотків. Переходимо до більш важкої артилерії, розрахунками на GPU.

Обчислення за допомогою NVIDIA CUDA

Якщо не вдаватися в подробиці, то процес обчислень, що виконується на відеокарті, можна подати як кілька паралельних апаратних блоків (blocks), кожен з яких виконує кілька процесів (threads).

Для прикладу можна навести приклад функції додавання 2х векторів документації CUDA:

__global__
void add(int n, float *x, float *y)
{
 int index = threadIdx.x;
 int stride = blockDim.x;
 for (int i = index; i < n; i += stride)
 y[i] = x[i] + y[i];
}

Масиви x і y — загальні для всіх блоків, а сама функція таким чином виконується одночасно відразу на декількох процесорах. Ключ тут у паралелізм — процесори відеокарти набагато простіше ніж звичайний CPU, зате їх багато і вони орієнтовані саме на обробку числових даних.

Це те, що нам потрібно. Ми маємо матрицю чисел X11, X12,..,X44. Запустимо процес з 16 блоків, кожен з яких буде виконувати 16 процесів. Номером блоку буде відповідати число X11, номером процесу число X12, а сам код буде обчислювати всі можливі квадрати з вибраних X11 і X12. Все просто, але тут є одна тонкість — дані потрібно не тільки знайти, але і передати з відеокарти назад, для цього в нульовому елементі масиву будемо зберігати кількість знайдених квадратів.

Код основний виходить дуже простим:

#define N 4
#define SQ_MAX 8*1024
#define BLOCK_SIZE (SQ_MAX*N*N + 1)

int main(int argc,char *argv[])
{
 const clock_t begin_time = clock();

 int *results = (int*)malloc(BLOCK_SIZE*sizeof(int));
 results[0] = 0;

 int *gpu_out = NULL;
 cudaMalloc(&gpu_out, BLOCK_SIZE*sizeof(int));
 cudaMemcpy(gpu_out, results, BLOCK_SIZE*sizeof(int), cudaMemcpyHostToDevice);
 squares<<<MAX, MAX>>>(gpu_out);
 cudaMemcpy(results, gpu_out, BLOCK_SIZE*sizeof(int), cudaMemcpyDeviceToHost);

 // Print results
 int squares = results[0];
 for(int p=0; p<SQ_MAX && p<squares; p++) {
 int i = MAX*p + 1;
 printf("[%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d]n", 
 results[i], results[i+1], results[i+2], results[i+3], 
 results[i+4], results[i+5], results[i+6], results[i+7],
 results[i+8], results[i+9], results[i+10], results[i+11],
 results[i+12], results[i+13], results[i+14], results[i+15]);
}
 printf ("CNT: %dn", squares);

 float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
 printf("T = %.2fsn", diff_t);

cudaFree(gpu_out);
free(results);

 return 0;
}

Ми виділяємо блок пам’яті на відеокарті з допомогою cudaMalloc, запускаємо функцію squares, вказавши їй 2 параметра 16,16 (число блоків і число потоків), відповідні перебираемым чисел 1..16, потім копіюємо дані назад через cudaMemcpy.

Читайте також  OpenSceneGraph: Основи роботи з геометрією сцени

Сама функція squares по суті повторює код з попередньої частини, з тією різницею, що приріст кількості знайдених квадратів робиться з допомогою atomicAdd — це гарантує, що змінна буде коректно змінюватися при одночасних зверненнях.
Вихідний код повністю

// Compile:
// nvcc -o magic4_gpu.exe magic4_gpu.cu

#include <stdio.h>
#include <ctime>

#define N 4
#define MAX (N*N)
#define SQ_MAX 8*1024
#define BLOCK_SIZE (SQ_MAX*N*N + 1)
#define S 34

// Magic square:
// x11 x12 x13 x14 
// x21 x22 x23 x24 
// x31 x32 x33 x34 
// x41 x42 x43 x44 

__global__ void squares(int *res_array) {
 int index1 = blockIdx.x, index2 = threadIdx.x;
 if (index1 + 1 > MAX || index2 + 1 > MAX) return;

 const int x11 = index1+1, x12 = index2+1;
 for(int x13=1; x13<=MAX; x13++) { 
 if (x13 == x11 || x13 == x12)
continue;
 int x14 = S - x11 - x12 - x13;
 if (x14 < 1 || x14 > MAX) continue;
 if (x14 == x11 || x14 == x12 || x14 == x13)
continue;
 for(int x21=1; x21<=MAX; x21++) { 
 if (x21 == x11 || x21 == x12 || x21 == x13 || x21 == x14)
continue;
 for(int x22=1; x22<=MAX; x22++) {
 if (x22 == x11 || x22 == x12 || x22 == x13 || x22 == x14 || x22 == x21)
continue;
 for(int x23=1; x23<=MAX; x23++) {
 int x24 = S - x21 - x22 - x23;
 if (x24 < 1 || x24 > MAX) continue;
 if (x23 == x11 || x23 == x12 || x23 == x13 || x23 == x14 || x23 == x21 || x23 == x22)
continue;
 if (x24 == x11 || x24 == x12 || x24 == x13 || x24 == x14 || x24 == x21 || x24 == x22 || x24 == x23)
continue;
 for(int x31=1; x31<=MAX; x31++) { 
 if (x31 == x11 || x31 == x12 || x31 == x13 || x31 == x14 || x31 == x21 || x31 == x22 || x31 == x23 || x31 == x24)
continue;
 for(int x32=1; x32<=MAX; x32++) {
 if (x32 == x11 || x32 == x12 || x32 == x13 || x32 == x14 || x32 == x21 || x32 == x22 || x32 == x23 || x32 == x24 || x32 == x31)
continue;
 for(int x33=1; x33<=MAX; x33++) {
 int x34 = S - x31 - x32 - x33;
 if (x34 < 1 || x34 > MAX) continue;
 if (x33 == x11 || x33 == x12 || x33 == x13 || x33 == x14 || x33 == x21 || x33 == x22 || x33 == x23 || x33 == x24 || x33 == x31 || x33 == x32)
continue;
 if (x34 == x11 || x34 == x12 || x34 == x13 || x34 == x14 || x34 == x21 || x34 == x22 || x34 == x23 || x34 == x24 || x34 == x31 || x34 == x32 || x34 == x33)
continue;

 const int x41 = S - x11 - x21 - x31, x42 = S - x12 - x22 - x32, x43 = S - x13 - x23 - x33, x44 = S - x14 - x24 - x34;
 if (x41 < 1 || x41 > MAX || x42 < 1 || x42 > MAX || x43 < 1 || x43 > MAX || x44 < 1 || x44 > MAX) 
continue;
 if (x41 == x11 || x41 == x12 || x41 == x13 || x41 == x14 || x41 == x21 || x41 == x22 || x41 == x23 || x41 == x24 ||
 x41 == x31 || x41 == x32 || x41 == x33 || x41 == x34)
continue;
 if (x42 == x11 || x42 == x12 || x42 == x13 || x42 == x14 || x42 == x21 || x42 == x22 || x42 == x23 || x42 == x24 ||
 x42 == x31 || x42 == x32 || x42 == x33 || x42 == x34 || x42 == x41)
continue;
 if (x43 == x11 || x43 == x12 || x43 == x13 || x43 == x14 || x43 == x21 || x43 == x22 || x43 == x23 || x43 == x24 || 
 x43 == x31 || x43 == x32 || x43 == x33 || x43 == x34 || x43 == x41 || x43 == x42)
continue;
 if (x44 == x11 || x44 == x12 || x44 == x13 || x44 == x14 || x44 == x21 || x44 == x22 || x44 == x23 || x44 == x24 || 
 x44 == x31 || x44 == x32 || x44 == x33 || x44 == x34 || x44 == x41 || x44 == x42 || x44 == x43)
continue;

 int sh1 = x11 + x12 + x13 + x14, sh2 = x21 + x22 + x23 + x24, sh3 = x31 + x32 + x33 + x34, sh4 = x41 + x42 + x43 + x44;
 int sv1 = x11 + x21 + x31 + x41, sv2 = x12 + x22 + x32 + x42, sv3 = x13 + x23 + x33 + x43, sv4 = x14 + x24 + x34 + x44;
 int sd1 = x11 + x22 + x33 + x44, sd2 = x14 + x23 + x32 + x41;
 if (sh1 != S || sh2 != S || sh3 != S || sh4 != S || sv1 != S || sv2 != S || sv3 != S || sv4 != S || sd1 != S || sd2 != S)
continue;

 // Square found: save in array (MAX numbers for each square)
 int p = atomicAdd(res_array, 1);
 if (p >= SQ_MAX) continue;

 int i = MAX*p + 1;
 res_array[i] = x11; res_array[i+1] = x12; res_array[i+2] = x13; res_array[i+3] = x14; 
 res_array[i+4] = x21; res_array[i+5] = x22; res_array[i+6] = x23; res_array[i+7] = x24; 
 res_array[i+8] = x31; res_array[i+9] = x32; res_array[i+10] = x33; res_array[i+11] = x34; 
 res_array[i+12]= x41; res_array[i+13]= x42; res_array[i+14] = x43; res_array[i+15] = x44; 

 // Warning: printf from kernel makes calculation 2-3x slower
 // printf("%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %dn", x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34, x41, x42, x43, x44);
}
}
}
}
}
}
}
}

int main(int argc,char *argv[])
{
 int *gpu_out = NULL;
 cudaMalloc(&gpu_out, BLOCK_SIZE*sizeof(int));

 const clock_t begin_time = clock();

 int *results = (int*)malloc(BLOCK_SIZE*sizeof(int));
 results[0] = 0;
 cudaMemcpy(gpu_out, results, BLOCK_SIZE*sizeof(int), cudaMemcpyHostToDevice);

 squares<<<MAX, MAX>>>(gpu_out);

 cudaMemcpy(results, gpu_out, BLOCK_SIZE*sizeof(int), cudaMemcpyDeviceToHost);

 // Print results
 int squares = results[0];
 for(int p=0; p<SQ_MAX && p<squares; p++) {
 int i = MAX*p + 1;
 printf("[%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d]n", results[i], results[i+1], results[i+2], results[i+3],
 results[i+4], results[i+5], results[i+6], results[i+7],
 results[i+8], results[i+9], results[i+10], results[i+11],
 results[i+12], results[i+13], results[i+14], results[i+15]);
}
 printf ("CNT: %dn", squares);

 float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
 printf("T = %.2fsn", diff_t);

cudaFree(gpu_out);
free(results);

 return 0;
}

Результат не потребує коментарів — час виконання склало 2.7, що приблизно в 30 разів краще початкового однопотокового варіанти:

Як підказали в коментарях, можна задіяти ще більше апаратних блоків відеокарти, так що був випробуваний варіант з 256 блоків. Зміна коду мінімально:

__global__ void squares(int *res_array) {
 int index1 = blockIdx.x/MAX, index2 = blockIdx.x%MAX;
...
}

squares<<<MAX*MAX, 1>>>(gpu_out);

Це скоротило час ще в 2 рази, до 1.2. Далі, на кожному блоці можна запустити 16 процесів, що дає найкращий час 0.44 з.
Остаточний варіант коду

#include <stdio.h>
#include <ctime>

#define N 4
#define MAX (N*N)
#define SQ_MAX 8*1024
#define BLOCK_SIZE (SQ_MAX*N*N + 1)
#define S 34

// Magic square:
// x11 x12 x13 x14 
// x21 x22 x23 x24 
// x31 x32 x33 x34 
// x41 x42 x43 x44 

__global__ void squares(int *res_array) {
 int index1 = blockIdx.x/MAX, index2 = blockIdx.x%MAX, index3 = threadIdx.x;
 if (index1 + 1 > MAX || index2 + 1 > MAX || index3 + 1 > MAX) return;

 const int x11 = index1+1, x12 = index2+1, x13 = index3+1;
 if (x13 == x11 || x13 == x12)
return;
 int x14 = S - x11 - x12 - x13;
 if (x14 < 1 || x14 > MAX) return;
 if (x14 == x11 || x14 == x12 || x14 == x13)
return;
 for(int x21=1; x21<=MAX; x21++) { 
 if (x21 == x11 || x21 == x12 || x21 == x13 || x21 == x14)
continue;
 for(int x22=1; x22<=MAX; x22++) {
 if (x22 == x11 || x22 == x12 || x22 == x13 || x22 == x14 || x22 == x21)
continue;
 for(int x23=1; x23<=MAX; x23++) {
 int x24 = S - x21 - x22 - x23;
 if (x24 < 1 || x24 > MAX) continue;
 if (x23 == x11 || x23 == x12 || x23 == x13 || x23 == x14 || x23 == x21 || x23 == x22)
continue;
 if (x24 == x11 || x24 == x12 || x24 == x13 || x24 == x14 || x24 == x21 || x24 == x22 || x24 == x23)
continue;
 for(int x31=1; x31<=MAX; x31++) { 
 if (x31 == x11 || x31 == x12 || x31 == x13 || x31 == x14 || x31 == x21 || x31 == x22 || x31 == x23 || x31 == x24)
continue;
 for(int x32=1; x32<=MAX; x32++) {
 if (x32 == x11 || x32 == x12 || x32 == x13 || x32 == x14 || x32 == x21 || x32 == x22 || x32 == x23 || x32 == x24 || x32 == x31)
continue;
 for(int x33=1; x33<=MAX; x33++) {
 int x34 = S - x31 - x32 - x33;
 if (x34 < 1 || x34 > MAX) continue;
 if (x33 == x11 || x33 == x12 || x33 == x13 || x33 == x14 || x33 == x21 || x33 == x22 || x33 == x23 || x33 == x24 || x33 == x31 || x33 == x32)
continue;
 if (x34 == x11 || x34 == x12 || x34 == x13 || x34 == x14 || x34 == x21 || x34 == x22 || x34 == x23 || x34 == x24 || x34 == x31 || x34 == x32 || x34 == x33)
continue;

 const int x41 = S - x11 - x21 - x31, x42 = S - x12 - x22 - x32, x43 = S - x13 - x23 - x33, x44 = S - x14 - x24 - x34;
 if (x41 < 1 || x41 > MAX || x42 < 1 || x42 > MAX || x43 < 1 || x43 > MAX || x44 < 1 || x44 > MAX) 
continue;
 if (x41 == x11 || x41 == x12 || x41 == x13 || x41 == x14 || x41 == x21 || x41 == x22 || x41 == x23 || x41 == x24 ||
 x41 == x31 || x41 == x32 || x41 == x33 || x41 == x34)
continue;
 if (x42 == x11 || x42 == x12 || x42 == x13 || x42 == x14 || x42 == x21 || x42 == x22 || x42 == x23 || x42 == x24 ||
 x42 == x31 || x42 == x32 || x42 == x33 || x42 == x34 || x42 == x41)
continue;
 if (x43 == x11 || x43 == x12 || x43 == x13 || x43 == x14 || x43 == x21 || x43 == x22 || x43 == x23 || x43 == x24 || 
 x43 == x31 || x43 == x32 || x43 == x33 || x43 == x34 || x43 == x41 || x43 == x42)
continue;
 if (x44 == x11 || x44 == x12 || x44 == x13 || x44 == x14 || x44 == x21 || x44 == x22 || x44 == x23 || x44 == x24 || 
 x44 == x31 || x44 == x32 || x44 == x33 || x44 == x34 || x44 == x41 || x44 == x42 || x44 == x43)
continue;

 int sh1 = x11 + x12 + x13 + x14, sh2 = x21 + x22 + x23 + x24, sh3 = x31 + x32 + x33 + x34, sh4 = x41 + x42 + x43 + x44;
 int sv1 = x11 + x21 + x31 + x41, sv2 = x12 + x22 + x32 + x42, sv3 = x13 + x23 + x33 + x43, sv4 = x14 + x24 + x34 + x44;
 int sd1 = x11 + x22 + x33 + x44, sd2 = x14 + x23 + x32 + x41;
 if (sh1 != S || sh2 != S || sh3 != S || sh4 != S || sv1 != S || sv2 != S || sv3 != S || sv4 != S || sd1 != S || sd2 != S)
continue;

 // Square found: save in array (MAX numbers for each square)
 int p = atomicAdd(res_array, 1);
 if (p >= SQ_MAX) continue;

 int i = MAX*p + 1;
 res_array[i] = x11; res_array[i+1] = x12; res_array[i+2] = x13; res_array[i+3] = x14; 
 res_array[i+4] = x21; res_array[i+5] = x22; res_array[i+6] = x23; res_array[i+7] = x24; 
 res_array[i+8] = x31; res_array[i+9] = x32; res_array[i+10] = x33; res_array[i+11] = x34; 
 res_array[i+12]= x41; res_array[i+13]= x42; res_array[i+14] = x43; res_array[i+15] = x44; 

 // Warning: printf from kernel makes calculation 2-3x slower
 // printf("%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %dn", x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34, x41, x42, x43, x44);
}
}
}
}
}
}
}

int main(int argc,char *argv[])
{
 int *gpu_out = NULL;
 cudaMalloc(&gpu_out, BLOCK_SIZE*sizeof(int));

 const clock_t begin_time = clock();

 int *results = (int*)malloc(BLOCK_SIZE*sizeof(int));
 results[0] = 0;
 cudaMemcpy(gpu_out, results, BLOCK_SIZE*sizeof(int), cudaMemcpyHostToDevice);

 squares<<<MAX*MAX, MAX>>>(gpu_out);

 cudaMemcpy(results, gpu_out, BLOCK_SIZE*sizeof(int), cudaMemcpyDeviceToHost);

 // Print results
 int squares = results[0];
 for(int p=0; p<SQ_MAX && p<squares; p++) {
 int i = MAX*p + 1;
 printf("[%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d]n", results[i], results[i+1], results[i+2], results[i+3],
 results[i+4], results[i+5], results[i+6], results[i+7],
 results[i+8], results[i+9], results[i+10], results[i+11],
 results[i+12], results[i+13], results[i+14], results[i+15]);
}
 printf ("CNT: %dn", squares);

 float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
 printf("T = %.2fsn", diff_t);

cudaFree(gpu_out);
free(results);

 return 0;
}

Швидше за все, це далеко не ідеал, наприклад, можна запустити ще більше блоків на GPU, але це зробить код більш заплутаним і складним для розуміння. І зрозуміло, розрахунки даються не «безкоштовно» — при завантаженому GPU інтерфейс Windows починає помітно пригальмовувати, а енергоспоживання комп’ютера збільшується практично в 2 рази, з 65 до 130Вт.

Читайте також  Імплементація катсцен і послідовностей дій в іграх

Висновок

Задача знаходження «магічних квадратів» виявилася технічно вельми цікавою, і в той же час непростий. Навіть з обчисленнями на GPU пошук всіх квадратів 5х5 може зайняти кілька годин, а оптимізацію для пошуку магічних квадратів розмірності 7х7 і вище, ще належить зробити.

Математично і алгоритмічно, теж є кілька невирішених моментів:
— Залежність кількості «магічних квадратів» від N. Відомо що квадрата 2х2 не існує, квадратів 3х3 існує всього 8 (але по суті 1, решта його відбиття або повороти), квадратів 4х4 як ми з’ясували, 7040, але виключення поворотів або відображень в алгоритм поки не додано. Для великих розмірностей питання поки що відкрите.
— Виключення квадратів, є поворотами або відбитками вже знайденого.
— Швидкість і оптимізація алгоритму. На жаль, потестувати код на суперкомп’ютері або хоча б на NVIDIA Tesla можливості немає, якщо хтось може запустити, було б цікаво. Якщо у когось є ідеї щодо самого алгоритму, їх теж можна спробувати реалізувати. При бажанні можна навіть запустити розподілений проект з пошуку квадратів, якщо звичайно набереться достатньо число читаталей 😉

Про аналіз і властивостях самих магічних квадратів можна написати окрему статтю, якщо буде інтерес.

PS: ДО питання, який напевно буде, «а навіщо це треба». З точки зору витрат електроенергії обчислення магічних квадратів нічим не краще або гірше обчислення биткоинов, так що чому б і ні? До того ж, це цікава розминка для розуму і цікаве завдання в області прикладного програмування.

Степан Лютий

Обожнюю технології в сучасному світі. Хоча частенько і замислююся над тим, як далеко вони нас заведуть. Не те, щоб я прям і знаюся на ядрах, пікселях, коллайдерах і інших парсеках. Просто приходжу в захват від того, що може в творчому пориві вигадати людський розум.

You may also like...

Залишити відповідь

Ваша e-mail адреса не оприлюднюватиметься. Обов’язкові поля позначені *