===============================================================================================

1.  Расстояния между парами образов
    -------------------------------

В n-мерном евклидовом пространстве расстояние между двумя точками и определяется как .

Ставится задача построить матрицу вещественных чисел, каждый элемент которой должен показывать расстояние между -м и -м образами в выборке размерностью .

Каждый элемент матрицы может быть рассчитан независимо от остальных, поэтому построение матрицы можно распараллелить по принципу SIMD.

Естественным и тривиальным решением для универсальной многоядерной архитектуры является равномерное распределение работы по вычислению столбцов (строк) матрицы между вычислительными потоками. Рекомендуется выбирать равным количеству ядер процессора, чтобы эффективно задействовать все его вычислительные возможности. В таком случае потоки вычислят по столбцов (строк) матрицы , а потоки – по . Здесь – остаток от деления на .

Ниже приведен код реализации алгоритма вычисления матрицы (dImIm) и соответствующей матрицы квадратов расстояний (DImIm) на C++. В качестве средства распараллеливания используется OpenMP 2.0. Для распараллеливания итераций цикла достаточно указания директивы #pragma omp parallel for.

#ifdef _OPENMP

#pragma omp parallel for

#endif

for (int i = 0; i < N; i++)

{

for (int j = 0; j < N; j++)

{

float dstSqr = 0;

if (fmasklen < n && fmasklen > 0)

{

for (int k = 0; k < fmasklen; k++)

{

float delta = x[i + fmask[k] * MAX_IMAGE_CNT] - x[j + fmask[k] * MAX_IMAGE_CNT];

dstSqr += delta * delta;

}

}

else

{

for (int k = 0; k < n; k++)

{

float delta = x[i + k * MAX_IMAGE_CNT] - x[j + k * MAX_IMAGE_CNT];

dstSqr += delta * delta;

}

}

int idx = i + j * MAX_IMAGE_CNT;

dImIm[idx] = sqrtf(dstSqr);

DImIm[idx] = dstSqr;

}

}

Здесь мы отказываемся от очевидной оптимизации «вычисления треугольником», вытекающей из свойства , чтобы не нарушать общность и иметь возможность корректно сравнить результаты с полученными на графическом процессоре архитектуры NVIDIA CUDA.

Обозначенный выше подход для универсальных многоядерных процессоров прозрачно переносится на графические процессоры архитектуры NVIDIA CUDA с учетом того, что мы вольны выбрать . В таком случае каждый поток займется вычислением единственного столбца (строки) матрицы . Естественность обозначенного поточно-параллельного подхода для графических процессоров архитектуры NVIDIA CUDA вытекает из следующих одновременно выполняющихся условий:

• каждый поток выполняет одинаковый объем вычислений, что исключает нежелательные ветвления внутри warp-ов при окончании какого-либо из циклов по элементам столбца (строки), влекущие к простоям аппаратуры;

• естественно организуется когерентная запись результатов вычислений (потоки с последовательными номерами записывают результаты в последовательные слова области памяти, выделенной под результирующую матрицу), что позволяет в разы снижать количество транзакций доступа к глобальной памяти;

• потоки внутри блока могут использовать общую память мультипроцессора для хранения данных об образах, что позволяет получить существенный выигрыш за счет уменьшения в десятки раз количества обращений к глобальной памяти из каждого потока за данными об образах.

Опишем особенности реализации обозначенного поточно-параллельного подхода (рисунок 1). Вычислительная сетка является линейной, каждый блок содержит потоков. В целях увеличения производительности каждым блоком потоков используется общая память для кэширования атрибутов образов. Единовременное обновление содержимого этого кэша различными блоками должно производиться различными подмассивами (длиной ) массива атрибутов образов, чтобы единовременный доступ к глобальной памяти всеми активными warp-ами приходился на различные сегменты глобальной памяти. Такой подход существенно повышает производительность при работе с глобальной памятью для устройств вычислительной совместимости 1.x.

{width="6.4875in" height="5.033333333333333in"}

Рисунок 1 – Параллельное вычисление матрицы расстояний между парами образов на CUDA GPU

Ниже приведен код реализации алгоритма вычисления матрицы (dImIm) и соответствующей матрицы квадратов расстояний (DImIm) на CUDA C.

__constant__ int fmask[MAX_ATTR_CNT];

template <bool usemask> __global__ void g_distImIm(float* dImIm, float* DImIm, float* X, int N, int n)

