3.1.1 Решение СЛАУ на графических ускорителях

В CUDA, названия данных команд заимствованы со ставшего стандартом библиотеки LAPACK. LAPACK - это библиотека на языке Fortran 90 содержащая набор подпрограмм для решения наиболее часто встречающихся задач линейной алгебры. Она разрабатывалась с учетом свойств большинства современных архитектур для достижения наивысшей производительности.

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

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

Таблица 3.1.1 - Обозначение типа подоперации в LAPACK.

S	REAL
D	DOUBLE PRECISION
C	COMPLEX
Z	COMPLEX*16

Буквы XX обозначают тип матрицы с которой работает подпрограмма:

Таблица 3.1.2 - Обозначение типа матрицы в LAPACK.

BD	двухдиагональная
DI	диагональная
GE	произвольная матрица
GT	трехдиагональная
PO	симметричная, положительно определенная
SY	симметричная
TG	треугольная

Буквы ZZZ обозначают тип проводимых вычислений, например TRF обозначает LU разложение, EV - поиск собственных значений. Более подробный перечень обозначений можно найти в справке LAPACK [14].

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

1. gerf.

2. ormqr.

3. trsm.

Рассмотрим каждый шаг.

Geqrf. Данная команда выполняет QR разложение матриц, возвращая верхнюю треугольную матрицу R в верхней треугольной части матрицы A и матрицу Q в форме Хаусхолдера, размещенной в нижней треугольной части матрицы A.

Код QR разложения матрицы на языке CUDA:

cusolveSafeCall(cusolverDnDgeqrf(solver_handle, Nrows, Ncols, d_A, Nrows, d_TAU, work, work_size, devInfo));

int devInfo_h = 0;

gpuErrchk(cudaMemcpy(&devInfo_h, devInfo, sizeof(int), cudaMemcpyDeviceToHost));

if (devInfo_h != 0) std::cout << "Unsuccessful gerf execution\n\n";

Листинг 3.1.1 - QR разложение матрицы на языке CUDA.

Интерфейс подпрограммы cusolverDnDgeqrf:



Таблица 3.1.3 - Интерфейс подпрограммы cusolverDnDgeqrf.

Параметр	Память	In/out	Описание
handle	host	input	Указатель на контекст библиотеки cuSolverDN.
m	host	input	Число строк матрицы A.
n	host	input	Число столбцов матрицы A.
A	device	in/out	<type> массив размерности lda * n, где lda не меньше, чем max(1,m).
lda	host	input	Ведущая размерность двумерного массива используемой для записи матрицы A.
TAU	device	output	<type> массив с размерностью минимум min(m,n).
Workspace	device	in/out	Рабочее пространство, <type> массив с размером L.
Lwork	host	input	Размер рабочего массива Workspace.
devInfo	device	output	if info = 0, LU разложение прошло успешно. if info = -i, i-й параметр неверен.

После выполнения данного кода верхняя треугольная часть матрицы A содержит элементы R. Также возвращает матрицу Q в форме Хаусхолдера в нижней треугольной части матрицы A, в то время как коэффициент масштабирования возвращается через параметр TAU.

Ormqr возвращает произведение матриц A и Q.

Код произведение матриц на языке CUDA:

cusolveSafeCall(cusolverDnDormqr(solver_handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_N, Nrows, Ncols, min(Nrows, Ncols), d_A, Nrows, d_TAU, d_Q, Nrows, work, work_size, devInfo));

Листинг 3.1.2 - произведение матриц на языке CUDA.

Интерфейс подпрограммы cusolverDnDormqr:

Таблица 3.1.4 - Интерфейс подпрограммы cusolverDnDormqr.

