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.
Ok first thanks for response and yes I need to dig into SH realms better.... but let's make it simple Stefano (well not so simple ): given a normal vector n (x,y,z) I can get second order SH coefficents c = (c0,c1,c2,c3)
c0 = sqrt(pi) / 2
c1 = -sqrt(pi/3) * y
c2 = sqrt(pi/3) * z
c3 = - sqrt(pi/3) * x
so for every VPL ( virtual point light ) I'll have a n vector and an associated SH coefficents c.
The question is how I can go back and compute reflected radiance from those second order SH coefficents?