{

extern __shared__ float subX[];

float* curX = (float*)(subX + (blockDim.x + threadIdx.x) * n);

int image = threadIdx.x + blockIdx.x * blockDim.x;

if (image < N)

{

//Скопировать атрибуты первичного образа в общую память.

for (int i = 0; i < n; i++) //Когерентное чтение.

curX[i] = X[image + (usemask ? fmask[i] : i) * MAX_IMAGE_CNT];

}

//Перебрать все «плитки» вторичных образов для первичного.

for (int i = 0; i < gridDim.x; i++)

{

//Скопировать «плитку» вторичных образов в общую память.

int offset = ((blockIdx.x + i) % gridDim.x) * blockDim.x;

if (threadIdx.x + offset < N)

for (int j = 0; j < n; j++) //Когерентное чтение.

subX[threadIdx.x + j * blockDim.x] = X[threadIdx.x + offset + (usemask ? fmask[j] : j) * MAX_IMAGE_CNT];

__syncthreads();

if (image < N)

{

//Рассчитать расстояния между первичным образом и всеми скопированными в общую память вторичными.

int copied = min(blockDim.x, N - offset);

for (int j = 0; j < copied; j++)

{

float dstSqr = 0;

for (int k = 0; k < n; k++)

{

float delta = curX[k] - subX[j + k * blockDim.x];

dstSqr += delta * delta;

}

int idx = image + (j + offset) * MAX_IMAGE_CNT;

dImIm[idx] = sqrtf(dstSqr); //Когерентная запись.

DImIm[idx] = dstSqr; //Когерентная запись.

}

}

__syncthreads();

}

}

Стремление организовать когерентный доступ к глобальной памяти вполне оправдано, поскольку некогерентный доступ в разы повышает количество дорогостоящих транзакций, способных свести на нет все преимущества параллельных вычислений на CUDA GPU.

Заметим, что для «вычислений треугольником» можно добиться когерентной записи результатов, только если каждый поток будет рассчитывать единственный элемент треугольной подматрицы . При этом warp-ы должны содержать потоки, вычисляющие элементы только одной строки (столбца) . В такой ситуации мы лишены возможности эффективно использовать общую память мультипроцессора для хранения данных об образах, поскольку каждому потоку требуется информация только о двух образах: -м (общем для всех потоков в блоке) и -м образе (различном для всех потоков в блоке). Если -й образ можно сохранить в общей памяти, то за -м образом каждый поток будет обращаться к глобальной памяти. Кроме того, нам потребуется организовать потоков для расчета всех элементов . Итак, если каждый поток будет обращаться за вторичным образом к глобальной памяти, будет запрошено образов. Первичный образ является общим на весь блок, обращение за его атрибутами происходит единожды. В итоге получаем обращений за атрибутами образов, где – размерность признакового пространства.

Затем нам потребуется заполнить вторичную треугольную подматрицу матрицы . Для этого потребуется чтение значений, рассчитанных ранее и расположенных в . Пусть все операции чтения выполнялись когерентно. Тогда – стоимость всех операций чтения (в тактах), где тактов – стоимость одной операции чтения в условиях когерентного доступа к памяти.

Для рассмотренного ранее подхода этот показатель равен , поскольку каждый блок читает данные всей выборки образов единожды. Здесь – количество потоков в блоке, – количество блоков в сетке.

Определим проигрыш подхода «вычислений треугольником» в аспекте чтений глобальной памяти. .

Ясно, что выигрыш подхода «вычислений треугольником» в аспекте вычисления клеток матрицы равен

Здесь тактов – стоимость вычитания, тактов – стоимость сложения, тактов – стоимость умножения, тактов – стоимость извлечения квадратного корня.

Выясним, может ли выигрыш от вычисления только подматрицы матрицы компенсировать дополнительные затраты на чтение глобальной памяти.

Согласно документации NVIDIA CUDA C Programming Guide для устройств вычислительной совместимости 2.0 имеем: такт, тактов. тактов.

Знаки коэффициентов при в выражениях для и не в пользу подхода «вычислений треугольником». Рассмотрим теперь коэффициенты при . Примем для определенности (рекомендованное значение NVIDIA).

. Таким образом, справедливо заключаем , то есть «вычисления треугольником» не дают выигрыша по сравнению с полными вычислениями и этот подход можно отвергнуть.

Расстояние от точки в признаковом пространстве до множества точек , относящихся к одному классу , может быть определено как расстояние от точки до ближайшей точки множества: .

Можно определить среднеквадратическое расстояние от точки до множества из элементов: .

Ставится задача построить матрицу вещественных чисел, каждый элемент которой должен показывать расстояние между -м классом и -м образом и матрицу , каждый элемент которой должен показывать среднеквадратическое расстояние между -м классом и -м образом. Размерность выборки образов – . Количество классов – .

Каждый столбец матриц может быть рассчитан независимо от остальных, поэтому построение матриц можно распараллелить по принципу SIMD.

