diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp index 57f68a8395d..e83c84c9dec 100644 --- a/Core/include/Acts/Propagator/EigenStepper.ipp +++ b/Core/include/Acts/Propagator/EigenStepper.ipp @@ -233,8 +233,41 @@ Acts::Result Acts::EigenStepper::step( return EigenStepperError::StepInvalid; } - // for moment, only update the transport part - state.stepping.jacTransport = D * state.stepping.jacTransport; + // See the documentation of Acts::blockedMult for a description of blocked + // matrix multiplication. However, we can go one further. Let's assume that + // some of these sub-matrices are zero matrices 0₈ and identity matrices + // I₈, namely: + // + // D₁₁ = I₈, J₁₁ = I₈, D₂₁ = 0₈, J₂₁ = 0₈ + // + // Which gives: + // + // K₁₁ = I₈ * I₈ + D₁₂ * 0₈ = I₈ + // K₁₂ = I₈ * J₁₂ + D₁₂ * J₂₂ = J₁₂ + D₁₂ * J₂₂ + // K₂₁ = 0₈ * I₈ + D₂₂ * 0₈ = 0₈ + // K₂₂ = 0₈ * J₁₂ + D₂₂ * J₂₂ = D₂₂ * J₂₂ + // + // Furthermore, we're constructing K in place of J, and since + // K₁₁ = I₈ = J₁₁ and K₂₁ = 0₈ = D₂₁, we don't actually need to touch those + // sub-matrices at all! + if ((D.topLeftCorner<4, 4>().isIdentity()) && + (D.bottomLeftCorner<4, 4>().isZero()) && + (state.stepping.jacTransport.template topLeftCorner<4, 4>() + .isIdentity()) && + (state.stepping.jacTransport.template bottomLeftCorner<4, 4>() + .isZero())) { + state.stepping.jacTransport.template topRightCorner<4, 4>() += + D.topRightCorner<4, 4>() * + state.stepping.jacTransport.template bottomRightCorner<4, 4>(); + state.stepping.jacTransport.template bottomRightCorner<4, 4>() = + (D.bottomRightCorner<4, 4>() * + state.stepping.jacTransport.template bottomRightCorner<4, 4>()) + .eval(); + } else { + // For safety purposes, we provide a full matrix multiplication as a + // backup strategy. + state.stepping.jacTransport = D * state.stepping.jacTransport; + } } else { if (!state.stepping.extension.finalize(state, *this, h)) { return EigenStepperError::StepInvalid;