vertexoperator/cgcuda.md

## cgcuda.md

      
    Raw
  

              cgcuda.md
            
          
    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が限界だったこともあると書かれている。
まとめ:計算量が少ない時は収束チェックを減らすのが最適っぽい。但し、カーネル起動しすぎないように注意する必要がある
実装	#iter	time(ms)
Eigen(12スレッド)	366	122
CUDAサンプル	677	440
CUDA実装1	378	304	前処理追加