0

我正在尝试反转由复数组成的矩阵，其中我使用矩阵求逆代码来表示“用户” cuda 矩阵逆高斯约旦在以下链接中发布的实数

代码编译，没有错误，但问题是输出错误！我不知道我哪里错了。任何人都可以，请，帮助。先感谢您！

这是完整的代码：

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#pragma comment(lib, "cuda.lib")
#pragma comment(lib, "cudart.lib")
#include <cuda.h>
#include <math.h>
#include <cuda_runtime.h>
#include <cuda_runtime_api.h>
#include "device_launch_parameters.h"
#include <cublas_v2.h>
#include "cuComplex.h"
#include <complex>

__device__ __host__ cuDoubleComplex  operator*(cuDoubleComplex a, cuDoubleComplex b) { return cuCmul(a,b); }
__device__ __host__ cuDoubleComplex  operator+(cuDoubleComplex a, cuDoubleComplex b) { return cuCadd(a,b); }
__device__ __host__ cuDoubleComplex  operator/(cuDoubleComplex a, cuDoubleComplex b) { return cuCdiv(a,b); }
__device__ __host__ cuDoubleComplex  operator-(cuDoubleComplex a, cuDoubleComplex b) { return cuCsub(a,b); }

using namespace std;

 __global__ void gaussjordan(cuDoubleComplex *A,  cuDoubleComplex *I,int n, int i)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    cuDoubleComplex P;

    if(x<n && y<n)
        if(x>i){ 
            P=A[x*n+i]/A[i*n+i];
            I[x*n+y] = I[x*n+y] - I[i*n+y]*P; 
            if(y>=i){ 
                A[x*n+y] = A[x*n+y] - A[i*n+y]*P;  
            }
        }
 }


 __global__ void dev(cuDoubleComplex *d_A,  cuDoubleComplex *dI, int h)
{
    cuDoubleComplex temp = make_cuDoubleComplex(0,0);
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if(x<h && y<h)
        if( cuCimag(d_A[x*h+x]) != cuCimag(temp)){
            if( cuCreal(d_A[x*h+x]) != cuCreal(temp)){

            dI[x*h+y]  = dI[x*h+y]/d_A[x*h+x];
            d_A[x*h+y] = d_A[x*h+y]/d_A[x*h+x];
            }
        }
    __syncthreads();

}

int main()
{
    int const n = 3;
// creating input
    cuDoubleComplex iL[n*n],L[n*n], I[n*n];

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j ) L[i*n+j] =make_cuDoubleComplex(0,1);
            else L[i*n+j] = make_cuDoubleComplex(0,0);

            printf("%.2f  ", cuCimag(L[i*n+j]));
        }
    printf("\n");
    }
printf("\n");

    cuDoubleComplex *d_A, *d_L, *dI;
    float time;
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    int ddsize = n*n*sizeof(cuDoubleComplex);

    dim3 threadsPerBlock(n/16,n/16); //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    dim3 numBlocks(16,16);     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// memory allocation    
    cudaMalloc( (void**)  &d_A, ddsize);   
    cudaMalloc( (void**)   &dI, ddsize); 

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j) I[i*n+i]=make_cuDoubleComplex(1,0);
                else I[i*n+j]=make_cuDoubleComplex(0,0);
        }
    }

 //copy data from GPU to CPU
    cudaMemcpy(  d_A,    L, ddsize, cudaMemcpyHostToDevice); 
    cudaMemcpy(   dI,    I, ddsize, cudaMemcpyHostToDevice); 
//timer start
    cudaEventRecord( start, 0);
// L^(-1)    
    for(int i=0;i<n;i++){
        gaussjordan<<<numBlocks,threadsPerBlock>>>(d_A, dI, n, i);
    }
    dev<<<numBlocks,  threadsPerBlock>>>(d_A, dI, n); 

    cudaMemcpy(iL, dI, ddsize, cudaMemcpyDeviceToHost ); 
    cudaMemcpy(L, d_A, ddsize, cudaMemcpyDeviceToHost ); 

    cudaEventRecord( stop, 0 );
    cudaEventSynchronize( stop );
    cudaEventElapsedTime( &time, start, stop );
    cudaEventDestroy( start );
    cudaEventDestroy( stop );


    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            printf("%.2f  ", cuCimag(iL[i*n+j]));
        }
    printf("\n");
    }
printf("\n");



    std::cout<<"Cuda Time - inverse: "<< time <<"ms\n";

    cudaFree(d_A);
    cudaFree(dI);

    system("Pause");
 return 0;
}

感谢@RobertCrovella 为您提供快速且非常有见地的建议！关于您对我的问题的回答：我更改了我的 threadsPerBlock(4,4) 和 numBlocks(1,1) 所以我将使用 1 个块和 16 个线程作为我的 4x4 矩阵。我的输入矩阵如下

1  0  0  0
0  2  0  0 
0  0  3  0
0  0  0  4

这里所有的数字都是实数，那么预期的倒矩阵应该看起来像

1   0    0   0
0   1/2  0   0 
0   0    1/3 0
0   0    0   1/4

