Skip to content

Instantly share code, notes, and snippets.

@vertexoperator vertexoperator/cgcuda.md
Last active Jul 24, 2018

Embed
What would you like to do?

CUDAの共役勾配法サンプルが遅かった話

Conjugate Gradient Method

  • 正定値対称行列$A$とベクトル$b$に対して$A x = b$を解くアルゴリズム

  • $A$が疎行列の時に使うことが多い。前処理付きでやるのが普通

  • 正定値対称疎行列に対する連立一次方程式ソルバーとしては、今でも(多分)最強の方法

  • 実は、正定値じゃなくても、大抵答えが出るらしい

cf)不定値対称行列に対する共役勾配法の収束について

https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/42071

Eigen3.3.4の実装

前処理付きConjugate Gradientメインループ

  while(i < maxIters)
  {
    tmp.noalias() = mat * p;                    // SpMV(sparse matrix-vector multiplication)

    Scalar alpha = absNew / p.dot(tmp);         // the amount we travel on dir
    x += alpha * p;                             // update solution
    residual -= alpha * tmp;                    // update residual
    
    residualNorm2 = residual.squaredNorm();
    if(residualNorm2 < threshold)
      break;
    
    z = precond.solve(residual);                // approximately solve for "A z = residual"

    RealScalar absOld = absNew;
    absNew = numext::real(residual.dot(z));     // update the absolute value of r
    RealScalar beta = absNew / absOld;          // calculate the Gram-Schmidt value used to create the new search direction
    p = z + beta * p;                           // update search direction
    i++;
  }

precondがEigen::DiagonalPreconditionerの場合

1ループ当たり浮動小数点演算数 = $2 N_{nz} + 14 N_{dim} + 2$

$N_{nz}$:疎行列の非零要素数

$N_{dim}$:行列の次元

普通、$N_{nz} >> N_{dim}$なので、演算数的にSpMVがボトルネック

SpMVは、iterationごとに間接メモリ参照が入るので、メモリ律速になる

CUDA7.5のサンプル実装

メインループ(前処理なし。Eigenとの対応が見やすいように変数名など微修正)

    while (absNew > threshold && k <= max_iter)
    {
        if (k > 1)
        {
            beta = absNew / absOld;
            //p = r + beta*p
            cublasStatus = cublasSscal(cublasHandle, N, &beta, d_p, 1);
            cublasStatus = cublasSaxpy(cublasHandle, N, &one , d_r, 1, d_p, 1);
        }
        else
        {
            cublasStatus = cublasScopy(cublasHandle, N, d_r, 1, d_p, 1);  //p=r
        }
        //Ax = A*p
        cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, Nnz, &one, descr, d_val, d_row, d_col, d_p, &zero, d_Ax);
        cublasStatus = cublasSdot(cublasHandle, N, d_p, 1, d_Ax, 1, &dot);        //dot = sdot(p,Ax);
        float alpha = absNew / dot;

        cublasStatus = cublasSaxpy(cublasHandle, N, &alpha, d_p, 1, d_x, 1);      //x += alpha*p;
        float nalpha = -alpha;
        cublasStatus = cublasSaxpy(cublasHandle, N, &nalpha, d_Ax, 1, d_r, 1);    //r -= alpha*Ax;

        absOld = absNew;
        cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &absNew);      //absNew = sdot(r,r);
        cudaDeviceSynchronize();
        k++;
    }

CUDA実装1

Eigenと処理の順序を合わせて、前処理を追加

    int B=256 , P=(N+B-1)/B;  //thread per block, number of blocks
    for(k=1;k<=max_iter;k++){
        cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, Nnz, &one, descr, d_val, d_row, d_col, d_p, &zero, d_Ax);
        cublasStatus = cublasSdot(cublasHandle, N, d_p, 1, d_Ax, 1, &dot);
        float alpha = absNew/dot;
        cublasStatus = cublasSaxpy(cublasHandle, N, &alpha, d_p, 1, d_x, 1);
        float nalpha = -alpha;
        cublasStatus = cublasSaxpy(cublasHandle, N, &nalpha, d_Ax, 1, d_r, 1);
        cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &residualNorm2);
        cudaDeviceSynchronize();
        if(residualNorm2 < threshold)
            break;

        precond<<<P,B>>>(d_Mdiag , d_r , d_z , N);             //DiagonalPreconditioner
        float absOld  = absNew;
        cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, d_z, 1, &absNew);
        float beta = absNew/absOld;
        cublasStatus = cublasSscal(cublasHandle, N, &beta, d_p, 1);
        cublasStatus = cublasSaxpy(cublasHandle, N, &one, d_z, 1, d_p, 1);
    }