Параметр	Память	In/out	Описание
handle	host	input	Указатель на контекст библиотеки cuSolverDN.
side	host	input	Обозначает, находится ли матрица Q левой или правой стороне C.
trans	host	input	Операция op(Q), которая не- или (conj.) транспонированная.
m	host	input	Число строк матрицы A.
n	host	input	Число столбцов матрицы A.
k	host	input	Число отражений элементов.
A	device	in/out	<type> массив размерности lda * n, где lda не меньше, чем max(1,m).
lda	host	input	Ведущая размерность двумерного массива используемой для записи матрицы A. if side == CUBLAS_SIDE_LEFT, lda >= max(1,m); if side == CUBLAS_SIDE_RIGHT, lda >= max(1,n).
tau	device	output	<type> массив с размерностью минимум min(m,n).
C	device	in/out	<type> массив с размером ldc * n. При завершении, C перезаписывается операцией op(Q)*C.
ldc	host	input	Ведущая размерность двумерного массива используемой для записи матрицы C. ldc >= max(1,m).
work	device	in/out	Рабочее пространство, <type> массив с размером L.
lwork	host	input	Размер рабочего массива Workspace.
devInfo	device	output	if info = 0, LU выполнение прошло успешно. if info = -i, i-й параметр неверен.

Trsm. Еще одной часто используемой BLAS-функцией 3-го уровня является TRSM(s,ul,tz,d,m,n,б,Z,lz,Y,ly). Она предназначена для решения левостороннего (ZЧX=б•Y) или правостороннего (XЧZ=б•Y) матричного уравнения для всевозможных вариантов представления треугольной матрицы Z. Установкой параметров можно задать не только, с какой стороны (s) в этом уравнении участвует матрица Z, но и уточнить, является она нижней или верхней треугольной (ul), транспонированная ли она (tz) или сопряж?нная (для комплексных), единичная ли у не? диагональ (d) или нет. Функция может применяться не только к матрицам целиком, но и к отдельным матричным блокам, поэтому при е? вызове задаются также параметры (lz,ly), уточняющие размеры строк используемых матриц.

const double alpha = 1.;

cublasSafeCall(cublasDtrsm(cublas_handle, CUBLAS_SIDE_LEFT, CUBLAS_FILL_MODE_UPPER, CUBLAS_OP_N, CUBLAS_DIAG_NON_UNIT, Ncols, Ncols,

&alpha, d_R, Ncols, d_B, Ncols));

Листинг 3.1.3 - Решение верхне-треугольной линейной системы на языке CUDA.

Таблица 3.1.5 - Интерфейс подпрограммы cublasDtrsm.

Параметр	Память	In/out	Описание
handle		input	Указатель на контекст cuBLAS.
side		input	Обозначает с какой стороны А от Х.
uplo		input	Обозначает, верхняя или няжняя часть матрицы А записана.
trans		input	Операция op(A), которая не- или (conj.) транспонирована.
diag		input	Обозначает, являются ли элементы главной диагонали матрицы 1.
m		input	Число строк В.
n		input	Число столбцов В.
alpha	host or device	input	<type> множитель
A	device	input	<type> массив размерности lda x m, где lda>=max(1,m) if side == CUBLAS_SIDE_LEFT и lda x n, где lda>=max(1,n) иначе.
lda		input	Ведущая размерность двумерного массива, в котором записана матрица А
B	device	in/out	<type> массив. Он имеет размерность ldb x n, где ldb>=max(1,m).
ldb		input	Ведущая размерность двумерного массива, в котором записана матрица В

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

3.1.2 Создание веб-службы с помощью gSOAP

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

Процесс формирования кода сервера в gSOAP почти полностью автоматизирован. Для него разработан свой язык описания сервисов, похожий по грамматике на C и C++. Благодаря чему, для разработки веб-службы не требуется знание структуры файла WSDL.

Разработка веб-службы с помощью gSOAP начинается с написания заголовочного файла с описанием функции, реализуемой на сервере:

//gsoap ns service name: spline

//gsoap ns service style: rpc

//gsoap ns service encoding: encoded

//gsoap ns service namespace: urn:spline

//gsoap ns service location: http://localhost:9999

//gsoap ns service method-action: spline urn:spline

struct Vector {

double *__ptr;

int __size;

};

