In optimization algorithms, at times it is useful to factorise a large sparse matrix ONCE ONLY and use the factors in different ways. For example, we can do some computation using the factors (eg. block parts of the L and P in PAQ=LU are used for null space methods in image processing and rendering). This then can be used in many computations without re-factoring the original matrix A. Unfortunately,  matrixLU() on a type of Eigen::SparseLU solver returns a reference to a type which only allows a triangular solve -- matrixL().solveInPlace(b);.

It would be great to implement other operations on this.

Thanks.