The whole point is that you don't generate a "cubemap" as output, you generate a set of SH coefficients. Basically you iterate over every texel of your input cubemap, compute the SH coefficients for that direction, multiply the coefficients by the cubemap value, multiply by a correction factor that accounts for cubemap->sphere distortion, and sum the results into a set of SH coefficients. Page 9 of this paper has some pseudo-code for the process.

The long answer:

To compute irradiance, you need to take your incoming radiance and integrate it with a cosine lobe oriented about the surface normal. With spherical harmonics you typically have your incoming radiance represented as a set of SH coefficients (this is what you're computing when you integrate your cubemap onto the SH basis functions), which means that it makes sense to also represent the cosine lobe with SH. If you do this, then computing the integral can be done using an SH convolution, which is essentially just a dot product of two sets of SH coefficients. This paper, which is also by Ravi Ramamoorthi, goes into the full details of how to represent your cosine lobe with SH coefficients. Basically you take the zonal harmonics coefficients for a cosine lobe oriented about the Z axis (This is just a set of constants, since it never changes), and then you rotate it so that it's now aligned with your normal direction. Once you've done that you perform the SH dot product with your incoming radiance, and you get irradiance as a result.

The short answer:

You compute a new set of coefficients from your normal direction, and perform a dot product with the SH coefficients that you got from integrating the cubemap. This is simple enough to do in a pixel shader, so you can just do it all at runtime per-pixel:

static const float Pi = 3.141592654f;
static const float CosineA0 = Pi;
static const float CosineA1 = (2.0f * Pi) / 3.0f;
static const float CosineA2 = Pi * 0.25f;

struct SH9
{
    float c[9];
};

struct SH9Color
{
    float3 c[9];
};

SH9 SHCosineLobe(in float3 dir)
{
    SH9 sh;

    // Band 0
    sh.c[0] = 0.282095f * CosineA0;

    // Band 1
    sh.c[1] = 0.488603f * dir.y * CosineA1;
    sh.c[2] = 0.488603f * dir.z * CosineA1;
    sh.c[3] = 0.488603f * dir.x * CosineA1;

    // Band 2
    sh.c[4] = 1.092548f * dir.x * dir.y * CosineA2;
    sh.c[5] = 1.092548f * dir.y * dir.z * CosineA2;
    sh.c[6] = 0.315392f * (3.0f * dir.z * dir.z - 1.0f) * CosineA2;
    sh.c[7] = 1.092548f * dir.x * dir.z * CosineA2;
    sh.c[8] = 0.546274f * (dir.x * dir.x - dir.y * dir.y) * CosineA2;

    return sh;
}

float3 ComputeSHIrradiance(in float3 normal, in SH9Color radiance)
{
    // Compute the cosine lobe in SH, oriented about the normal direction
    SH9 shCosine = SHCosineLobe(normal);

    // Compute the SH dot product to get irradiance
    float3 irradiance = 0.0f;
    for(uint i = 0; i < 9; ++i)
        irradiance += radiance.c[i] * shCosine.c[i];

    return irradiance;
}

float3 ComputeSHDiffuse(in float3 normal, in SH9Color radiance, in float3 diffuseAlbedo)
{
    // Diffuse BRDF is albedo / Pi
    return ComputeSHIrradiance(normal, radiance) * diffuseAlbedo * (1.0f / Pi);
}

I assembled and simplified this code from different places, so I apologize if there's a typo in there. But it should give you the general idea. Those cosine A0/A1/A2 terms are from Ravi's paper, and they're the zonal harmonics coefficients of a cosine lobe oriented around the Z axis. Also, you should notice how to compute the final diffuse value, I had to divide the result by Pi. This is because the actual diffuse BRDF is albedo / Pi, and if you forget that factor the result will be too bright. If you'd like, you can combine the 1 / Pi into the AO/A1/A2 terms, which simplifies nicely.

I had to change CosineA2 to be Pi * 0.25f, which corresponds with the 1, 2/3, 1/4 bands

Whoops, sorry about that! I corrected the code, in case anybody else looks at this thread later on.

Either way, I'm glad that you got it to work!