Противоречие cublasDgetrfBatched и cublasDtrsmBatched, когда нужно решать массив линейных систем с использованием cuBLAS

Question

У меня много плотных линейных систем, которые я хочу решить в пакетном формате cuBLAS. Так что мой план

1. использовать cublasDgetrfBatched для пакетной декомпозиции LU

2. Затем используйте cublasDtrsmBatched для пакетной нижней треугольной и пакетной верхней треугольной части один за другим.

Код дается как

 #include<stdio.h>
 #include<stdlib.h>
 #include<cuda_runtime.h> 
 #include<device_launch_parameters.h>
 #include<cublas_v2.h>

 const int N = 32;
 const int Nmatrices = N;

 __global__ void initiate(double *d_A, double *d_B)
    {
     int i = threadIdx.x;       int j = blockIdx.x;

     int id = j*N*N + i*N;      int idb = j*N + i;  

     for(int k = 0; k< N ; k++)
         {
        d_A[id + k] = 0.0;

        if(k == i-2)    d_A[id + k] = 1.0; 
        if(k == i-1)    d_A[id + k] = 2.0; 
        if(k == i)      d_A[id + k] = 8.0;
        if(k == i+1)    d_A[id + k] = 2.0;
        if(k == i+2)    d_A[id + k] = 1.0;  
        } 
     d_B[idb] = 8.0;   
    }


int main()
 {
    cublasHandle_t handle;      cublasSafeCall(cublasCreate(&handle));

// Allocate device space for the input matrices
  double *d_A_sys; cudaMalloc((void**)&d_A_sys, N*N*Nmatrices*sizeof(double));
  double *d_B_sys; cudaMalloc((void**)&d_B_sys, N*Nmatrices  *sizeof(double));

// Allocate host space for the solution
  double *h_B_sys = (double *)malloc(N*Nmatrices*sizeof(double));

// kernel for initiat d_A_sys and d_B_sys
  initiate<<<Nmatrices, N>>>(d_A_sys, d_B_sys);

//Creating the array of pointers needed as input/output to the batched getrf
  double **h_A_pointers = (double **)malloc(Nmatrices*sizeof(double *));
  for (int i = 0; i < Nmatrices; i++) h_A_pointers[i] = d_A_sys + i*N*N;

  double **h_b_pointers = (double **)malloc(Nmatrices*sizeof(double *));
  for (int i = 0; i < Nmatrices; i++) h_B_pointers[i] = d_B_sys + i*N;

  double **d_A_pointers;
  cudaMalloc((void**)&d_A_pointers, Nmatrices*sizeof(double *));
  cudaMemcpy(d_A_pointers, h_A_pointers, Nmatrices*sizeof(double *), cudaMemcpyHostToDevice);


  double **d_b_pointers;
  cudaMalloc((void**)&d_b_pointers, Nmatrices*sizeof(double *));
  cudaMemcpy(d_b_pointers, h_b_pointers, Nmatrices*sizeof(double *), cudaMemcpyHostToDevice);

  int *d_InfoArrays; cudaMalloc((void**)&d_InfoArrays,  Nmatrices*sizeof(int));
  int *h_InfoArrays = (int *)malloc(Nmatrices*sizeof(int));

//Batched LU decomposition
   cublasDgetrfBatched(handle, N, d_A_pointers, N, NULL, d_InfoArrays, Nmatrices));

  //Batched Lower triangular part
  cublasDtrsmBatched(handle, 
                   CUBLAS_SIDE_LEFT, 
                   CUBLAS_FILL_MODE_LOWER,
                   CUBLAS_OP_N,
                   CUBLAS_DIAG_UNIT,
                   N,
                   N,
                   &alpha,
                   d_A_pointers,
                   N,
                   d_b_pointers,
                   N,
                   Nmatrices);

  //Batched Upper triangular part
  cublasDtrsmBatched(handle,
                   CUBLAS_SIDE_LEFT,
                   CUBLAS_FILL_MODE_UPPER,
                   CUBLAS_OP_N,
                   CUBLAS_DIAG_NON_UNIT,
                   N,
                   N,
                   &alpha,
                   d_A_pointers,
                   N,
                   d_b_pointers,
                   N,
                   Nmatrices);

  cudaMemcpy(h_B_sys, d_B_sys, N*Nmatrices*sizeof(double), cudaMemcpyDeviceToHost);
  printf("Print out the solutions \n");

  cublasDestroy(handle);
  gpuErrchk(cudaDeviceReset());
  return 0;
 }

Требование cublasDgetrfBatched и cublasDtrsmBatched d_A_pointers должно быть в двойном типе, но когда я выполняю, более поздний, который дает мне ошибку компиляции как это, видит рис:

Как преодолеть проблему, любая помощь?

Answer 1

Вы можете решить проблему правильности const, выполнив что-то вроде этого:

const double **d_A_pointers;
cudaMalloc((void**)&d_A_pointers, Nmatrices*sizeof(double *));

....

//Batched LU decomposition
cublasDgetrfBatched(handle, N, const_cast<double**>(d_A_pointers), N, NULL, d_InfoArrays, Nmatrices);

то есть отбросить постоянство в cublasDgetrfBatched

Другая вещь, которую вы явно ошибаетесь в размещенном коде - это размеры входов в cublasDtrsmBatched звонки. Я считаю, что вы хотите что-то вроде:

//Batched Lower triangular part
cublasDtrsmBatched(handle, 
        CUBLAS_SIDE_LEFT, 
        CUBLAS_FILL_MODE_LOWER,
        CUBLAS_OP_N,
        CUBLAS_DIAG_UNIT,
        N,
        1,
        &alpha,
        d_A_pointers,
        N,
        d_b_pointers,
        N,
        Nmatrices);

т. е. входной размер матрицы RHS для вашего примера - Nx1, а не NxN (вы не решаете N проблем RHS на разложенную матрицу, только одну).

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