SUNMatScaleAdd_Sparse has time complexity O(M*N)

[phannebohm]

As the title says, adding two sparse matrices scales badly. I believe the problem are loops like the sparsity pattern overlap check here:

sundials/src/sunmatrix/sparse/sunmatrix_sparse.c

Lines 827 to 832 in 7c31824

	/* determine if A already contains the sparsity pattern of B */
	newvals = 0;
	for (j=0; j<N; j++) {
	/* clear work array */
	for (i=0; i<M; i++) w[i] = 0;

Just filling the work array w with zeros takes M*N operations. This happens in a few places inside this function with w and x. If I replace the for loops with calls to memset I save a significant amount of computation time (maybe 75% on my machine, don't trust that number though). However it still scales badly. I would expect addition to scale like O(NNZ) instead.

On a side note, for my particular use case (OpenModelica/OpenModelica#10178) I can take a modified version of SUNMatScaleAddI_Sparse, I call it SUNMatScaleIAdd_Sparse where the result is $A + c I$. This scales linearly as expected. Perhaps this function could be useful for other people.

[balos1]

Thanks for noting this. Indeed this is because the current implementation support a very general case. Perhaps we can add an option (controlled with a Set function) to use less general, but faster implementations of ScaleAdd and ScaleAddI.

[phannebohm]

I think I have a general solution that avoids these scaling issues, at least for ScaleAddI. I can open a PR for it and see if it can also be ported to ScaleAdd.

[balos1]

Closed by #257.