int ns__BuildSpline(int size, struct Vector X, struct Vector Y, struct Vector F, struct Vector* result);

Листинг 3.1.4 - Содержание заголовочного файла CudaSpline.h.

С помощью команды soapcpp2.exe -2 -i -e CudaSpline.h создаем несколько файлов с исходным кодом, который нам необходимо лишь дополнить. Параметры команды значат следующее:

· -2 - использовать стандарт SOAP v1.2;

· -i - сгенерировать код классов, наследуемых от структуры soap;

· -e - использовать стиль SOAP RPC при определении типов.

· В результате данной команды создаются следующие файлы:

· spline.wsdl - файл с описанием сервиса в формате WSDL.

· t.xsd - пользовательские типы хранятся в этом файле.

· spline.req.xml и spline.res.xml - примеры запроса и ответа веб-службы.

· splineProxy.cpp и splineProxy.h - код для клиентской стороны. Поскольку мы используем MATLAB в качестве клиента, эти файлы нам не нужны.

· splineService.cpp и spline.h - код для серверной части.

· spline.nsmap - пространства имен используемые сторонами.

· soapStub.h - описание структур, используемых на сервере и клиенте.

Благодаря данному интерфейсу, к веб-службе можно обращаться из любого клиента поддерживающего интерфейс SOAP. На веб-службе на вычислительном узле с графическим ускорителем реализована единственная функция, которая находит коэффициенты сплайна: BuildSpline(int size, struct Vector X, struct Vector Y, struct Vector F, struct Vector *result). На вход данной функции подается массив точек. Функция возвращает массив коэффициентов сплайна.

Рисунок 12 - Схема WSDL описания веб-службы.



Рисунок 13 - Состав проекта веб-службы с использованием технологии CUDA.

Нам необходимо переопределить виртуальный метод класса spline: int spline::BuildSpline(int size, struct Vector X, struct Vector Y, struct Vector F, struct Vector *result). Согласно правилам gSOAP, данный метод должен возвращать значение через последний параметр. Сам метод должен возвращать значение SOAP_OK, если выполнение прошло без ошибок. Данная константа определена в файле stdsoap2.h.

3.2 Клиент на MATLAB

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

Данная среда содержит функции, которые позволяют обращаться напрямую к веб-службам SOAP:

· callSoapService - посылает SOAP сообщение сервису;

· createSoapMessage - создает SOAP сообщение;

function result = Utilits(X, Y, F)

values = { size(X, 2), X, Y, F };

names = { 'size', 'X', 'Y', 'F' };

types = { ...

'{http://www.w3.org/2001/XMLSchema}int', ...

'{http://www.w3.org/2001/XMLSchema}double[]', ...

'{http://www.w3.org/2001/XMLSchema}double[]', ...

'{http://www.w3.org/2001/XMLSchema}double[]' ...

};

soapMessage = createSoapMessage( 'urn:spline', 'BuildSpline', values, names, types, 'rpc');

response = callSoapService( 'http://localhost:9999', '', soapMessage); temp = parseSoapResponse(response); result = temp.item;

Листинг 3.2.1 - Скрипт MATLAB для обращения к нашей веб-службе.

В ходе выполнения данной работы выяснилась неполная совместимость этих функций с реализацией веб-службы с помощью gSOAP. При передаче массивов на веб-службу, не указывалась длина массива в запросе. В связи с этим была сделана корректировка вышеупомянутых стандартных функций для правильной работы с веб-службой, созданной с помощью gSOAP. Была изменена функция из стандартного набора MATLAB: функция populate(dom,node,values,names,types). Данная функция участвует в формировании XML документа, отправляемого в запросе на веб-службу. В стандартной реализации данной функции в описании массивов добавляются XML-атрибуты xsi:type="SOAP-ENC:Array" SOAP-ENC:arrayType="xsd:double[]". Чтобы веб-служба на gSOAP обрабатывала данный запрос, требуется добавлять число элементов массива в XML-атрибут SOAP-ENC:arrayType. Например: SOAP-ENC:arrayType ="xsd:double[2]"