/*

__global__ void precond(float *A, float *B, float *C, int maxN){
   int i=blockIdx.x*blockDim.x + threadIdx.x;

  if(i<maxN){
      C[i]=A[i]*B[i];
  }

}
*/

$N_{dim}$が3万程度。$N_{nz}$は100万程度のデータ。EigenもCUDAも全部floatを使用

実装 #iter time(ms)
Eigen(12スレッド) 366 122
CUDAサンプル 677 440
CUDA実装1 378 304 前処理追加

profilerによると、cusparseScsrmvはGTX1080上で一回あたり68μsかかっていた(〜30 GFlops)。677回呼ぶと46ms

EigenのSpMVは測定すると一回あたり200μs前後かかっているっぽい。366回呼ぶと73msec

CUDAのサンプルが遅すぎる

CUDA実装2

residualNorm2以外、内積計算の結果をCPUへ戻す必要はない

residualNorm2も、毎回チェックしなくてもよさそう。試しに20回に1回だけチェックするようにしてみる

    cusparseSetPointerMode(cusparseHandle, CUSPARSE_POINTER_MODE_DEVICE);
    cublasSetPointerMode(cublasHandle, CUBLAS_POINTER_MODE_DEVICE);
    for(k=1;k<=max_iter;k++){
        cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, Nnz, addr_d_one, descr, d_val, d_row, d_col, d_p, addr_d_zero, d_Ax);
        cublasStatus = cublasSdot(cublasHandle, N, d_p, 1, d_Ax, 1, addr_d_dot);
        compute_alpha_kernel<<<1,1>>>();
        cublasStatus = cublasSaxpy(cublasHandle, N, addr_d_alpha, d_p, 1, d_x, 1);
        cublasStatus = cublasSaxpy(cublasHandle, N, addr_d_nalpha, d_Ax, 1, d_r, 1);
        cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, addr_d_residualNorm2);
        if(k%20==0){
            cudaDeviceSynchronize();
            cudaMemcpyFromSymbol(&residualNorm2 , d_residualNorm2 , sizeof(float) , 0 ,  cudaMemcpyDeviceToHost);
            if(residualNorm2 < threshold)
                break;
        }
        precond<<<P,B>>>(d_Mdiag , d_r , d_z , N);
        cudaMemcpyToSymbol(d_absOld , addr_d_absNew , sizeof(float) , 0 , cudaMemcpyDeviceToDevice); 
        cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, d_z, 1, addr_d_absNew);
        compute_beta_kernel<<<1,1>>>();
        cublasStatus = cublasSscal(cublasHandle, N, addr_d_beta, d_p, 1);
        cublasStatus = cublasSaxpy(cublasHandle, N, addr_d_one, d_z, 1, d_p, 1);
    }
実装 #iter time(ms)
Eigen(12スレッド) 366 122
CUDAサンプル 677 440
CUDA実装1 378 304 前処理追加
CUDA実装2 380 191 データ移動削減

まだCPUより遅い

CUDA実装3

カーネル呼び出しを減らそうと思って、cublasのカーネルをfusionしようとした

