I am not familiar with RSMs but, judging from the code you posted, I believe SH suit your needs.
Remember that SH represent a signal defined over a spherical domain. In this case the signal is the total incoming light from all the secondary sources. A SH probe has an associated position in world space. For a <light, probe> pair, the signal is (cosThetaI * lR * Flux) in the direction defined by R (normalized). The cosThetaJ term depends on the receiver's normal and will be taken into account later at runtime when shading a surface.
for each probe
SH = 0
for each light
signal = R / | R | * (cosThetaI * lR * Flux)
SH += signalToSH(signal)
signalToSH is a function that encodes a directional light in the SH basis. See D3DXSHEvalDirectionalLight for an implementation.
At runtime, when shading a pixel, pick the probe spatially closest to it. The pixel has a normal associated. Assuming that the BRDF is purely diffuse, you have to convolve the cosine lobe around the pixel's normal with the signal encoded in the probe, achieved with a SH rotation and dot product. See the paper http://www.cs.berkeley.edu/~ravir/papers/envmap/envmap.pdf
I suggest you read more material on SH in order to have a better understanding of their properties and use cases. As for myself, I will read the RSM paper asap and edit my reply if I've made wrong assumptions.
Please let me know if you have more questions.