4. Тестирование

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

(4.1)

4.1 Эллиптический параболоид

Эллиптический параболоид -- поверхность, задаваемая функцией вида:

(4.1.1)

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

Если a=b, то эллиптический параболоид представляет собой поверхность вращения, образованную вращением параболы вокруг её оси симметрии.

Построим несколько сплайн-поверхностей на основе эллиптического параболоида и измерим отклонение для каждого из них. a = b = 1.



Рисунок 4.1.1 - Сплайн-поверхность, построенная на 20 точках эллиптического параболоида

Цветом обозначена разница между сплайном и исходной поверхностью. Число узлов интерполяции - 20.

При N = 20, мера S = 1,2763.

Рисунок 4.1.2 - Сплайн-поверхность, построенная на 40 точках эллиптического параболоида



Цветом обозначена разница между сплайном и исходной поверхностью. Число узлов интерполяции - 40.

При N = 40, мера S = 0,6546.

Рисунок 4.1.3 - Сплайн-поверхность, построенная на 60 точках эллиптического параболоида

Цветом обозначена разница между сплайном и исходной поверхностью. Число узлов интерполяции - 60.

При N = 60, мера S = 0,1801.

Зависимость меры S от N изображено на рисунке 4.1.4.



Рисунок 4.1.4 - зависимость меры S от N для эллиптического параболоида.

4.2 Трансцендентная поверхность 1

Данная поверхность задается уравнением:

z = sin(x) + cos(y) (4.2.1).

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



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

При N = 20, мера S = 0,6383.

Рисунок 4.2.2 - Сплайн-поверхность с числом узлов интерполяции - 40. Цветом обозначена разница между сплайном и исходной поверхностью.

При N = 40, мера S = 0,1251.



Рисунок 4.2.3 - Сплайн-поверхность с числом узлов интерполяции - 60. Цветом обозначена разница между сплайном и исходной поверхностью.

При N = 60, мера S = 0,0642.

График зависимости S от N изображен на рисунке 4.2.4.

Рисунок 4.2.4 - график зависимости S от N для первой трансцендентной поверхности.



4.3 Трансцендентная поверхность 2

В качестве второй трансцендентной поверхности возьмем поверхность, задающейся функцией:

z = sin(x) /(x² + y² + 0.3) (4.3.1).

Рисунок 4.3.1 - Сплайн-поверхность с числом узлов интерполяции - 20. Цветом обозначена разница между сплайном и исходной поверхностью.

При N = 20, мера S = 0,2043.



Рисунок 4.3.2 - Сплайн-поверхность с числом узлов интерполяции - 40. Цветом обозначена разница между сплайном и исходной поверхностью.

При N = 40, мера S = 0,0616.

Рисунок 4.3.3 - Сплайн-поверхность с числом узлов интерполяции - 60. Цветом обозначена разница между сплайном и исходной поверхностью.



При N = 60, мера S = 0,0818.

Рисунок 4.3.4 - график зависимости S от N для второй трансцендентной поверхности.

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



Заключение

В ходе данной работы был разработан метод построения интерполяционного сплайна на нерегулярной сетке с использованием технологии CUDA. Также, в частности, был разработан алгоритм решения системы линейных алгебраических уравнений на графических ускорителях с поддержкой CUDA для нахождения коэффициентов сплайн-поверхности Смоляка с помощью библиотек cuBLAS, cuSPARSE и cuSOLVER.

Разработан веб-сервис для предоставления доступа клиентам, у которых нет необходимых графических ускорителей, к ресурсам вычислительного узла с требуемой аппаратной поддержкой, посредством открытого протокола прикладного уровня SOAP.

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

Сделана визуализация сплайн-поверхностей на языке MATLAB с цветовым кодированием отклонения сплайна от оригинальной поверхности.