我一点也不明白。我输入了 cuda memcheck 工具来查看我的内核是否没有在午餐，但它没有显示任何错误消息。我最近开始学习 CUDA，没有太多经验。谁能给出更详细的答复？谢谢你！

这是我修改后的代码。

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#pragma comment(lib, "cuda.lib")
#pragma comment(lib, "cudart.lib")
#include <cuda.h>
#include <math.h>
#include <cuda_runtime.h>
#include <cuda_runtime_api.h>
#include "device_launch_parameters.h"
#include <cublas_v2.h>
#include "cuComplex.h"
#include <complex>

__device__ __host__ cuDoubleComplex  operator*(cuDoubleComplex a, cuDoubleComplex b) { return cuCmul(a,b); }
__device__ __host__ cuDoubleComplex  operator+(cuDoubleComplex a, cuDoubleComplex b) { return cuCadd(a,b); }
__device__ __host__ cuDoubleComplex  operator/(cuDoubleComplex a, cuDoubleComplex b) { return cuCdiv(a,b); }
__device__ __host__ cuDoubleComplex  operator-(cuDoubleComplex a, cuDoubleComplex b) { return cuCsub(a,b); }

using namespace std;

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
   if (code != cudaSuccess) 
   {
      fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}




 __global__ void gaussjordan(cuDoubleComplex *A,  cuDoubleComplex *I,int n, int i)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    cuDoubleComplex P;

    if(x<n && y<n)
        if(x>i){ 
            P=A[x*n+i]/A[i*n+i];
            I[x*n+y] = I[x*n+y] - I[i*n+y]*P; 
            if(y>=i){ 
                A[x*n+y] = A[x*n+y] - A[i*n+y]*P;  
            }
        }
 }


 __global__ void dev(cuDoubleComplex *d_A,  cuDoubleComplex *dI, int h)
{
    cuDoubleComplex temp = make_cuDoubleComplex(0,0);
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if(x<h && y<h)
        if( cuCimag(d_A[x*h+x]) != 0 ){
            if( cuCreal(d_A[x*h+x]) != 0 ){

            dI[x*h+y]  = dI[x*h+y]/d_A[x*h+x];
            d_A[x*h+y] = d_A[x*h+y]/d_A[x*h+x];
            }
        }
    __syncthreads();
}

int main()
{
    int const n= 4;
// creating input
    cuDoubleComplex iL[n*n],L[n*n], I[n*n];

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j ) L[i*n+j] =make_cuDoubleComplex(i+1,0);
            else L[i*n+j] = make_cuDoubleComplex(0,0);

            printf("%.2f ", cuCreal(L[i*n+j]));
        }
    printf("\n");
    }
printf("\n");

    cuDoubleComplex *d_A, *dI;
    float time;
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    int ddsize = n*n*sizeof(cuDoubleComplex);

    dim3 threadsPerBlock(n,n); //!!!!!!!!!!!!!!!!!!
    dim3 numBlocks(1,1);       //!!!!!!!!!!!!!!!!!!

// memory allocation    
    cudaMalloc( (void**)  &d_A, ddsize);   
    cudaMalloc( (void**)   &dI, ddsize); 

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j) I[i*n+i]=make_cuDoubleComplex(1,0);
                else I[i*n+j]=make_cuDoubleComplex(0,0);
        }
    }

 //copy data from GPU to CPU
    cudaMemcpy(  d_A,    L, ddsize, cudaMemcpyHostToDevice); 
    cudaMemcpy(   dI,    I, ddsize, cudaMemcpyHostToDevice); 
//timer start
    cudaEventRecord( start, 0);
// L^(-1)    
    for(int i=0;i<n;i++){
        gaussjordan<<<numBlocks,threadsPerBlock>>>(d_A, dI, n, i);
        gpuErrchk( cudaPeekAtLastError() );
    }
    dev<<<numBlocks,  threadsPerBlock>>>(d_A, dI, n); 

    gpuErrchk( cudaPeekAtLastError() );

    gpuErrchk(cudaMemcpy(iL, dI, ddsize, cudaMemcpyDeviceToHost )); 
    gpuErrchk(cudaMemcpy(L, d_A, ddsize, cudaMemcpyDeviceToHost )); 

    cudaEventRecord( stop, 0 );
    cudaEventSynchronize( stop );
    cudaEventElapsedTime( &time, start, stop );
    cudaEventDestroy( start );
    cudaEventDestroy( stop );


    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            printf("%.2f ", cuCreal(iL[i*n+j]));
        }
    printf("\n");
    }
printf("\n");



    std::cout<<"Cuda Time - inverse: "<< time <<"ms\n";

    cudaFree(d_A);
    cudaFree(dI);

    system("Pause");
 return 0;
}
4

1 回答 1

1

免责声明：我不是矩阵求逆的专家。我还没有详细研究实矩阵求逆和复矩阵求逆之间的差异（我不认为应该有很多差异）。正如已经建议的那样，可能有更好/更快的方法来反转矩阵。

