Solution 1:

To settle this question: one can use the Rodrigues rotation formula to construct the rotation matrix that rotates by an angle $\varphi$ about the unit vector $\mathbf{\hat u}=\langle u_x,u_y,u_z\rangle$ (and if your vector is not a unit vector, normalization does the trick). Letting

$$\mathbf W=\begin{pmatrix}0&-u_z&u_y\\u_z&0&-u_x\\-u_y&u_x&0\end{pmatrix}$$

the Rodrigues rotation matrix is constructed as

$$\mathbf I+\left(\sin\,\varphi\right)\mathbf W+\left(2\sin^2\frac{\varphi}{2}\right)\mathbf W^2$$

where $\mathbf I$ is an identity matrix.

(For those of an advanced bent, one constructs $\mathbf W$ from $\mathbf{\hat u}$ through a premultiplication with the Levi-Civita tensor.)

Conventionally, the scalar multiplying the $\mathbf W^2$ term above is written as $1-\cos\,\varphi$, but this version is more prone to subtractive cancellation when $\varphi$ is near $2k\pi$ ($k$ is an integer), so the expression with the sine is more numerically sound.