How do I parallelize this matrix by vector multiplication?

I need to parallelize this code. As you can see, the q array contains values that are summed at every loop.

/// <summary>
/// q = A * d, with A stored in CSR format.
/// </summary>
public static double[] Multiplication(int order, double[] d, List<double> A, List<int> J, int[] I)
{ 
    double[] q = new double[order];

    for (int i = 0; i < order; i++)
    {                    
       int indexFirst = I[i];

       q[i] += A[indexFirst] * d[J[indexFirst]];

       for (int j = indexFirst + 1; j < I[i + 1]; j++)
       {
          int col = J[j];
          double a = A[j];

          q[i] += a * d[col];

          q[col] += a * d[i];
       }
   }

   return q;
}

I was able to parallelize the code in this way (splitting the problem in a CPU count number of sections and summing up resulting q arrays in the end) but it seems overcomplicated. Any better idea?

Parallel.For(0, cpuCount, delegate(int i)
{
    MultiplySection(startIndices[i], endIndices[i], d, A, J, I, qSection[i]);
});

double[] q = new double[order];

for (int j = 0; j < cpuCount; j++)
{
    for (int i = startIndices[j]; i < order; i++)
    {
        q[i] += qSection[j][i];
    }
}

private static void MultiplySection(int startIndex, int endIndex, double[] d, List<double> A, List<int> J, int[] I, double[] q)
{
    for (int i = startIndex; i < endIndex; i++)
    {
        int indexFirst = I[i];

        q[i] += A[indexFirst]*d[J[indexFirst]];

        for (int j = indexFirst + 1; j < I[i + 1]; j++)
        {
            int col = J[j];
            double a = A[j];

            q[i] += a*d[col];

            q[col] += a*d[i];
        }
    }
}

Your current implementation is coarse grained, so if a CPU core is busy with something else you risk having poor utilization at the end of your process.

I would at least try using the built in partitioning by Parallel.For. This should give a much fine grained concurrency that might help with efficiency. You can also use localInit/localFinally methods to create and merge your arrays:

public static double[] Multiplication(int order, double[] d, List<double> A, List<int> J, int[] I)
{
    double[] qTotal = new double[order];

    Parallel.For(0,
        order,
        () => new double[order], // LocalInit
        (i, loopState, q) => // Main Body
        {
            int indexFirst = I[i];

            q[i] += A[indexFirst] * d[J[indexFirst]];

            for (int j = indexFirst + 1; j < I[i + 1]; j++)
            {
                int col = J[j];
                double a = A[j];

                q[i] += a * d[col];

                q[col] += a * d[i];
            }
            return q;
        },
        q => // Local Finally 
        {
            lock (qTotal)
            {
                for (int i = 0; i < q.Length; i++)
                {
                    qTotal[i] += q[i];
                }
            }
        });
    return qTotal;
}

If the amount of work inside the loop body is small you might still benefit from batching, but you might want to use a fixed batch size instead of basing it of the number of processors. This should help a bit with various types of overheads. You will likely need to experiment a bit to find your optimal batch size.

The following code, from an answer to Parallel vs Serial results aren't understandable results in 9x faster performance in a 6-core hyper-threaded machine. That's 1.5 the number of cores. It could be improved perhaps with SIMD operations for another 2x-4x speedup, depending on the CPU's SSE/AVX register size

static double[,] MatrixProductPVect(double[,] matrixA,
    double[,] matrixB)
{
    int aRows = matrixA.GetLength(0); int aCols = matrixA.GetLength(1);
    int bRows = matrixB.GetLength(0); int bCols = matrixB.GetLength(1);
    if (aCols != bRows)
        throw new Exception("xxx");

    double[,] result = new double[aRows, bCols];

    Parallel.For(0, aRows, i =>
    {
        double[] row = new double[aCols];
        System.Buffer.BlockCopy(matrixA, i*aCols*sizeof(double), row, 0, aCols*sizeof(double));
        // each row of A
        for (int j = 0; j < bCols; ++j) // each col of B
        {
            double temp=0;
            
            for (int k = 0; k < aCols; ++k) // could use k < bRows
                temp += row[k] * matrixB[k, j];
            result[i, j] = temp;
        }
    });
    return result;
}

The main performance killers in parallel processing is synchronization and non-local memory access.

• It's faster to access values that are already in the CPU's cache than going to RAM.
• It's also faster to have each core work with is own partition of the
  data than trying to synchronize access to the same data.
• With the cached data, it's faster to access values along a row instead of
  jumping N items to get to the next column on the next row.
• Finally, when calculating totals it's faster to calculate a local total,
  that's stored in a CPU register, than trying to increment a shared
  variable.

This code takes advantage of these rules in the following way:

• Parallel.For will partition the rows among all the cores, with each core accessing its own consecutive rows
• Processing works along rows, which means the data are close together and usually already in cache
• System.Buffer.BlockCopy makes a local copy of each row so the cores don't have to synchronize access to RAM.
• The temp total is calculated locally, inside the column loop.

The multiplication code was borrowed from this article.

The answer uses BenchmarkDotNet to compare multiple multiplication versions, from sequential, to naive parallel, to using local row copies

BenchmarkDotNet=v0.13.1, OS=Windows 10.0.22621
Intel Core i7-10850H CPU 2.70GHz, 1 CPU, 12 logical and 6 physical cores
.NET SDK=6.0.302
  [Host]     : .NET 6.0.7 (6.0.722.32202), X64 RyuJIT
  DefaultJob : .NET 6.0.7 (6.0.722.32202), X64 RyuJIT


           Method |      Mean |    Error |   StdDev | Ratio |
----------------- |----------:|---------:|---------:|------:|
       SerialMult | 288.63 ms | 5.736 ms | 8.585 ms |  1.00 |
     ParallelMult |  75.92 ms | 1.311 ms | 1.162 ms |  0.26 |
 ParallelMultTemp |  49.63 ms | 0.873 ms | 0.817 ms |  0.17 |
 ParallelMultVect |  33.04 ms | 0.588 ms | 0.521 ms |  0.11 |