Solving a linear system in matrix notation I'm trying to implement the paper "Conformation Constraints for Efficient Viscoelastic Fluid Simulation": https://www.gmrv.es/Publications/2017/BGAO17/SA2017_BGAO.pdf
In implementing eqn. 6, they solve the following linear system for $Q$:
$$
\frac{Q - Q_0}{\Delta t} = Q \nabla \vec{u} + (\nabla \vec{u})^T Q - \frac{1}{\tau} (Q - I)
$$
$Q_0$, $\nabla \vec{u}$, $\Delta t$ and $\tau$ are all knowns. $Q$ is a symmetric 3x3 matrix.
In the paper, they describe reformulating $Q$ as a 6-component vector: $$q = [Q_{xx},Q_{yy},Q_{zz},Q_{xy},Q_{xz},Q_{yz}]^T$$
With $I$ becoming the vector: $$I = \bar{q} = [1,1,1,0,0,0]^T$$
To be able to solve this in practice I think I need to get it in the form $Aq = b$. 
I'd love some help figuring out how to reformulate the above equation into a matrix-vector product that I can use a linear solver on, as I can't seem to just naïvely rearrange terms.
 A: The authors are using Mandel-Voigt notation, but 
you can use ordinary vectorization, e.g.
$${\rm vec}(ABC) = (C^T\otimes A)\,{\rm vec}(B)$$
to flatten the matrices into vectors
For typing convenience, define the variables
$$\eqalign{
M &= \nabla u,\quad
B = -\bigg(\frac{I}{\tau} + \frac{Q_0}{\Delta t}\bigg),\quad
\lambda = -\bigg(\frac{\tau+\Delta t}{\tau\Delta t}\bigg) 
\\
}$$
Write the equation in terms of these variables. Then vectorize it.
$$\eqalign{
B &= QM + M^TQ + \lambda Q \\
B &= IQM + M^TQI + \lambda IQI \\
{\rm vec}(B) &= (M^T\otimes I + I\otimes M^T + \lambda I\otimes I)\,{\rm vec}(Q) \\
b &= Aq \\
}$$
which is in the form that you wanted.
For symmetric matrices, another standard technique is half-vectorization. Duplication and elimination matrices can be used to convert between vec and vech representations.
$$\eqalign{
 {\rm vech}(B) &= L_n\;{\rm vec}(B) &\implies b_h = L_nb \\
 {\rm vec}(Q) &= D_n\;{\rm vech}(Q) &\implies q = D_nq_h \\
}$$
These matrices can be used to write a half-vec version of the above
$$b_h = (L_nAD_n)\,q_h$$
