This statement has bothered me ever since I read it; after a quick read of George Hart’s “Multidimensonal Analysis”, I have been persuaded that a dimensionful Jacobian can be well defined and should also have a well-defined inverse. Any matrix that represents a set of dimensionally-consistent equations can be given a well-defined structure of its dimensions, as noted here and implemented by one person here. Certainly it is true that many algorithms may not be dimensionally consistent as currently implemented, but I am optimistic that many of these are not unresolvable issues, and maybe with wrappers that carry around the dimensional parts, we can perhaps even still use BLAS, etc.
On the other hand, discussion in Quantity currently incompatible with QuadGK · Issue #40 · SymbolicML/DynamicQuantities.jl · GitHub makes me inclined to believe that encoding dimensions in types hews closer to a mathematically rigorous definition of dimension.