RFC: Lyapunov / Sylvester solver via LAPack xTRSYL

As discussed in #5814, solver for the Lyapunov equation
``AX + XA' +C = 0``
and the Sylvester equation
``AX + XB + C = 0``

Function name and arguments is oriented at the Matlab neighbors.
lyapchol(A,B) for the cholesky factor of ``A*X + X*A' + B*B' = 0`` is still missing.


Yes, I was also looking for the right way to compute A*X*A' for lyap()