とりあえず、cublasSdotを自前のsdotカーネルに変更

    cusparseSetPointerMode(cusparseHandle, CUSPARSE_POINTER_MODE_DEVICE);
    cublasSetPointerMode(cublasHandle, CUBLAS_POINTER_MODE_DEVICE);
    for(k=1;k<=max_iter;k++){
        cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, Nnz, addr_d_one, descr, d_val, d_row, d_col, d_p, addr_d_zero, d_Ax);
        cudaMemcpyToSymbol(d_dot , addr_d_zero , sizeof(float) , 0 , cudaMemcpyDeviceToDevice); 
        sdot<<<P,B>>>(d_p , d_Ax , addr_d_dot , N);
        compute_alpha_kernel<<<1,1>>>();
        cublasStatus = cublasSaxpy(cublasHandle, N, addr_d_alpha, d_p, 1, d_x, 1);
        cublasStatus = cublasSaxpy(cublasHandle, N, addr_d_nalpha, d_Ax, 1, d_r, 1);
        cudaMemcpyToSymbol(d_residualNorm2 , addr_d_zero , sizeof(float) , 0 , cudaMemcpyDeviceToDevice); 
        sdot<<<P,B>>>(d_r , d_r , addr_d_residualNorm2 , N);
        if(k%20==0){
           cudaDeviceSynchronize();
           cudaMemcpyFromSymbol(&residualNorm2 , d_residualNorm2 , sizeof(float) , 0 ,  cudaMemcpyDeviceToHost);
           if(residualNorm2 < threshold)
              break;
        }
        precond<<<P,B>>>(d_Mdiag , d_r , d_z , N);
        cudaMemcpyToSymbol(d_absOld , addr_d_absNew , sizeof(float) , 0 , cudaMemcpyDeviceToDevice); 
        cudaMemcpyToSymbol(d_absNew , addr_d_zero , sizeof(float) , 0 , cudaMemcpyDeviceToDevice); 
        sdot<<<P,B>>>(d_r , d_z , addr_d_absNew , N);
        compute_beta_kernel<<<1,1>>>();
        cublasStatus = cublasSscal(cublasHandle, N, addr_d_beta, d_p, 1);
        cublasStatus = cublasSaxpy(cublasHandle, N, addr_d_one, d_z, 1, d_p, 1);
    }


/*
__global__ void sdot(float *A , float *B , float *res , int maxN ){
   int i = blockIdx.x*blockDim.x + threadIdx.x;
   __shared__ float cache[256];
   //assert(blockDim.x==256);
   cache[threadIdx.x] = 0.0f;

   while(i < maxN) {
       cache[threadIdx.x] += A[i]*B[i];
       i += gridDim.x * blockDim.x;
   }
   __syncthreads();
   
   for(unsigned int s=1;s<blockDim.x;s*=2) {
       if(threadIdx.x%(2*s)==0) {
           cache[threadIdx.x] += cache[threadIdx.x + s];
       }
       __syncthreads();
   }
   if(threadIdx.x==0)atomicAdd(res , cache[0]);
}
*/
実装 #iter time(ms)
Eigen(12スレッド) 366 122
CUDAサンプル 677 440
CUDA実装1 378 304 前処理追加
CUDA実装2 380 191 データ移動削減
CUDA実装3 380 86 cublasSdot置き換え

cublasSdotが遅かったことが分かった。理由は不明

CUDA実装4

一回のカーネル呼び出しにまとめられるものを、全部まとめた。cublasはいなくなった

    cusparseSetPointerMode(cusparseHandle, CUSPARSE_POINTER_MODE_DEVICE);
    for(k=1;k<=max_iter;k++){
        cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, Nnz, addr_d_one, descr, d_val, d_row, d_col, d_p, addr_d_zero, d_Ax);
        sdot<<<P,B>>>(d_p , d_Ax , addr_d_dot , N);

        //alpha = absNew/dot , absOld = absNew , dot = absNew = residualNorm2 = 0
        cg_update_params<<<1,1>>>();

        //x += alpha*p , r -= alpha*Ax , z = M*r
        cgstep<<<P,B>>>(d_p , d_x, d_r , d_Ax, d_Mdiag , d_z , N);

        //residualNorm2 , absNew = sdot(r,r) , sdot(r,z)
        cg_sdot_fusion<<<P,B>>>(d_r , d_z , N);

        if(k%20==0){
           cudaDeviceSynchronize();
           cudaMemcpyFromSymbol(&residualNorm2 , d_residualNorm2 , sizeof(float) , 0 , cudaMemcpyDeviceToHost);
           if(residualNorm2 < threshold)break;
        }

        //p = z + (absNew/absOld)*p
        sxpby<<<P,B>>>(d_z , d_p , N);
    }