Для вычисления -го столбца матрицы потребуется -й столбец рассчитанной ранее матрицы расстояний между образами , а для вычисления -го столбца матрицы потребуется -й столбец рассчитанной ранее матрицы квадратов расстояний между образами .

При вычислении столбцов и будут последовательно перебираться элементы соответствующих столбцов матриц и . В зависимости от класса, к которому принадлежит каждый элемент, будет обновляться соответствующий элемент результирующего столбца. Такой подход справедлив в силу ассоциативности операций нахождения минимума и сложения.

Как и в случае расстояний между образами, воспользуемся поточно-параллельным подходом для независимого вычисления столбцов матриц. Ниже приведен код реализации алгоритма вычисления матрицы (dImCl) и соответствующей матрицы среднеквадратических расстояний (DImCl) на C++. В качестве средства распараллеливания используется OpenMP 2.0. Для распараллеливания итераций цикла достаточно указания директивы #pragma omp parallel for.

#ifdef _OPENMP

#pragma omp parallel for

#endif

for (int i = 0; i < N; i++)

{

for (int j = 0; j < c; j++) //Проинициализировать начальные значения.

{

int idx = i + j * MAX_IMAGE_CNT;

dImCl[idx] = FLT_MAX;

DImCl[idx] = 0;

}

for (int j = 0; j < N; j++) //Определить минимумы и суммы квадратов расстояний от фиксированного образа до каждого класса.

{

int idx1 = i + clsInfo[j] * MAX_IMAGE_CNT;

int idx2 = i + j * MAX_IMAGE_CNT;

dImCl[idx1] = min(dImCl[idx1], dImIm[idx2]);

DImCl[idx1] += DImIm[idx2];

}

for (int j = 0; j < c; j++) //Учесть поправку на размер класса.

DImCl[i + j * MAX_IMAGE_CNT] /= clsDims[j];

}

Обозначенный выше подход для универсальных многоядерных процессоров прозрачно переносится на графические процессоры архитектуры NVIDIA CUDA с учетом того, что мы вольны выбрать количество потоков равным . В таком случае каждый поток займется вычислением единственного столбца матрицы и матрицы .

Алгоритм переносится лишь с одним изменением: промежуточные результаты инициализируются и накапливаются в общей памяти, а затем записываются в глобальную память. Это позволяет избежать некогерентного доступа к медленной глобальной памяти в процессе накопления результатов и организовать их когерентный перенос в глобальную память из общей по окончании вычислений.

Ниже приведен код реализации алгоритма вычисления матрицы (dImCl) и соответствующей матрицы среднеквадратических расстояний (DImCl) на CUDA C.

#define WideToFold(N, d) (N + (d - N % d) % d)

__global__ void g_distImCl(float* dImCl, float* DImCl, float* dImIm, float* DImIm, unsigned char* clsInfo, unsigned short* clsDims, int c, int N)

{

extern __shared__ int shared[];

unsigned char* sharedClsInfo = (unsigned char*)(shared);

unsigned short* sharedClsDims = (unsigned short*)(shared);

float* d = (float*)(shared + WideToFold(N * sizeof(unsigned char), sizeof(int)) / sizeof(int) + threadIdx.x * c * sizeof(float) / sizeof(int));

float* D = (float*)(d + blockDim.x * c);

//Скопировать информацию о принадлежности образов классам в общую память.

for (int i = threadIdx.x; i < N; i += blockDim.x)

sharedClsInfo[i] = clsInfo[i];

__syncthreads();

//Проинициализировать начальные значения.

for (int i = 0; i < c; i++)

{

d[i] = FLT_MAX;

D[i] = 0;

}

//Определить минимумы и суммы квадратов расстояний от фиксированного образа до каждого класса.

int image = threadIdx.x + blockIdx.x * blockDim.x;

if (image < N)

{

for (int i = 0; i < N; i++)

{

int idx = image + i * MAX_IMAGE_CNT;

int cl = sharedClsInfo[i];

d[cl] = min(d[cl], dImIm[idx]);

D[cl] += DImIm[idx];

}

}

//Скопировать размеры классов в общую память, полагая sizeof(unsigned short) * c <= sizeof(unsigned char) * N.

for (int i = threadIdx.x; i < c; i += blockDim.x)

sharedClsDims[i] = clsDims[i];

__syncthreads();

//Перенести результаты в глобальную память когерентно.

if (image < N)

{

for (int cl = 0; cl < c; cl++)

{

int idx = image + cl * MAX_IMAGE_CNT;

dImCl[idx] = d[cl];

//Учесть поправку на размер класса.

DImCl[idx] = D[cl] / sharedClsDims[cl];

}

}

}

Среднеквадратическое расстояние между фиксированной точкой и остальными точками множества определяется следующим образом: .

Это позволяет определить внутриклассовое расстояние как .

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

