Implement A := c*A+B method in `MatrixHandler`

[pelesh]

Summary

Implement following operation in MatrixHandler class for all three suported hardware backends:

• $A := cA + I$ - Multiply all matrix element by a constant and add identity matrix to it.
• $A := cA + B$ - Multiply all matrix element by a constant and add another matrix to it.

Rationale

These methods are needed for interfacing with SUNDIALS, which requires matrix objects to implement following methods

• SUNMatScaleAdd,
• SUNMatScaleAddI.

Description

Add methods to the existing MatrixHandler implementations. Test on CPU, CUDA and HIP backends. Add unit tests to the MatrixHandlerTests class.

Additional information

Please keep in mind that this operation will be called repeatedly for matrices with the same sparsity pattern. Make sure the analysis part to compute the sparsity pattern of the resulting matrix is done once and available for reuse whenever the same operation is called.

[pelesh]

CC @shakedregev

[quantumsteve]

@pelesh Where might I find a test or an example of adding two matrices in Resolve?

[shakedregev]

You won't, but I implemented in HyKKT (cuda only) [here].(https://github.com/ORNL/hykkt/blob/develop/src/matrix_matrix_ops.cpp#L219)
This uses the raw pointers for the CSR arrays, but it should be easy enough to just pass a ReSolve CSR matrix and call the same functions by accessing those arrays through that class.

[pelesh]

Additional information: These operations are used when solving systems of differential equations of form $M\dot y + f(y) = 0$. An approximation at each time step is of form $\alpha M y + A y = r$, where $r$ is a residual vector. Solving for $y$ would require an operation