Add rules for solutions to Sylvester and Lyapunov equations #384
Conversation
FiniteDifferences only consistently works during overflow, but locally it works sometimes, which is enough to know that the scaling is correctly handled.
Codecov Report
@@ Coverage Diff @@
## master #384 +/- ##
==========================================
+ Coverage 97.72% 97.80% +0.07%
==========================================
Files 19 20 +1
Lines 1496 1547 +51
==========================================
+ Hits 1462 1513 +51
Misses 34 34
Continue to review full report at Codecov.
|
| trans = T <: Complex ? 'C' : 'T' | ||
| ∂D, scale2 = LAPACK.trsyl!(trans, trans, RA, RB, ∂Y) | ||
| ∂Z = rmul!(QA * (∂D * QB'), -inv(scale2)) | ||
| return NO_FIELDS, @thunk(∂Z * Ω'), @thunk(Ω' * ∂Z), @thunk(∂Z * inv(scale)) |
There was a problem hiding this comment.
I wonder if we should have a shorthand for inplaceable thunking multiplication.
Maybe @thunk should be smart in that way?
Not for this PR but just something to consider
There was a problem hiding this comment.
That would be convenient. I haven't really put any effort into writing InplaceableThunks, since nothing supports them yet. Simple cases like this would be easy to transform but maybe could break? e.g. for the last thunk here, the macro would not know that scale is a number not a matrix. And is the case that we thunk just the multiplication common enough to make the code complexity worth it?
Co-authored-by: Lyndon White <[email protected]>
This PR adds rules for the following functions:
fruleforLAPACK.trsyl!(solution to the Sylvester equation for (quasi)-triangular inputs)frule/rruleforsylvester(solution to the Sylvester equation)frule/rruleforlyap(solution to the Lyapunov equation)All of the rules can be written in terms of the primal function itself. However,
sylvesterandlyapboth use a Schur decomposition to solve an easier equation usingLAPACK.trsyl!, so these rules then reuse that Schur decomposition. As a result, the included rule forlyapis much more efficient than the one in Zygote forStridedMatrixinputs.