The current C code reduces the number of parallel threads to stabilize convergence when the positivity constraint is turned off. If we implement the relaxation factor for voxel update, then there will be no need to reduce the parallel threads. I tested the C code with relaxation factor set to be 0.8 on NERSC cori cluster. That stabilizes the algorithm convergence.