Определение детерминанта с помощью CUDA [закрыто]

Существует ли какая-либо библиотека или свободно доступный код, который вычислит определитель small (6x6) матрица двойной точности целиком на GPU?

4 голоса | спросил Paul Caheny 2 PM00000050000001731 2012, 17:57:17

3 ответа


0

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

Я не собираюсь писать реальный код, но это должно быть полезно.

1) Launch # blocks = # matrices. Каждый блок вычисляет определитель каждой матрицы.

2) det (A) = det (A11 * A22 - A21 * A12); где A - 6x6, A11, A12, A21, A22 - подматрицы 3x3 матрицы A.

3) Напишите функцию device , которая умножает матрицы для матриц 3x3

4) det матрицы 3x3 легко вычислить: используйте формулу отсюда .

РЕДАКТИРОВАТЬ : очевидно (2) работает, только если A21 * A12 == A12 * A21

Альтернативой может быть следующее

1) Декомпозиция LU путем исключения Гаусса для каждого 6x6 матрица

2) Умножьте диагональные элементы U, чтобы получить определитель.

ответил Pavan Yalamanchili 2 PM000000100000004631 2012, 22:18:46
0

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

Я думаю, что проблема вычисления определителя матриц является per se интересной проблемой, которая может встречаться в приложениях несколько раз. В настоящее время я не знаю ни одной библиотеки, предоставляющей процедуры для вычисления определителя с помощью CUDA (ни cuBLAS, ни cuSOLVER есть такая подпрограмма), поэтому у вас есть две возможности:

  1. Внедрение вашего собственного подхода, как указал Паван;
  2. Управляйте другими доступными подпрограммами.

Что касается последнего пункта, одной из возможностей является использование разложения Холецкого , а затем вычисление определитель как квадрат произведения диагональных элементов матрицы Холецкого. В Matlab это будет выглядеть так:

prod(diag(chol(A)))^2 

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

Приведенный ниже код предназначен для больших матриц, поэтому он сразу пригодится людям, которым необходимо вычислить определитель больших матриц. Но как адаптировать его к нескольким маленьким матрицам? Одна из возможностей - использовать потоки cuSOLVER для разложения Холецкого, а затем использовать новую функцию динамического параллелизма в Thurst 1.8. Обратите внимание, что в CUDA 7.0 #include "cuda_runtime.h" #include "device_launch_paraMeters.h" #include<iostream> #include<iomanip> #include<stdlib.h> #include<stdio.h> #include<assert.h> #include<ostream> #include <cusolverDn.h> #include <cublas_v2.h> #include <cuda_runtime_api.h> #include "Utilities.cuh" #include <thrust/iterator/counting_iterator.h> #include <thrust/iterator/transform_iterator.h> #include <thrust/iterator/permutation_iterator.h> #include <thrust/functional.h> #include <thrust/fill.h> #include <thrust/device_vector.h> #include <thrust/host_vector.h> #include <thrust/copy.h> /*************************/ /* STRIDED RANGE FUNCTOR */ /*************************/ template <typename Iterator> class strided_range { public: typedef typename thrust::iterator_difference<Iterator>::type difference_type; struct stride_functor : public thrust::unary_function<difference_type,difference_type> { difference_type stride; stride_functor(difference_type stride) : stride(stride) {} __host__ __device__ difference_type operator()(const difference_type& i) const { return stride * i; } }; typedef typename thrust::counting_iterator<difference_type> CountingIterator; typedef typename thrust::transform_iterator<stride_functor, CountingIterator> TransformIterator; typedef typename thrust::permutation_iterator<Iterator,TransformIterator> PermutationIterator; // type of the strided_range iterator typedef PermutationIterator iterator; // construct strided_range for the range [first,last) strided_range(Iterator first, Iterator last, difference_type stride) : first(first), last(last), stride(stride) {} iterator begin(void) const { return PermutationIterator(first, TransformIterator(CountingIterator(0), stride_functor(stride))); } iterator end(void) const { return begin() + ((last - first) + (stride - 1)) / stride; } protected: Iterator first; Iterator last; difference_type stride; }; int main(void) { const int Nrows = 5; const int Ncols = 5; const int STRIDE = Nrows + 1; // --- Setting the host, Nrows x Ncols matrix double h_A[Nrows][Ncols] = { { 2., -2., -2., -2., -2.,}, {-2., 4., 0., 0., 0.,}, {-2., 0., 6., 2., 2.,}, {-2., 0., 2., 8., 4.,}, {-2., 0., 2., 4., 10.,} }; // --- Setting the device matrix and moving the host matrix to the device double *d_A; gpuErrchk(cudaMalloc(&d_A, Nrows * Ncols * sizeof(double))); gpuErrchk(cudaMemcpy(d_A, h_A, Nrows * Ncols * sizeof(double), cudaMemcpyHostToDevice)); // --- cuSOLVE input/output parameters/arrays int work_size = 0; int *devInfo; gpuErrchk(cudaMalloc(&devInfo, sizeof(int))); // --- CUDA solver initialization cusolverDnHandle_t solver_handle; cusolverDnCreate(&solver_handle); // --- CUDA CHOLESKY initialization cusolveSafeCall(cusolverDnDpotrf_bufferSize(solver_handle, CUBLAS_FILL_MODE_LOWER, Nrows, d_A, Nrows, &work_size)); // --- CUDA POTRF execution double *work; gpuErrchk(cudaMalloc(&work, work_size * sizeof(double))); cusolveSafeCall(cusolverDnDpotrf(solver_handle, CUBLAS_FILL_MODE_LOWER, Nrows, d_A, Nrows, work, work_size, devInfo)); int devInfo_h = 0; gpuErrchk(cudaMemcpy(&devInfo_h, devInfo, sizeof(int), cudaMemcpyDeviceToHost)); if (devInfo_h != 0) std::cout << "Unsuccessful potrf execution\n\n"; cusolverDnDestroy(solver_handle); // --- Strided reduction of the elements of d_A: calculating the product of the diagonal of the Cholesky factorization thrust::device_ptr<double> dev_ptr = thrust::device_pointer_cast(d_A); typedef thrust::device_vector<double>::iterator Iterator; strided_range<Iterator> pos(dev_ptr, dev_ptr + Nrows * Ncols, STRIDE); double det = thrust::reduce(pos.begin(), pos.end(), 1., thrust::multiplies<double>()); det = det * det; printf("determinant = %f\n", det); return 0; } не позволяет использовать динамический параллелизм.

Вот код:

----+:=8=:+----
ответил JackOLantern 7 AMpTue, 07 Apr 2015 10:15:14 +030015Tuesday 2015, 10:15:14
0

Вы можете использовать OpenCL или CUDA в качестве библиотеки и написать короткую программу (ядро в OpenCL), которая вычисляет определитель на GPU

CUDA http://www.nvidia.de/object/cuda_home_new_de.html

OpenCL http://www.khronos.org/opencl/

http://www.csd.uwo.ca/~ Moreno /Публикации /DetHpcsPaper-proceedings.pdf

Этот документ должен содержать псевдокод для CUDA

ответил glethien 2 PM00000060000002031 2012, 18:00:20

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132