直接的问题似乎出在您的dev内核中，尤其是在这里：

    if( cuCimag(d_A[x*h+x]) != cuCimag(temp)){
        if( cuCreal(d_A[x*h+x]) != cuCreal(temp)){

这要求所讨论的矩阵元素的实部和虚部都非零d_A，以便dev内核进行任何工作。但是我不认为这个条件应该是必要的。对于除法，我们可能只要求实部或虚部非零。我认为在复数域中，只有当实部和虚部都为零时，我们实际上才除以零。如果您检查 中cuCdiv提供的功能cuComplex.h，您可以自己确定在什么条件下它会“爆炸”，因此需要测试和避免什么条件。我相信你的测试不正确。

对于您的简单测试用例，以下修改后的代码对我来说可以正常工作：

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include "cuComplex.h"
#include <complex>

__device__ __host__ cuDoubleComplex  operator*(cuDoubleComplex a, cuDoubleComplex b) { return cuCmul(a,b); }
__device__ __host__ cuDoubleComplex  operator+(cuDoubleComplex a, cuDoubleComplex b) { return cuCadd(a,b); }
__device__ __host__ cuDoubleComplex  operator/(cuDoubleComplex a, cuDoubleComplex b) { return cuCdiv(a,b); }
__device__ __host__ cuDoubleComplex  operator-(cuDoubleComplex a, cuDoubleComplex b) { return cuCsub(a,b); }

using namespace std;

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
   if (code != cudaSuccess)
   {
      fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}




 __global__ void gaussjordan(cuDoubleComplex *A,  cuDoubleComplex *I,int n, int i)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    cuDoubleComplex P;

    if(x<n && y<n)
        if(x>i){
            P=A[x*n+i]/A[i*n+i];
            I[x*n+y] = I[x*n+y] - I[i*n+y]*P;
            if(y>=i){
                A[x*n+y] = A[x*n+y] - A[i*n+y]*P;
            }
        }
 }


 __global__ void dev(cuDoubleComplex *d_A,  cuDoubleComplex *dI, int h)
{
    cuDoubleComplex temp = make_cuDoubleComplex(0,0);
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if(x<h && y<h)
        if(( cuCimag(d_A[x*h+x]) != 0 ) || ( cuCreal(d_A[x*h+x]) != 0 )){

            dI[x*h+y]  = dI[x*h+y]/d_A[x*h+x];
            d_A[x*h+y] = d_A[x*h+y]/d_A[x*h+x];

        }
    __syncthreads();
}

int main()
{
    int const n= 4;
// creating input
    cuDoubleComplex iL[n*n],L[n*n], I[n*n];

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j ) L[i*n+j] =make_cuDoubleComplex(i+1,0);
            else L[i*n+j] = make_cuDoubleComplex(0,0);

            printf("%.2f ", cuCreal(L[i*n+j]));
        }
    printf("\n");
    }
printf("\n");

    cuDoubleComplex *d_A, *dI;
    float time;
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    int ddsize = n*n*sizeof(cuDoubleComplex);

    dim3 threadsPerBlock(n,n); //!!!!!!!!!!!!!!!!!!
    dim3 numBlocks(1,1);       //!!!!!!!!!!!!!!!!!!

// memory allocation
    cudaMalloc( (void**)  &d_A, ddsize);
    cudaMalloc( (void**)   &dI, ddsize);

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j) I[i*n+i]=make_cuDoubleComplex(1,0);
                else I[i*n+j]=make_cuDoubleComplex(0,0);
        }
    }

 //copy data from GPU to CPU
    cudaMemcpy(  d_A,    L, ddsize, cudaMemcpyHostToDevice);
    cudaMemcpy(   dI,    I, ddsize, cudaMemcpyHostToDevice);
//timer start
    cudaEventRecord( start, 0);
// L^(-1)
    for(int i=0;i<n;i++){
        gaussjordan<<<numBlocks,threadsPerBlock>>>(d_A, dI, n, i);
        gpuErrchk( cudaPeekAtLastError() );
    }
    dev<<<numBlocks,  threadsPerBlock>>>(d_A, dI, n);

    gpuErrchk( cudaPeekAtLastError() );

    gpuErrchk(cudaMemcpy(iL, dI, ddsize, cudaMemcpyDeviceToHost ));
    gpuErrchk(cudaMemcpy(L, d_A, ddsize, cudaMemcpyDeviceToHost ));

    cudaEventRecord( stop, 0 );
    cudaEventSynchronize( stop );
    cudaEventElapsedTime( &time, start, stop );
    cudaEventDestroy( start );
    cudaEventDestroy( stop );


    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            printf("%.2f ", cuCreal(iL[i*n+j]));
        }
    printf("\n");
    }
printf("\n");



    std::cout<<"Cuda Time - inverse: "<< time <<"ms\n";

    cudaFree(d_A);
    cudaFree(dI);

 return 0;
}

最后免责声明：我并不是说这是一种完全验证的任意维度矩阵求逆方法。我只是指出一个严重的错误，它似乎使它在您的简单测试用例中失败。我对你所链接的上一个问题也表达了一些保留意见。

于 2014-08-29T21:43:22.407 回答