bertiniteam
computational inefficiency
I found some pessimizations in explicit_predictors.hpp, in FullStep(). A temp vector is set to 0 when it doesn't need to be. It can just be set to the first of them, and then the loops can run from 1 to bla instead of from 0 to bla. this is probably a fairly important improvement to make.
I think there's a mathematical error in this code, from explicit_predictors.hpp 741:
// TODO this random vector should not be made fresh every time. especiallyif the numeric type is mpfr_complex!
Vec<ComplexType> randy = RandomOfUnits<ComplexType>(numVariables_);
Vec<ComplexType> temp_soln = LUref.solve(randy);
norm_J = NumErrorT(dhdxref.norm());
norm_J_inverse = NumErrorT(temp_soln.norm());
No, this is actually correct. But doing that extra LU is potentially VERY expensive. do we really need this? Is there a better way?
I think that the ODE predictor information, the numbers for the Butcher tables, etc, should be written in at compile time. Refactor to use a typetraits system?
EvalRHS() overwrites the columns of K. It does not add to them. Is it pointless to set K to 0 before calling this function?
In FullStep, we tell EvalRHS the stage. Does EvalRHS already have the logic for "on stage 0, just overwrite; other stages add"?
Why do stages later than 0 NOT re-use the LU decomposition member variable? Are we saving the one for stage 0 for something special?
S.SetAndReset<ComplexType>(space, time);
Mat<ComplexType>& dhdxtempref = std::get< Mat<ComplexType> >(dh_dx_temp_);
S.JacobianInPlace(dhdxtempref);
auto LU = dhdxtempref.lu(); //<<<<---- here!!!!
line 886, explicit_predictors in EvalRHS(
Inside all stages, EvalRHS overwrites the stage-th column of K. They both have this line of code:
K.col(stage) = LUref.solve(-dhdtref);
See https://eigen.tuxfamily.org/dox/group__InplaceDecomposition.html, where they describe how to do in-place LU decompositions. This could really be nice.
Todo: rewrite the LU stuff in explicit_predictors.hpp to use In-place. Warning: if there are later things wanting to use the temp matrices and they depend on the previous state of the matrices, using LU will make that code erroneous. Need to be better about describing when a thing is expected to remain something special.
Actually, I think @jbcolli2 did a decent job of making variables for temps that have some purpose beyond their first use, by adding _0_ to their names. So maybe a naming convention? I'd prefer something expressed at compile time, but hey, just names are a good place to start.
in LUPartialPivotDecompositionSuccessful, why doesn't the (0,0) entry come first?