Dear developers,

I need to solve the following linear system:

AX+XA^dag=M

where A and M are known "n x n" matricies and X is the unknown "n x n" matrix

There are two possible strategy to solve such problem:

a) Recast it as AA vec[X] = vec[M]

where

AA = (1 x A) + (A* x 1),

i.e. an "n^2 x n^2" matrix and vec[X] and vec[M] are two "n^2 x 1" vectors

b) Diagonalize A and rewrite the problem in the rotated basis set

a_i = eigenalues of A

Y=U^dag X U new unknown matrix

M'=U^dag M U

Y_{ij}(a_i+a_j) = M_{ij}

However they are both problematic for me, because:

a) works fine, but it is too slow, needing to deal with an "n^2 x n^2" matrix

b) is not numerically sensitive, since the A matrix is of the form 1+x, with |x|<<1, and doing the diagonalization there is a lot of noise

I was thus wondering if it exist, or if it would be possible to have, a lapack subroutine which solves the commutator linear system of equations in a faster way than a) without needing to diagonalize (or invert) the A (or A^dag) matrix.

Kind regards,

D.