Тестирование показало ожидаемое снижение отклонения сплайна с увеличением числа узлов интерполяции. Также тестирование подтвердило корректность реализации алгоритма нахождения коэффициентов сплайна Смоляка, и показало согласованное взаимодействие сетевых узлов по SOAP и модулей системы в целом.

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

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

Созданная библиотека построения сплайн-поверхности Смоляка является органичным дополнением к существующим библиотекам CUDA и может использоваться сообществом разработчиков CUDA для построения сплайнов на нерегулярной интерполяционной сетке. Отсутствие зависимости от сторонних библиотек и использование только стандартных библиотек самой CUDA (cuBLAS, cuSPARSE, cuSOLVER) способствует быстрому и легкому применению данной разработки.

Интерфейс SOAP позволяет использовать данную веб-службу с любым клиентом, поддерживающим данный стандарт, без изменений кода.

Использование инструментария gSOAP позволяет без труда вносить изменения в существующий веб-сервис не прибегая к программированию сетевого кода. Компилятор "soapcpp2" автоматически обновляет схему типов данных и функций XML веб-службы (WSDL) согласно коду программы на C/C++.



Список литературы

1. Козлов В. А. Открытые информационные системы. М.: ФиС, 1999.

2. Д. Каханер, К.Моулер, С.Нэш, Численные методы и программное обеспечение М.: Мир, 1998. 575 с

3. Калиткин Н.Н. Численные методы. / Калиткин Н.Н. - М.: Наука, 1978. - 512с.

4. Турчак Л.И. Основы численных методов. / Турчак Л.И. - М.: Наука, 1987. - 320с.

5. Самарский А.А. Численные методы / Самарский А.А., Гулин А.В. - M.: Наука, 1989. - 431c.

6. Стечкин С.Б. Сплайны в вычислительной математике. / Стечкин С.Б., Субботин Ю.Н. - М.: Наука, 1976. - 349с.

7. Шуп. Т.Е. Прикладные численные методы в физике и технике. / Т.Е.Шуп - М.: Высш. шк., 1990. - 255с.

8. Смоляк С.А. Оптимальное восстановление функций и связанные с ним геометрические характеристики множеств // Труды третьей зимней школы по математическому программированию и смежным вопросам (24 января - 3 февраля 1970 г.,г.Дрогобыч), вып.III. - М.: 1970. - С. 509-557.

9. Смоляк С.А. Сплайны и их применение // Экономика и математические методы. - 1971. - Т.7, вып.3. - . 419-431.

10. Harder R.L., Desmarais R.N. Interpolation using surface splines // Journal of Aircraft. - 1972. - Vol.9, № 2. - P. 189-191.

11. Ашкеназы В. О. А98 Сплайн-поверхности: Основы теории и вычислительные алгоритмы: Учебное пособие. - Тверь: Тверской гос. ун-т, 2003. - 82 с.

12. MATLAB - MathWorks [Электронный ресурс]. URL: http://www.mathworks.com/products/matlab/?s_tid=hp_fp_ml (дата обращения: 12.05.2016).

13. System Requirements & Platform Availability [Электронный ресурс]. URL: http://www.mathworks.com/products/availability/index.html#ML (дата обращения: 12.05.2016).

14. LAPACK Routines [Электронный ресурс]. URL: http://www.netlib.org/lapack/lawn41/node111.html (дата обращения: 12.05.2016).

15. 

ПРИЛОЖЕНИЕ А

Код веб-службы.

/* Includes, system */

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

/* Includes, cuda */

#include <cuda_runtime.h>

#include <cublas_v2.h>

#include <helper_cuda.h>

#include "LinSolver.cuh"

#include "spline.nsmap"

#include "soapsplineService.h"

int main(int argc, char **argv)

{

const int SERVICE_PORT = 9999;

printf_s("Started");

std::auto_ptr<splineService> srv(new splineService());

if (srv->run(SERVICE_PORT))

{

srv->soap_stream_fault(std::cerr);

exit(-1);

}

return 0;

}

