CVS MatrixProduct

CVS is coordinate-wise in many ways. However, some operations are challenging. Performing CVS-to-CVS multiplication is possible but with a catch (at least as of this writing). CVS's own rule about keeping distinct non-zero values and aligning their lists of linear indexes creates a challenge when performing operations that can create zeroes or repeating values.

CVSExtensions' MatrixProduct() performs multiplication on two CVS instances (and, therefore, two matrices). In its current form, it performs loops within loops (total of 4: a loop starts a loop that starts a loop that starts a loop). The loops eventually fill a LinearIndexKeyedConcurrentSparseMatrix instance. The results are saved in the result sparse matrix, which is then compressed to CVS and returned.

The algorithm works. It's even inherently concurrent, which is fantastic! But even though LinearIndexKeyedConcurrentSparseMatrix theoretically uses 30-50% of memory that CoordinatesKeyedConcurrentSparseMatrix may use, the operation to multiply still requires a huge memory jump.

I have yet to investigate the possibility of performing the operation as a CVS-pure multi-step process. CVS-pure meaning there is no "cheating" or, specifically, no internal usage of DOK sparse matrix. Multi-step meaning that the process may include data clean-up, integrity checks and what not before returning the result CVS. I can imagine this operation to be non-concurrent though. I can also imagine it to be slow. But it can solve the memory space problem. Compromises!

So now for the code. In CVSExtensions.cs, add the following:

        public static CVS<P> MatrixProduct<T, M, P>(
            this CVS<T> cvs1, CVS<M> cvs2)
        {
            if (cvs1.Size.rows != cvs2.Size.columns)
            {
                throw new InvalidOperationException(
                    "Input matrix number of columns must match matrix number of rows.");
            }
            var r = new LinearIndexKeyedConcurrentSparseMatrix<P>(
                (cvs1.Size.rows, cvs2.Size.columns), cvs1.LinearIndexMode);
            Parallel.For(0, cvs1.Value.Count, (i) =>
            {
                T v1 = cvs1.Value[i];
                Parallel.ForEach(cvs1.LinearIndex[i], (li1) =>
                {
                    (int row, int column) c1 = MatrixCoordinates.ToCoordinates(
                        cvs1.Size, li1, cvs1.LinearIndexMode);
                    Parallel.For(0, cvs2.Value.Count, (j) =>
                    {
                        M v2 = cvs2.Value[j];
                        Parallel.ForEach(cvs2.LinearIndex[j], (li2) =>
                        {
                            (int row, int column) coordinates = (c1.row, 
                                MatrixCoordinates.ToCoordinates(
                                    cvs2.Size, li2, cvs2.LinearIndexMode).column);
                            int key = MatrixCoordinates.ToLinearIndex(
                                r.Size, coordinates, r.LinearIndexMode);
                            P value = r[coordinates.row, coordinates.column] 
                                + (v1 * (dynamic)v2);
                            if (value.Equals(default(P)))
                            {
                                r.Matrix.TryRemove(key, out P v);
                            }
                            else
                            {
                                r.SetEntry(key, value);
                            }
                        });
                    });
                });
            });
            return r.ToCVS();
        }

Have fun!