Method MultTrans

ngscuda/cuda_linalg.cpp:745–803

  }

  void DevBlockDiagonalMatrixSoA :: MultTrans (const BaseVector & x, BaseVector & y) const
  {
    static Timer t("DevBlockDiagonalMatrixSoA::MultTrans"); RegionTimer reg(t);
    
    UnifiedVectorWrapper ux(x);
    UnifiedVectorWrapper uy(y);
    ux.UpdateDevice();
    uy.UpdateDevice();
      
    FlatMatrix<Dev<double>> a(dimx*dimy, blocks, (Dev<double>*)dev_data);
    FlatMatrix<Dev<double>> b(dimy, blocks,  (Dev<double>*)ux.DevData());
    FlatMatrix<Dev<double>> res(dimx, blocks,  (Dev<double>*)uy.DevData());


    {
      static Timer t("DevBlockDiagonalMatrixSoA::MultTrans");
      CudaRegionTimer rt(t);

      /*
    DeviceParallelFor
      (res.Width(),
       [a,b,res,inds=FlatArray(indices)] DEVICE_LAMBDA (auto i)
       {
         for (int j = 0; j < res.Height(); j++)
           res(j,i) = 0;
         
         for (int j = 0; j < inds.Size(); j+=3)
           {
             int rowa = inds[j];
             int rowb = inds[j+2];
             int rowres = inds[j+1];
             res(rowres,i) += a(rowa,i) * b(rowb,i);
           }
       });
      */
      
    DeviceParallelFor
      (res.Width(),
       [a,b,res,dimx=this->dimx,sparseT=FlatTable<int>(sparseT)] DEVICE_LAMBDA (auto i)
       {
         for (int j = 0; j < sparseT.Size(); j++)
           {
             double sum = 0;
             for (int k = 0; k < sparseT[j].Size(); k++)
               {
                 int ind = sparseT[j][k];
                 sum += a(ind*dimx+j,i) * b(ind,i);
               }
             res(j,i) = sum;
           }
       });


    
    }
    
    if (synckernels) cudaDeviceSynchronize();
    uy.InvalidateHost();