Ставится задача построить матрицу вещественных чисел такую, что каждый элемент главной диагонали должен показывать внутриклассовое расстояние для класса (), а все остальные элементы – расстояние между классами и (; ; ). Количество классов - .

Заметим, что , где – элемент построенной ранее матрицы среднеквадратических расстояний между образами и классами. Здесь – номер класса в общем списке классов, – номер образа в общем списке образов. Следовательно, . Кроме того, , где – также является элементом . Следовательно, .

Разумно разделить просмотр столбцов матрицы равномерно между всеми вычислительными потоками, как и в предыдущих случаях. При просмотре -го столбца матрицы будет сначала определено – номер класса -го образа, а затем каждое значение будет увеличено на в предположении, что было проинициализировано нулевым (). По окончании вычислений следует умножить каждый элемент главной диагонали на , а остальные элементы на , где – размер класса (,).

Ниже приведен код реализации алгоритма вычисления матрицы (DInInterCl) на C++. В качестве средства распараллеливания используется OpenMP 2.0. Избежать конфликтного доступа нескольких потоков к одной ячейке памяти позволяет директива OpenMP 2.0 #pragma omp atomic.

#ifdef _OPENMP

#pragma omp parallel

//Вход в параллельный регион, создание потоков.

{

#pragma omp for

#endif

//Параллельная инициализация нулевыми значениями.

for (int i = 0; i < c; i++)

{

for (int j = 0; j < c; j++)

DInInterCl[j + i * MAX_CLASS_CNT] = 0;

}

#ifdef _OPENMP

#pragma omp for

#endif

//Параллельный просмотр столбцов матрицы среднеквадратических расстояний между образами и классами.

for (int i = 0; i < N; i++)

{

int row = clsInfo[i];

for (int col = 0; col < c; col++)

{

#ifdef _OPENMP

#pragma omp atomic

#endif

//Атомарная операция инкремента.

DInInterCl[col + row * MAX_CLASS_CNT] += DImCl[i + col * MAX_IMAGE_CNT];

}

}

#ifdef _OPENMP

#pragma omp for

#endif

//Параллельная поправка на величину (размер класса – 1) или (размер класса)

for (int i = 0; i < c; i++)

{

for (int j = 0; j < c; j++)

{

bool inClass = (i == j);

DInInterCl[j + i * MAX_CLASS_CNT] /= max(1, clsDims[i] - (int)(inClass));

}

}

#ifdef _OPENMP

}

#endif

Строки (или столбцы) матрицы могут рассчитываться независимо. Для расчета каждой строки требуется просмотр строки матрицы среднеквадратических расстояний между образами и классами. Пусть каждый блок линейной вычислительной сетки будет рассчитывать соответствующую строку матрицы .

Поскольку графические процессоры архитектуры NVIDIA CUDA не имеют аппаратного решения проблемы конфликтного доступа многих потоков к одной ячейке памяти, разумно организовать вычисление индивидуального промежуточного вектора-результата каждым потоком в отдельности, а затем произвести редукцию к результирующему вектору. По-прежнему просмотр столбцов матрицы будет равномерно распределен между всеми вычислительными потоками.

Каждый поток единственного блока сконструирует промежуточный вектор-результат в общей памяти, просмотрев лишь часть столбцов . Затем, потоков осуществят редукцию промежуточных векторов результатов к результирующему вектору (каждый поток вычислит для соответствующего класса ).

Ниже приведен код реализации алгоритма вычисления матрицы (DInInterCl) на CUDA C.

__global__ void g_distInInterCl(float* DInInterCl, float* DImCl, unsigned char* clsInfo, unsigned short* clsDims, int c, int N)

{

extern __shared__ int shared[];

unsigned char* sharedClsInfo = (unsigned char*)(shared);

float* D = (float*)(shared + WideToFold(N * sizeof(unsigned char), sizeof(int)) / sizeof(int));

//Скопировать информацию о принадлежности образов классам в общую память.

for (int i = threadIdx.x; i < N; i += blockDim.x)

sharedClsInfo[i] = clsInfo[i];

__syncthreads();

//Проинициализировать промежуточные векторы-результаты.

for (int i = 0; i < c; i++)

D[i + threadIdx.x * c] = 0;

//Вычислить промежуточные векторы-результаты.

for (int i = threadIdx.x; i < N; i += blockDim.x)

D[sharedClsInfo[i] + threadIdx.x * c] += DImCl[i + blockIdx.x * MAX_IMAGE_CNT];

__syncthreads();

//Произвести редукцию к результирующему вектору с учетом поправки на величину (размер класса – 1) или (размер класса).

if (threadIdx.x < c) //Здесь мы полагаем blockDim.x >= c.

{

float Di = 0;

for (int i = 0; i < blockDim.x; i++)

Di += D[threadIdx.x + i * c];

bool inClass = (blockIdx.x == threadIdx.x);

DInInterCl[threadIdx.x + blockIdx.x * MAX_CLASS_CNT] = Di / max(1, clsDims[threadIdx.x] - (int)(inClass));

}

}

