Compare computation speed of Matrix:: versus base::

[ ] compare diag(C1[, i]) to Diagonal(x = C1[, i])to Matrix(diag(C1[, i]), sparse = TRUE) since Diagonal keeps 0's
[ ] compare rhs = t(X) %*% diag(C1[, i] + 1) %*% P[, i] to rhs = t(X) %*% CP[, i]
[ ] compare Lambda = Matrix::Diagonal(f, x = lambda) to Lambda = diag(x = lambda, f, f)
[ ] compare X = Matrix(...) to X = matrix(...)

andland / implicitcf