如何在CUDA/cublas中转置矩阵？

Question

假设我有一个矩阵在GPU上，其维度为A*B，其中B（列数）是以C风格为前导维度。在CUDA（或cublas）中是否有一种方法将此矩阵转置为FORTRAN风格，其中A（行数）成为前导维度？

如果能在主机->设备传输期间进行转置并保持原始数据不变，那就更好了。

- Hailiang Zhang

---

由于CUBLAS可以在转置和正常矩阵上操作，即使在使用行主序的矩阵时，您可能不需要显式计算矩阵转置。 - talonmies

Answer 1

如标题所述，要转置一个按行存储的矩阵 A[m][n]，可以按照以下方式进行：

    float* clone = ...;//copy content of A to clone
    float const alpha(1.0);
    float const beta(0.0);
    cublasHandle_t handle;
    cublasCreate(&handle);
    cublasSgeam( handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, &alpha, clone, n, &beta, clone, m, A, m );
    cublasDestroy(handle);

为了将两个行主元矩阵 A[m][k] 和 B[k][n] 相乘，得到 C=A*B。

    cublasSgemm( handle, CUBLAS_OP_N, CUBLAS_OP_N, n, m, k, &alpha, B, n, A, k, &beta, C, n );

其中C也是行优先矩阵。

- Feng Wang

---

不使用克隆技术，能否完成这个任务？ - bge0

显然，文档中写了一些关于“原地”计算的指令，但我没有找到一种不使用临时矩阵的方法来执行它。希望在整个过程中我的矩阵都比较小。 - Romain Laneuville

Answer 2

CUDA SDK包括一个矩阵转置（matrix transpose），你可以在这里（here）看到实现它的代码示例，从朴素的实现到优化版本。

例如：

朴素的转置

__global__ void transposeNaive(float *odata, float* idata,
int width, int height, int nreps)
{
    int xIndex = blockIdx.x*TILE_DIM + threadIdx.x;
    int yIndex = blockIdx.y*TILE_DIM + threadIdx.y;
    int index_in = xIndex + width * yIndex;
    int index_out = yIndex + height * xIndex;

    for (int r=0; r < nreps; r++)
    {
        for (int i=0; i<TILE_DIM; i+=BLOCK_ROWS)
        {
          odata[index_out+i] = idata[index_in+i*width];
        }
    }
}

正如talonmies所指出的，您可以指定是否要将矩阵作为转置进行操作，在cublas矩阵操作中，例如：对于cublasDgemm()，其中C = a * op(A) * op(B) + b * C，假设您想要将A作为转置（A^T）进行操作，则可以在参数上指定它是（'N'普通或'T'转置）

- dreamcrash

---

谢谢您提供的论文！ - bge0

Answer 3

使用CUDA 5工具包捆绑的CUBLAS版本包含类似BLAS的方法（cublasgeam），可用于转置矩阵。文档在这里查看。

- Talal

Answer 4

这是一个可用的示例：

#include "cublas_v2.h"
#include <vector>
#include <iostream>
using std::cout;

void print_matrix(float *data, int rows, int cols) {
    cout << "[";
    for( int row=0; row < rows; row++) {
        cout << "[";
        for( int col=0; col < cols; col++) {
            cout << data[row*cols+col] << ",";
        }
        cout << "]";
    }
    cout << "]";
}

int main() {
    // allocate host vector
    std::vector<float> h_a = {1,2,3,4,5,6,7,8,9,10};
    int nbytes=h_a.size()*sizeof(*h_a.data());
    std::vector<float> h_b(h_a.size());

    // define the number or rows and the number of columns
    int m=2,n=5;

    // allocate device vectors
    float *d_a, *d_b;
    cudaMalloc(&d_a, nbytes);
    cudaMalloc(&d_b, nbytes);

    // copy host vector to device
    cudaMemcpy(d_a,h_a.data(), nbytes, cudaMemcpyHostToDevice);

    // perform a transpose
    {

        float alpha=1;
        float *A=d_a;
        int lda=n;

        float beta=0;
        float *B=NULL;
        int ldb=n;

        float *C=d_b;
        int ldc=m;

        cublasHandle_t handle;
        cublasCreate(&handle);
        cublasStatus_t success=cublasSgeam( handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, &alpha, A, lda, &beta, B, ldb, C, ldc);
        if ( success != CUBLAS_STATUS_SUCCESS)
            cout << "\33[31mError: " << success << "\33[0m\n";
        cublasDestroy(handle);
    }

    // copy back to host
    cudaMemcpy(h_b.data(),d_b,nbytes,cudaMemcpyDeviceToHost);

    cout << "origional:  ";
    print_matrix(h_a.data(),m,n);
    cout << "\n";

    cout << "transposed: ";
    print_matrix(h_b.data(),n,m);
    cout << "\n";

    cudaFree(d_a);
    cudaFree(d_b);
    return 0;
}

- philn