実装 #iter time(ms)
Eigen(12スレッド) 366 122
CUDAサンプル 677 440
CUDA実装1 378 304 前処理追加
CUDA実装2 380 191 データ移動削減
CUDA実装3 380 86 cublasSdot置き換え
CUDA実装4 380 46 カーネルfusion

cusparseScsrmvはGTX1080上で一回あたり68μsかかっていた(〜30 GFlops)。380回呼ぶと25.8ms。

CUDA kernelのoverheadが5〜10μsといわれる。実装4では、SpMV以外にループあたり5個のカーネルがあるので、380回呼ぶと、9.5~19ms

saxpyやsdotも、演算が少なすぎてメモリ律速になる。例えばsaxpyは、2演算あたりfloatのロード2回ストア1回で、演算強度1/6[Flops/Byte]だから、GTX1080の場合、バンド幅320GB/sから、概算で50GFlopsオーダーのはず。

sdotとかも、ざっくり50GFlopsと仮定すると、$N_{dim}$が3万程度だから、ループ当たり8.4μsかかる計算。380回呼ぶと、3.2msとなって、全体への影響は小さそう

全部足すと、38.5~48msくらいで、他にresidualNorm2コピーの時間があるので、実装4は妥当な速度といえる

GTX-1080一枚で、メモリ的に$N_{nz}$は一億以上いけるので、それくらいの規模になったら、サンプル実装でも問題ないと思われる

CUDA実装5

20回に1回終了条件をチェックするのは、計算量がもっと大きい場合にはよくない

といって、ループごとにcudaMemcpyFromSymbolすると、CPUより遅くなってしまう(大体150msくらい)。

pinned memoryを使って折衷案を試みる

    //cudaHostAlloc( (void**) &mapped_mem , sizeof(float) , cudaHostAllocMapped );
    mapped_mem[0] = FLT_MAX;
    cusparseSetPointerMode(cusparseHandle, CUSPARSE_POINTER_MODE_DEVICE);
    for(k=1;k<=max_iter;k++){
        //Ax = A*p;
        cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, Nnz, addr_d_one, descr, d_val, d_row, d_col, d_p, addr_d_zero, d_Ax);
        sdot<<<P,B>>>(d_p , d_Ax , addr_d_dot , N);
        cg_update_params<<<1,1>>>( mapped_mem );
        cudaDeviceSynchronize();
        if(mapped_mem[0] < threshold)break;
        cgstep<<<P,B>>>(d_p , d_x , d_r , d_Ax , d_Mdiag , d_z , N );
        cg_sdot_fusion<<<P,B>>>(d_r , d_z , N );
        sxpby<<<P,B>>>(d_z , d_p , N);
    }

/*
__global__ void cg_update_params(float *buf){
  d_alpha = d_absNew/d_dot;
  d_absOld = d_absNew;
  buf[0] = d_residualNorm2;
  d_dot = d_absNew = d_residualNorm2 = 0.0f;
}

*/
実装 #iter time(ms)
Eigen(12スレッド) 366 122
CUDAサンプル 677 440
CUDA実装1 378 304 前処理追加
CUDA実装2 380 191 データ移動削減
CUDA実装3 380 86 cublasSdot置き換え
CUDA実装4 380 46 カーネルfusion
CUDA実装5 378 98 pinned memory使用

一応、CPUには勝ってるけど、大分遅くなった

補足

CUDA実装4では、最大120個のCUDAカーネルが、GPUのキューにたまる。最大何個まで大丈夫か

#include <stdio.h>

__global__ void forever(){
    for(;;){}
}

int main(){

   for(int k=0;k<100000;++k){
      forever<<<1,1>>>();
      auto err = cudaGetLastError();
      if(err!=cudaSuccess){
        fprintf(stderr, "Failed to launch kernel@k=%d (error code '%s')! \n", k,cudaGetErrorString(err));
        cudaDeviceReset();
        break;
      }
   }
   cudaDeviceSynchronize();
   return 0;
}

安定して、k=1024で失敗する。環境依存かも

Kernel launches blocking when 1024 kernels in a queueでは、かつて、32か64が限界だったこともあると書かれている。

まとめ:計算量が少ない時は収束チェックを減らすのが最適っぽい。但し、カーネル起動しすぎないように注意する必要がある

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.