int splineService::BuildSpline(int size, struct Vector X, struct Vector Y, struct Vector F, struct Vector *result) SOAP_PURE_VIRTUAL

{

int size_ext = size + 3;

double* A = new double[size_ext * size_ext];

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

{

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

{

A[i * size_ext + j] = 0.0;

if (i != j)

{

A[i* size_ext + j] = (X.__ptr[i] - X.__ptr[j]) * (X.__ptr[i] - X.__ptr[j]) + (Y.__ptr[i] - Y.__ptr[j]) * (Y.__ptr[i] - Y.__ptr[j]);

A[i* size_ext + j] = A[i* size_ext + j] * log(A[i* size_ext + j]);

}

}

}

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

{

A[i * size_ext + size + 0] = 1.0;

A[i * size_ext + size + 1] = X.__ptr[i];

A[i * size_ext + size + 2] = Y.__ptr[i];

A[(size + 0) * size_ext + i] = 1.0;

A[(size + 1) * size_ext + i] = X.__ptr[i];

A[(size + 2) * size_ext + i] = Y.__ptr[i];

}

for (int i = size; i < size + 3; ++i)

{

for (int j = size; j < size + 3; ++j)

{

A[i * size_ext + j] = 0.0;

}

}

double* C = new double[size_ext * size_ext];

LinSolver(A, size_ext, size_ext, F.__ptr, C);

//result = new Vector();

result->__ptr = new double[size_ext];

result->__size = size_ext;

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

result->__ptr[i] = C[i];

free(A);

free(C);

return SOAP_OK;

}



ПРИЛОЖЕНИЕ Б

Код вычисления СЛАУ с использованием CUDA.

#include "cuda_runtime.h"

#include "device_launch_paraMeters.h"

#include <iostream>

#include <iomanip>

#include <stdlib.h>

#include <stdio.h>

#include <assert.h>

#include <cusolverDn.h>

#include <cublas_v2.h>

#include <cuda_runtime_api.h>

#include "Utilities.cuh"

#define BLOCK_SIZE 32

__global__ void copy_kernel(const double * __restrict d_in1, double * __restrict d_out1, const double * __restrict d_in2, double * __restrict d_out2, const int M, const int N)

{

const int i = blockIdx.x * blockDim.x + threadIdx.x;

const int j = blockIdx.y * blockDim.y + threadIdx.y;

if ((i < N) && (j < N))

{

d_out1[j * N + i] = d_in1[j * M + i];

d_out2[j * N + i] = d_in2[j * M + i];

}

}

void LinSolver(double *h_A, int Nrows, int Ncols, double *h_C, double *h_B)

