Covariance of 3 dimensional sparse boolean array in Cublas/cuSPARSE

Hello I have sparse 3 dimensional boolean arrays (data - ones - are concentrated in relatively small area somwhere in the array - and array has traditional dense array format) and I need to calculate its covariance.
How to achieve it in this sparse, 3 dimensional setting, using for example CuSparse and CuBLAS?

For example How It is calculated in python

       gt_n = np.count_nonzero(self.reference)
        gt_indices = np.flip(np.where(self.reference == 1), axis=0)
        gt_mean = gt_indices.mean(axis=1)
        gt_cov = np.cov(gt_indices)