support for `Base.expm` (need advice implementing)

It would be very nice to have support for AD in Base.expm. I would like to contribute to it, but since I am unfamiliar with the AD architecture in Julia I would appreciate some advice.

For Frechet derivatives, the recommended approach is Al-Mohy and Higham (2009). But I am unsure that this is the best approach for AD.

I came across an article by Mike Giles on AD and matrix functions. Section 2.3.5 in this paper suggests that if one uses the squaring and scaling algorithm, AD calculations just amount to the solution of linear systems and matrix powers, both of which go through nicely.

Base.expm! uses Julia functions, except for two LAPACK ones, gebal! and gesv!, for balancing and solving linear systems, respectively. I am under the impression that AD already goes through linear systems, so it could be replaced. So all one would need to do is implement gebal in Julia, and replace gesv! by \, perhaps at the cost of some extra allocation.

Am I on the right track with this?

Start by reading Base.expm! and the Julia AD handling for linear systems, especially the gebal! and gesv! calls mentioned in the issue. Compare the squaring and scaling approach with the cited Frechet-derivative and Giles references. Done means establishing a decided implementation path and adding AD support for Base.expm.