Diagonals of a Matrix

Hi!

Working with @JPelamatti and @vchabri on the development of a new feature for HSIC indices, we have a small issue at porting some code from numpy matrices to OpenTURNS matrices.
The following code removes the diagonal of the covariance matrix:

import numpy as np
covarianceMatrix = covarianceModel.discretize(sample)
covarianceMatrix = covarianceMatrix - np.diag(np.diag(covarianceMatrix))

What we would need is:

# Use case 1
A = ot.Matrix([[1 2 3], [4 5 6]])
d = A.getDiagonal() # d = [1 5] > Point

# Use case 2
d = ot.Point([1.0, 5.0])
A = ot.IdentityMatrix(2)
A.setDiagonal(d) # fill_diagonal in numpy

# Use case 3
A = ot.Matrix([[1 2 3], [4 5 6], [7, 8, 9]])
d = ot.Point(3)
A.setDiagonal(d)

# Use case 4
A = ot.Matrix([[1 2 3], [4 5 6], [7, 8, 9]])
d = A.getDiagonal() # d = [1 5, 9] > Point
D = ot.DiagonalMatrix(d)
A -= D

But these methods do not exist.
Do you agree that these would be useful / necessary?
Regards,
Michaƫl

Sure, why not? You can tell it would be useful in some places of the C++ code:

If @JPelamatti wants to try his hand at C++ programming, this could be a good exercise.