{

int work_size = 0;

int *devInfo;

gpuErrchk(cudaMalloc(&devInfo, sizeof(int)));

cusolverDnHandle_t solver_handle;

cusolverDnCreate(&solver_handle);

cublasHandle_t cublas_handle;

cublasSafeCall(cublasCreate(&cublas_handle));

double *d_A;

gpuErrchk(cudaMalloc(&d_A, Nrows * Ncols * sizeof(double)));

gpuErrchk(cudaMemcpy(d_A, h_A, Nrows * Ncols * sizeof(double), cudaMemcpyHostToDevice));

double *d_TAU;

gpuErrchk(cudaMalloc((void**)&d_TAU, min(Nrows, Ncols) * sizeof(double)));

cusolveSafeCall(cusolverDnDgeqrf_bufferSize(solver_handle, Nrows, Ncols, d_A, Nrows, &work_size));

double *work; gpuErrchk(cudaMalloc(&work, work_size * sizeof(double)));

cusolveSafeCall(cusolverDnDgeqrf(solver_handle, Nrows, Ncols, d_A, Nrows, d_TAU, work, work_size, devInfo));

int devInfo_h = 0;

gpuErrchk(cudaMemcpy(&devInfo_h, devInfo, sizeof(int), cudaMemcpyDeviceToHost));

if (devInfo_h != 0)

std::cout << "Unsuccessful gerf execution\n\n";

double *h_Q = (double *)malloc(Nrows * Nrows * sizeof(double));

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

{

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

{

if (j == i) h_Q[j + i*Nrows] = 1.;

else h_Q[j + i*Nrows] = 0.;

}

}

double *d_Q;

gpuErrchk(cudaMalloc(&d_Q, Nrows * Nrows * sizeof(double)));

gpuErrchk(cudaMemcpy(d_Q, h_Q, Nrows * Nrows * sizeof(double), cudaMemcpyHostToDevice));

cusolveSafeCall(cusolverDnDormqr(solver_handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_N, Nrows, Ncols, min(Nrows, Ncols), d_A, Nrows, d_TAU, d_Q, Nrows, work, work_size, devInfo));

double *d_C;

gpuErrchk(cudaMalloc(&d_C, Nrows * Nrows * sizeof(double)));

gpuErrchk(cudaMemcpy(d_C, h_C, Nrows * Nrows * sizeof(double), cudaMemcpyHostToDevice));

cusolveSafeCall(cusolverDnDormqr(solver_handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_T, Nrows, Ncols, min(Nrows, Ncols), d_A, Nrows, d_TAU, d_C, Nrows, work, work_size, devInfo));

double *d_R;

gpuErrchk(cudaMalloc(&d_R, Ncols * Ncols * sizeof(double)));

double *d_B;

gpuErrchk(cudaMalloc(&d_B, Ncols * Ncols * sizeof(double)));

dim3 Grid(iDivUp(Ncols, BLOCK_SIZE), iDivUp(Ncols, BLOCK_SIZE));

dim3 Block(BLOCK_SIZE, BLOCK_SIZE);

copy_kernel << <Grid, Block >> >(d_A, d_R, d_C, d_B, Nrows, Ncols);

const double alpha = 1.;

cublasSafeCall(cublasDtrsm(cublas_handle, CUBLAS_SIDE_LEFT, CUBLAS_FILL_MODE_UPPER, CUBLAS_OP_N, CUBLAS_DIAG_NON_UNIT, Ncols, Ncols,

&alpha, d_R, Ncols, d_B, Ncols));

gpuErrchk(cudaMemcpy(h_B, d_B, Ncols * Ncols * sizeof(double), cudaMemcpyDeviceToHost));

cusolverDnDestroy(solver_handle);

}



ПРИЛОЖЕНИЕ В

Код клиента на MATLAB.

clear;

colormap gray;

% 20 40 60

n = 60;

for i = 1 : n

theta = (rand() - 0.5) * 6;

phi = (rand() - 0.5) * 6;

X(i) = theta;

Y(i) = phi;

F(i) = SurfaceEq(theta, phi);

end

scatter3(X, Y, F, 40, [1, 0, 0], 'filled');

hold on;

F(size(X, 2) + 1) = 0.0;

F(size(X, 2) + 2) = 0.0;

F(size(X, 2) + 3) = 0.0;

C = Utilits(X, Y, F);

% [mX, mY] = meshgrid(min(X) : 0.5 : max(X), min(Y) : 0.5 : max(Y));

[mX, mY] = meshgrid(-3 : 0.2 : 3, -3 : 0.2 : 3);

S(size(mX, 1), size(mX, 2)) = 0.0;

for i = 1 : size(X, 2)

Q = (mX - X(i)) .^ 2 + (mY - Y(i)) .^ 2;

S = S + C(i) .* Q .* log(Q);

end

S = S + C(size(X, 2) + 1) + C(size(X, 2) + 2) .* mX + C(size(X, 2) + 3) .* mY;

for i = 1 : size(mX, 1)

for j = 1 : size(mX, 2)

O(i, j) = SurfaceEq(mX(i, j), mY(i, j));

end

end

mesh(mX, mY, S, abs(S - O));

alpha(0.1);

Diff = S - O;

r = sqrt(mean(Diff(:).^2));

disp(r);

disp(n);

Размещено на Allbest.ru