Заметим, что каждый поток с номером меньшим выполнит сложений, где первое слагаемое – количество сложений при вычислении промежуточного результата, а второе слагаемое – количество сложений при редукции. Время работы вычислительного ядра можно оптимизировать, уменьшив это количество сложений. Нетрудно убедиться в том, что минимума можно достичь при количестве потоков в блоке, равном .

Для оценки информативности множеств признаков применяется отношение среднего межклассового Евклидова расстояния в признаковом пространстве к среднему внутриклассовому расстоянию: . Здесь – элемент матрицы внутримножественных и межмножественных расстояний. Количество классов - .

Подход к распараллеливанию расчетов числителя и знаменателя будет идентичным для обеих рассматриваемых архитектур. Организуем потоки для просмотра каждой строки матрицы внутримножественных и межмножественных расстояний (для GPU – по одному потоку на строку, для CPU – равномерное распределение). Для каждой просмотренной стоки будет вычислен свой промежуточный результат. Затем произведем редукцию результатов каждого потока, пользуясь, как и ранее, ассоциативностью операции сложения.

Ниже приведен код реализации алгоритма вычисления матрицы информативности Q на C++. В качестве средства распараллеливания используется OpenMP 2.0. Осуществить редукцию индивидуальных промежуточных результатов позволяет директива OpenMP 2.0 #pragma omp parallel for reduction.

float interClass = 0;

float inClass = 0;

#pragma omp parallel for reduction(+:interClass,inClass)

for (int i = 0; i < c; i++)

{

for (int j = 0; j < c; j++)

{

if (i == j)

inClass += DInInterCl[j + i * MAX_CLASS_CNT];

else

interClass += DInInterCl[j + i * MAX_CLASS_CNT];

}

}

float Q = interClass / ((c - 1) * inClass);

Ниже приведен код реализации алгоритма вычисления матрицы информативности Q на CUDA C. Редукция осуществляется нулевым потоком в блоке. Алгоритм исполняется на одном мультипроцессоре, поскольку организация редукции результатов, вычисленных различными мультипроцессорами, затруднительна в силу отсутствия мультипроцессорных барьеров.

__global__ void g_calcInformativeness(float* DInInterCl, int c, float* Q)

{

extern __shared__ float dst[];

float* interClass = dst;

float* inClass = dst + blockDim.x;

interClass[threadIdx.x] = 0;

inClass[threadIdx.x] = 0;

for (int i = 0; i < c; i++)

{

bool interCl = (i != threadIdx.x);

bool inCl = (i == threadIdx.x);

interClass[threadIdx.x] += (int)(interCl) * DInInterCl[threadIdx.x + i * MAX_CLASS_CNT];

inClass[threadIdx.x] += (int)(inCl) * DInInterCl[threadIdx.x + i * MAX_CLASS_CNT];

}

__syncthreads();

if (threadIdx.x == 0)

{

float interCl = 0;

float inCl = 0;

for (int i = 0; i < blockDim.x; i++)

{

interCl += interClass[i];

inCl += inClass[i];

}

Q[0] = interCl / ((c - 1) * inCl);

}

}

1. Список ближайших соседей, матрица ошибок, классификация по критерию минимума расстояния

   1. Список ближайших соседей

Ставится задача построить вектор длиной , -й компонент которого будет содержать номер ближайшего соседа -го образа. Поскольку каждый компонент может быть вычислен независимо от остальных, существует возможность распараллелить построение по принципу SIMD. Предположим, что матрица расстояний между парами образов уже построена. Тогда вычисление номера ближайшего соседа -го образа сводится к нахождению минимума в -м столбце (или -й строке) матрицы расстояний (очевидно, что нулевые элементы на главной диагонали не могут определять ближайшего соседа).

Подход к распараллеливанию алгоритма будет идентичным для обеих рассматриваемых архитектур. Организуем потоки для просмотра каждого столбца матрицы и определения номера строки, содержащей минимальный элемент (для GPU – по одному потоку на столбец, для CPU – равномерное распределение). Номер этой строки и будет номером ближайшего соседа -го образа.

Ниже приведен код реализации алгоритма вычисления вектора номеров ближайших соседей (nearest) на C++. В качестве средства распараллеливания используется OpenMP 2.0. Для распараллеливания итераций цикла достаточно указания директивы #pragma omp parallel for.

#ifdef _OPENMP

#pragma omp parallel for

#endif

for (int i = 0; i < N; i++)

