6

小さな（6x6）倍精度行列の行列式を完全にGPUで計算するライブラリまたは無料で利用できるコードはありますか？

4

3 に答える 3

4

これが計画です。これらの小さな行列を数百個バッファリングし、カーネルを1回起動して、それらすべての行列式を一度に計算する必要があります。

実際のコードを書くつもりはありませんが、これは役立つはずです。

1）＃ブロック=＃行列を起動します。各ブロックは、各行列の行列式を計算します。

2）det（A）= det（A11 * A22-A21 * A12）; ここで、Aは6x6、A11、A12、A21、A22はAの3x3サブ行列です。

3）3x3行列の行列乗算を行うデバイス関数を記述します

4）3x3行列の行列式は簡単に計算できます。ここからの式を使用してください。

編集：どうやら（2）はA21 * A12 == A12*A21の場合にのみ機能します

代替案は次のようになります

1）各6x6行列のガウス消去法によるLU分解

2）Uの対角要素を乗算して、行列式を取得します。

于 2012-08-02T18:18:46.593 に答える
4

上記のコメントで特にBartによってすでに指摘されているように、GPUを使用して小さな行列の行列式の計算を計算しても（それらの多くの場合でも）、他のコンピューティングプラットフォームよりもゲインが保証されません。

行列式の計算の問題は、アプリケーションで数回発生する可能性のある興味深い問題だと思います。現在、CUDAを使用した行列式計算用のルーチンを提供しているライブラリを認識していません（そのようなルーチンcuBLASもcuSOLVERありません）。そのため、2つの可能性があります。

  1. Pavanが指摘したように、独自のアプローチを実装する。
  2. 他の利用可能なルーチンを使用するために管理します。

最後の点に関して、1つの可能性は、コレスキー分解を使用してから、コレスキー行列の対角要素の積の2乗として行列式を計算することです。Matlabでは、次のようになります。

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

以下に、このアイデアを利用したコードを提供します。特に、コレスキー分解はcuSOLVERのpotrf関数を使用して実行されますが、コレスキー行列の対角線上の要素の積は、推力によるストライドリダクションのアプリケーションです。

以下のコードは大きな行列用ですので、大きな行列の行列式を計算する必要がある人にとってはすぐに役立ちます。しかし、それをいくつかの小さな行列に適応させる方法は？1つの可能性はcuSOLVER、コレスキー分解にのストリームを使用してから、Thurst1.8の新しい動的並列処理機能を使用することです。CUDA 7.0の時点でcuSOLVERは、動的並列処理の使用が許可されていないことに注意してください。

コードは次のとおりです。

#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;
}
于 2015-04-07T07:15:14.800 に答える
-2

OpenCLまたはCUDAをライブラリとして使用し、GPUで行列式を計算する短いプログラム（OpenCLのカーネル）を作成できます。

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

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

http://www.csd.uwo.ca/~moreno/Publications/DetHpcsPaper-proceedings.pdf

このペーパーには、CUDAの擬似コードが含まれている必要があります

于 2012-08-02T14:00:20.737 に答える