Function DeviceSparseCholeskySolveLTransKernel

ngscuda/dev_sparsecholesky.cpp:190–331

  __global__ void DeviceSparseCholeskySolveLTransKernel (
        FlatTable<int> dependency, 
        FlatArray<Dev<int>> incomingdep, 
        FlatVector<Dev<double>> hy,
        int & cnt,
        FlatArray<Dev<MicroTask>> microtasks,
        FlatArray<Dev<int>> blocks,
        FlatArray<Dev<int>> rowindex2,
        FlatArray<Dev<size_t>> firstinrow_ri,
        FlatArray<Dev<size_t>> firstinrow,
        FlatArray<Dev<double>> lfact
      )
  {
    __shared__ int myjobs[16]; // max blockDim.y;   
    
    while (true)
       {
          if (threadIdx.x == 0)
             myjobs[threadIdx.y] = atomicAdd(&cnt, 1);
          __syncwarp();

          int myjob = dependency.Size()-1 - myjobs[threadIdx.y];
          if (myjob < 0)
            break;

          volatile int * n_deps = (int*)&incomingdep[myjob];
          if (threadIdx.x == 0)
            while(*n_deps);
          __syncwarp();

          // do the work for myjob....
          
          MicroTask task = microtasks[myjob];
          size_t blocknr = task.blocknr;
          auto range = T_Range<int>(blocks[blocknr], blocks[blocknr+1]);
          auto base = firstinrow_ri[range.First()] + range.Size()-1;
          auto ext_size =  firstinrow[range.First()+1]-firstinrow[range.First()] - range.Size()+1;
          auto extdofs = rowindex2.Range(base, base+ext_size);
    __threadfence();
          // TODO: needs transpose:
          if ((task.type == MicroTask::B_BLOCK) || (task.type == MicroTask::LB_BLOCK))
            {
              if (extdofs.Size() != 0)
                  {
                    auto myr = Range(extdofs).Split (task.bblock, task.nbblocks);
                    auto my_extdofs = extdofs.Range(myr);
                                          
                    // for (auto i : range)
                    for (int i = range.First()+threadIdx.x; i < range.Next(); i += blockDim.x)
                        {
                            size_t first = firstinrow[i] + range.end()-i-1;
                            auto ext_lfact = lfact.Range(first, first+extdofs.Size());
                            
                            double val = 0.0;
                            for (auto j : Range(my_extdofs))
                                val += ext_lfact[myr.begin()+j] * hy(my_extdofs[j]);
                                
                            atomicAdd ((double*)(&hy(i)), -val);