{

float min_val = FLT_MAX;

int min_idx = -1;

for (int j = 0; j < N; j++)

{

if ((j != i) && (DistMatr[j + i * MAX_IMAGE_CNT] < min_val))

{

min_val = DistMatr[j + i * MAX_IMAGE_CNT];

min_idx = j;

}

}

nearest[i] = (unsigned short)(min_idx);

}

Ниже приведен код реализации алгоритма вычисления вектора номеров ближайших соседей (nearest) на CUDA C.

__global__ void g_findNearest(float* DistMatr, unsigned short* nearest, int N)

{

int image = threadIdx.x + blockIdx.x * blockDim.x;

if (image < N)

{

float min_val = FLT_MAX;

int min_idx = -1;

//Перебрать все элементы столбца и найти номер минимального.

for (int i = 0; i < N; i++)

{

bool better_fact = (i != image) && (DistMatr[image + i * MAX_IMAGE_CNT] < min_val);

float better_val = min(min_val, DistMatr[image + i * MAX_IMAGE_CNT]);

min_val -= min_val * (int)(better_fact);

min_val += better_val * (int)(better_fact);

min_idx = min_idx + (i - min_idx) * (int)(better_fact);

}

nearest[image] = (unsigned short)(min_idx);

}

}

Матрица ошибок служит для обозначения ситуаций, когда образы и их ближайшие соседи не принадлежат одному классу. Элемент матрицы определяется следующим образом (здесь – номер класса, которому принадлежит образ с номером ):

.

Поскольку каждый элемент матрицы может быть вычислен независимо от остальных, существует возможность распараллелить построение по принципу SIMD. Предположим, что вектор номеров ближайших соседей образов уже построен.

Подход к распараллеливанию алгоритма будет идентичным для обеих рассматриваемых архитектур. Организуем потоки для расчета каждого столбца матрицы ошибок (для GPU – по одному потоку на столбец, для CPU – равномерное распределение).

Ниже приведен код реализации алгоритма вычисления матрицы ошибок (errors) на C++. В качестве средства распараллеливания используется OpenMP 2.0. Для распараллеливания итераций цикла достаточно указания директивы #pragma omp parallel for.

#ifdef _OPENMP

#pragma omp parallel for

#endif

for (int i = 0; i < N; i++)

{

for (int j = 0; j < N; j++)

{

bool error = (j == nearest[i]) && (clsInfo[i] != clsInfo[j]);

errors[i + j * MAX_IMAGE_CNT] = (unsigned char)(error);

}

}

Ниже приведен код реализации алгоритма вычисления матрицы ошибок (errors) на CUDA C.

__global__ void g_findErrors(unsigned short* nearest, unsigned char* errors, unsigned char* clsInfo, int N)

{

extern __shared__ int shared[];

unsigned char* sharedClsInfo = (unsigned char*)(shared);

unsigned short* sharedNearest = (unsigned short*)(shared + WideToFold(N * sizeof(unsigned char), sizeof(int)) / sizeof(int));

//Скопировать в общую память информацию о принадлежности образов классам и номера ближайших соседей.

for (int i = threadIdx.x; i < N; i += blockDim.x)

{

sharedClsInfo[i] = clsInfo[i];

sharedNearest[i] = nearest[i];

}

__syncthreads();

int image = threadIdx.x + blockIdx.x * blockDim.x;

if (image < N)

{

//Вычислить все элементы столбца.

for (int i = 0; i < N; i++)

{

bool error = (i == sharedNearest[image]) && (sharedClsInfo[image] != sharedClsInfo[i]);

errors[image + i * MAX_IMAGE_CNT] = (unsigned char)(error);

}

}

}

Если класс представлен эталонами , то функция расстояния между образом и классом записывается в виде: , , то есть расстоянием между классифицируемым образом и классом полагается наименьшее из расстояний между образом и эталонами этого класса. Правило классификации: образ зачисляется в класс , если условие выполняется для всех .

Рассмотрим выборку образов с известной классификацией по классам . Правило классификации, основанное на принципе ближайшего соседа: образ зачисляется в класс , к которому принадлежит его ближайший сосед с расстоянием .

При разработке алгоритмов классификации будем предполагать следующие условия:

• системе предъявляется 1 или более образов для классификации;

• эталоны классов задаются списком своих номеров в общем списке образов.

Оба правила классификации могут быть реализованы одним алгоритмом отыскания минимума расстояния. Для правила ближайшего эталона будет просматриваться упомянутый список номеров эталонов, а для правила ближайшего соседа – базовая выборка с известной классификацией. Если вычислена матрица расстояний между парами образов, то ближайший сосед (или ближайший эталон) может быть определен с ее помощью. В противном случае будет проведен расчет расстояний между просматриваемыми и классифицируемыми образами.

