First of all, for a rotation matrix the two known columns have to be an orthonormal pair, meaning that there are conditions on the initial six variables. If these variables are relabeled as two 3x1 column vectors c1,c2 then
norm(c1) = norm(c2) = 1; dot(c1,c2) = 0
Let's assume that's true.
From your eqns, at first glance it might seem that there could be eight solutions, since you can have +-1 for each of r13,r23,r33. However, the rows have to be orthogonal, so e.g.
r11*r21 + r12*r22 + r13*r23 = 0
If you have a solution, then changing one of r13 or r23 will change the sign of the last term and will not provide a solution. You can change both signs or neither, but not just one. Similarly for the other pairs of rows. The net result is that there are two solutions. One solution is third column equaling the cross product of the first two columns (thus being perpendicular to both of them), and the other solution is simply that solution times an overall minus sign.
However, changing the overall sign in the third column also changes the sign of the determinant of the matrix. As Bruno pointed out, a rotation matrix must have determinant +1. So the only solution is
in which case
has determinant +1.
Note that for this to come out correctly there is a sign convention present (or one might say, hidden) in the cross function.