матричная инверсия для комплексных чисел по методу Гаусса-Йорда в cuda

Question

Я пытаюсь инвертировать матрицу, состоящую из комплексных чисел, где я использую код инверсии матрицы для действительных чисел, размещенный по следующей ссылке 'user' cuda matrix inverse gaussian jordan

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

здесь полный код:

#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 

и я не получаю это вообще. Я ввел инструмент mudcheck cuda, чтобы узнать, не работает ли мое ядро ​​ , но он не показывал никаких массовых сообщений об ошибках. Я начал изучать 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; 
} 

Answer 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; 
} 

ОКОНЧАТЕЛЬНЫЙ ОТКАЗ: Я не говорю, что это полностью подтверждено подход к инверсии матриц произвольной размерности. Я просто указываю на критическую ошибку, которая, кажется, делает ее неудачной для вашего простого теста. Я также высказал некоторые оговорки в предыдущем вопросе, который вы связали.