Для обеих рассматриваемых архитектур будет применяться следующая идея: первоначально отыскиваются частичные решения задачи поиска ближайшего соседа (или ближайшего эталона) многими потоками, а затем производится редукция частичных решений.

Пусть просмотр базовой выборки с известной кластеризацией (для правила ближайшего соседа) или списка номеров эталонов (для правила ближайшего эталона) равномерно распределен между потоками. Для каждого просматриваемого образа перебираются все предъявляемые образы. Таким образом, будут найдены частичные решения задачи поиска ближайшего соседа (или ближайшего эталона) для каждого предъявляемого образа в количестве, равном количеству потоков. Затем производится редукция частичных результатов для каждого предъявленного образа и причисление его к соответствующему классу.

Ниже приведен код реализации алгоритма классификации по обоим правилам на C++. В качестве средства распараллеливания используется OpenMP 2.0. Для распараллеливания итераций цикла достаточно указания директивы #pragma omp parallel for. Избежать конфликтного доступа нескольких потоков к одной ячейке памяти позволяет директива OpenMP 2.0 #pragma omp atomic.

//Выделить память для частичных (вытекающих из анализа фрагмента базовой выборки) решений задачи нахождения номера ближайшего соседа (или ближайшего эталона) для каждого предъявленного образа.

#ifdef _OPENMP

int threadsCnt = omp_get_max_threads();

#else

int threadsCnt = 1;

#endif

float* min_val = new float[threadsCnt * imagesFeatures.size()];

int* min_idx = new int[threadsCnt * imagesFeatures.size()];

#ifdef _OPENMP

#pragma omp parallel

#endif

//Вход в параллельный регион, создание потоков.

{

//Получить номера потоков и проинициализировать частичные решения.

#ifdef _OPENMP

int threadIdx = omp_get_thread_num();

#else

int threadIdx = 0;

#endif

for (int l = 0; l < (int)(imagesFeatures.size()); l++)

{

min_val[threadIdx + threadsCnt * l] = FLT_MAX;

min_idx[threadIdx + threadsCnt * l] = -1;

}

//Найти частичные решения для каждого предъявленного образа.

#ifdef _OPENMP

#pragma omp for

#endif

for (int i = 0; i < ((baseSize > 0) ? baseSize : fgaugeslen); i++)

{

int j = (baseSize > 0) ? i : fgauges[i];

for (vector< vector<float> >::iterator imageFeatures = imagesFeatures.begin(); imageFeatures != imagesFeatures.end(); imageFeatures++)

{

int l = (int)(imageFeatures - imagesFeatures.begin());

float dstVal = 0;

if (DistMatr == NULL)

{

if (fmasklen > 0)

{

for (int k = 0; k < fmasklen; k++)

{

float delta = (*imageFeatures)[fmask[k]] - x[j + fmask[k] * MAX_IMAGE_CNT];

dstVal += delta * delta;

}

}

else

{

for (int k = 0; k < n; k++)

{

float delta = (*imageFeatures)[k] - x[j + k * MAX_IMAGE_CNT];

dstVal += delta * delta;

}

}

}

else

{

dstVal = DistMatr[startImg + l + j * MAX_IMAGE_CNT];

}

if (dstVal < min_val[threadIdx + threadsCnt * l])

{

min_val[threadIdx + threadsCnt * l] = dstVal;

min_idx[threadIdx + threadsCnt * l] = j;

}

}

}

//Для каждого предъявленного образа: получение номера ближайшего соседа (или ближайшего эталона) и причисление к соответствующему классу.

#ifdef _OPENMP

//Редукция частичных решений.

#pragma omp for

for (int l = 0; l < (int)(imagesFeatures.size()); l++)

{

int min_idx_scalar = -1;

float min_val_scalar = FLT_MAX;

for (int i = 0; i < threadsCnt; i++)

{

if (min_val[i + threadsCnt * l] < min_val_scalar)

{

min_val_scalar = min_val[i + threadsCnt * l];

min_idx_scalar = min_idx[i + threadsCnt * l];

}

}

int imageClass = clsInfo[min_idx_scalar];

clsInfo[startImg + l] = imageClass;

//Атомарное увеличение размера класса (образы необязательно распределяются по различным классам).

#pragma omp atomic

clsDims[imageClass]++;

}

#else

//Однопоточный случай – редукция не требуется.

for (int l = 0; l < (int)(imagesFeatures.size()); l++)

{

int imageClass = clsInfo[min_idx[l]];

clsInfo[startImg + l] = imageClass;

clsDims[imageClass]++;

}

#endif

}

delete[] min_val;

delete[] min_idx;

