Лучший способ добиться диагональности вектора CUDA

Question

Что я хочу сделать, это подать в мою матрицу mxn и параллельно построить n квадратных диагональных матриц для каждого столбца матрицы, выполнить операцию на каждой квадратной диагональной матрице и затем рекомбинировать результат. Как мне это сделать?

До сих пор я начинал с матрицы m x n; результат предыдущего вычисления матрицы, где каждый элемент вычисляется с использованием функции y = f (g (x)).

Это дает мне матрицу с n элементами столбца [f1, f2 ... fn], где каждый fn представляет собой вектор-столбец высоты m.

Отсюда я хочу различать каждый столбец матрицы относительно g (x). Дифференцируя fn (x) w.r.t. g (x) приводит к квадратной матрице с элементами f '(x). При ограничении эта квадратная матрица сводится к якобиану с элементами каждой строки по диагонали квадратной матрицы и равна fn ', а все остальные элементы равны нулю.

Отсюда причина, по которой необходимо построить диагональ для каждой из векторных строк fn.

Для этого я беру целевой вектор, определенный как A (hA x 1), который был извлечен из большей матрицы A (m x n). Затем я подготовил нулевую матрицу, определенную как C (hA x hA), которая будет использоваться для хранения диагоналей.

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

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

Код ядра (который работает), чтобы выполнить это показано здесь:

_cudaDiagonalizeTest << <5, 1 >> >(d_A, matrix_size.uiWA, matrix_size.uiHA, d_C, matrix_size.uiWC, matrix_size.uiHC); 

__global__ void _cudaDiagonalizeTest(float *A, int wA, int hA, float *C, int wC, int hC) 
{ 
    int ix, iy, idx; 

    ix = blockIdx.x * blockDim.x + threadIdx.x; 
    iy = blockIdx.y * blockDim.y + threadIdx.y; 

    idx = iy * wA + ix; 

    C[idx * (wC + 1)] = A[idx]; 

} 

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

а) снижение

б) тяги

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

На рисунке ниже показан желаемый результат.

Я прочитал NVIDIA's article on reduction, но не смог достичь желаемых результатов.

Любая помощь или объяснение будет очень приветствоваться.

Спасибо.

Матрица А мишень с 4-мя колоннами. Я хочу взять каждый столбец и скопировать его элементы в матрицу B как диагональ, итерацию через каждый столбец.

Answer 1

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

Другой подход может быть основан на thrust strided_range example.

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

#include <thrust/device_vector.h> 
#include <thrust/scatter.h> 
#include <thrust/sequence.h> 
#include <thrust/iterator/transform_iterator.h> 
#include <thrust/iterator/counting_iterator.h> 
#include <thrust/functional.h> 
#include <iostream> 


template<typename V> 
void print_matrix(const V& mat, int rows, int cols) 
{ 
    for(int i = 0; i < rows; ++i) 
    { 
    for(int j = 0; j < cols; ++j) 
    { 
     std::cout << mat[i + j*rows] << "\t"; 
    } 
    std::cout << std::endl; 
    } 
} 

struct diag_index : public thrust::unary_function<int,int> 
{ 
    diag_index(int rows) : rows(rows){} 

    __host__ __device__ 
    int operator()(const int index) const 
    { 
     return (index*rows + (index%rows)); 
    } 

    const int rows; 
}; 

int main() 
{ 
    const int rows = 5; 
    const int cols = 4; 

    // allocate memory and fill with demo data 
    // we use column-major order 
    thrust::device_vector<int> A(rows*cols); 
    thrust::sequence(A.begin(), A.end()); 

    thrust::device_vector<int> B(rows*rows*cols, 0); 

    // fill diagonal matrix 
    thrust::scatter(A.begin(), A.end(), thrust::make_transform_iterator(thrust::make_counting_iterator(0),diag_index(rows)), B.begin()); 

    print_matrix(A, rows, cols); 
    std::cout << std::endl; 
    print_matrix(B, rows, rows*cols); 
    return 0; 
} 

Этот пример выведет:

0 5 10 15  
1 6 11 16  
2 7 12 17  
3 8 13 18  
4 9 14 19  

0 0 0 0 0 5 0 0 0 0 10 0 0 0 0 15 0 0 0 0  
0 1 0 0 0 0 6 0 0 0 0 11 0 0 0 0 16 0 0 0  
0 0 2 0 0 0 0 7 0 0 0 0 12 0 0 0 0 17 0 0  
0 0 0 3 0 0 0 0 8 0 0 0 0 13 0 0 0 0 18 0  
0 0 0 0 4 0 0 0 0 9 0 0 0 0 14 0 0 0 0 19  

Answer 2

Альтернативный ответ, который не использует тягу выглядит следующим образом:

_cudaMatrixTest << <5, 5 >> >(d_A, matrix_size.uiWA, matrix_size.uiHA, d_C, matrix_size.uiWC, matrix_size.uiHC); 

__global__ void _cudaMatrixTest(float *A, int wA, int hA, float *C, int wC, int hC) 
{ 
    int ix, iy, idx; 

    ix = blockIdx.x * blockDim.x + threadIdx.x; 
    iy = blockIdx.y * blockDim.y + threadIdx.y; 

    idx = iy * wA + ix; 

    C[idx * wC + (idx % wC)] = A[threadIdx.x * wA + (ix/wC)]; 
} 

где D_A является

0 5 10 15  
1 6 11 16  
2 7 12 17  
3 8 13 18  
4 9 14 19  

Оба ответа являются жизнеспособными решениями. Вопрос в том, что лучше/быстрее?