Parsing to wrap Level 3 BLAS routines syr2k/her2k

[jiahao]

PR #4813 provides a BLAS wrapper to the symmetric rank 2k update function syr2k, which provides the updates:

C := alpha*A*B' + alpha*B*A' + beta*C,
C := alpha*A'*B + alpha*B'*A + beta*C,

Is it worth parsing Julia expressions of this form and turning them into calls to syr2k/her2k?

[StefanKarpinski]

Doing this in the parser is not really the direction we should be heading here, IMO, but it's not entirely clear how else to approach it.

[kmsquire]

It would be nice if, at some point, the parser provided hooks to allow user code to analyze and rewrite expressions like this.

Kevin

[StefanKarpinski]

Yes. We've discussed this forever and haven't come up with anything satisfactory yet. It's really hard to make something like this not super brittle.

[JeffBezanson]

The only thing that comes to mind is something like Haskell rewrite rules.

[StefanKarpinski]

Honestly, I'm not even sure we should do this. It's getting into "sufficiently clever compiler" territory. Having the syr2k routines there for those who know about them and want to call them is already pretty good. This fits with our whole philosophy of doing something simple, understandable and reasonably fast by default but allowing experts to easily tune their codes to get even more performance with a bit of work. Looking for opportunities to use calls like syr2k and her2k seems like the kind of things human experts are best cut out for.

[JeffBezanson]

Yeah, it leaves me a bit cold too.

[ViralBShah]

I wouldn't mind trying out rewrite rules and it would be better if they can be expressed in julia rather than being part of the parser. Rewrite rules combined with type inference could make it possible for us to efficiently deal with a lot of array expressions.

rewrite
   C::Matrix = A::Matrix'*A::Matrix
as
   someblascall(A,C)
end

The other big thing is reusing memory in something like eigs or anything iterative, where repeated matrix-vector products happen.

[JeffBezanson]

That approach is certainly better than mucking with the parser.

[jiahao]

@ViralBShah's idea was sort of what I had in mind. I take it that this is implemented on top of the base Julia parser?

[JeffBezanson]

Yes, there would be no parser changes. These would be patterns that the compiler can match and apply if it is able to. They are advisory and for performance only.

[jiahao]

Yeah, that's the clearly the way to go then.

[andreasnoack]

@JeffBezanson do you think anything like #4950 (comment) will ever happen? If not, we should probably close this issue.

[KristofferC]

Adding this to the parser seems quite unlikely due to how dynamic Julia has become.