Пусть каждый мультипроцессор будет выполнять классификацию одного из предъявленных образов. В общей памяти будут вычислены частичные решения задачи поиска ближайшего соседа (или ближайшего эталона) каждым потоком в блоке. Затем нулевой поток произведет редукцию частичных решений и причислит образ к соответствующему классу. Поскольку графические процессоры NVidia не обладают возможностями выполнения атомарных операций, мы имеем право увеличить размер класса только в том случае, когда предъявлен всего один образ для классификации. В противном случае размеры классов придется обновить отдельно.

Ниже приведен код реализации алгоритма классификации по обоим правилам на CUDA C.

__constant__ int fgauges[gaugeslimit];

template <bool usemask, bool usegauges, bool ignorematr> __global__ void g_classify(float* DistMatr, unsigned char* clsInfo, unsigned short* clsDims, float* x, int fgaugeslen, int fmasklen, int startImg)

{

extern __shared__ float newX[];

int* index = (int*)(newX + fmasklen);

float* D = (float*)(index + blockDim.x);

if (ignorematr)

{

//Скопировать атрибуты предъявленного образа в общую память.

for (int i = threadIdx.x; i < fmasklen; i += blockDim.x)

newX[i] = x[startImg + blockIdx.x + (usemask ? fmask[i] : i) * MAX_IMAGE_CNT];

__syncthreads();

}

//Инициализация частичных решений.

float min_val = FLT_MAX;

int min_idx = -1;

//Определение частичных решений.

for (int i = threadIdx.x; i < fgaugeslen; i += blockDim.x)

{

int j = (usegauges ? fgauges[i] : i);

float dstVal = 0;

if (ignorematr)

{

for (int k = 0; k < fmasklen; k++)

{

float delta = newX[k] - x[j + (usemask ? fmask[k] : k) * MAX_IMAGE_CNT];

dstVal += delta * delta;

}

}

else

{

dstVal = DistMatr[j + (startImg + blockIdx.x) * MAX_IMAGE_CNT];

}

bool better_fact = (dstVal < min_val);

float better_val = min(min_val, dstVal);

min_val -= min_val * (int)(better_fact);

min_val += better_val * (int)(better_fact);

min_idx = min_idx + (j - min_idx) * (int)(better_fact);

}

index[threadIdx.x] = min_idx;

D[threadIdx.x] = min_val;

__syncthreads();

//Редукция частичных решений.

if (threadIdx.x == 0)

{

float min_val = FLT_MAX;

int min_idx = -1;

for (int i = 0; i < blockDim.x; i++)

{

bool better_fact = (D[i] < min_val);

float better_val = min(min_val, D[i]);

min_val -= min_val * (int)(better_fact);

min_val += better_val * (int)(better_fact);

min_idx = min_idx + (i - min_idx) * (int)(better_fact);

}

min_idx = index[min_idx];

int imageClass = clsInfo[min_idx];

clsInfo[startImg + blockIdx.x] = (unsigned char)(imageClass);

if (gridDim.x == 1) clsDims[imageClass]++;

}

}

Отдельная операция обновления размеров классов в случае, когда предъявляется более одного образа для классификации. Размеры классов увеличиваются в соответствии с проведенной классификацией. Список принадлежности образов классам просматривается начиная с первого образа, предъявленного для классификации. Этот «хвост» равномерно распределен между потоками в блоке. Каждый поток вычислит свою частичную добавку к размеру каждого класса. А затем потоков произведут редукцию частичных добавок к размерам соответствующих классов. Ниже приведен код реализации алгоритма обновления размеров классов на CUDA C.

__global__ void g_upd_clsDims(unsigned char* clsInfo, unsigned short* clsDims, int c, int N, int startImg)

{

extern __shared__ int shared[];

//Инициализировать частичные суммы (размеры классов).

for (int i = 0; i < c; i++)

shared[i + threadIdx.x * c] = 0;

//Определить частичные суммы (размеры классов).

for (int i = startImg + threadIdx.x; i < N; i += blockDim.x)

shared[clsInfo[i] + threadIdx.x * c]++;

__syncthreads();

//Выполнить редукцию частичных сумм (размеров) для каждого класса.

if (threadIdx.x < c) //Здесь мы полагаем blockDim.x >= c.

{

float add = 0;

for (int i = 0; i < blockDim.x; i++)

add += shared[threadIdx.x + i * c];

clsDims[threadIdx.x] += add;

}

}

Заметим, что каждый поток с номером меньшим выполнит сложений, где первое слагаемое – количество сложений при вычислении промежуточного результата, а второе слагаемое – количество сложений при редукции. Время работы вычислительного ядра можно оптимизировать, уменьшив это количество сложений. Нетрудно убедиться в том, что минимума можно достичь при количестве потоков в блоке, равном .