Direct way of computing clockwise angle between 2 vectors

2D case

Just like the dot product is proportional to the cosine of the angle, the determinant is proprortional to its sine. So you can compute the angle like this:

dot = x1*x2 + y1*y2      # dot product between [x1, y1] and [x2, y2]
det = x1*y2 - y1*x2      # determinant
angle = atan2(det, dot)  # atan2(y, x) or atan2(sin, cos)

The orientation of this angle matches that of the coordinate system. In a left-handed coordinate system, i.e. x pointing right and y down as is common for computer graphics, this will mean you get a positive sign for clockwise angles. If the orientation of the coordinate system is mathematical with y up, you get counter-clockwise angles as is the convention in mathematics. Changing the order of the inputs will change the sign, so if you are unhappy with the signs just swap the inputs.

3D case

In 3D, two arbitrarily placed vectors define their own axis of rotation, perpendicular to both. That axis of rotation does not come with a fixed orientation, which means that you cannot uniquely fix the direction of the angle of rotation either. One common convention is to let angles be always positive, and to orient the axis in such a way that it fits a positive angle. In this case, the dot product of the normalized vectors is enough to compute angles.

dot = x1*x2 + y1*y2 + z1*z2    #between [x1, y1, z1] and [x2, y2, z2]
lenSq1 = x1*x1 + y1*y1 + z1*z1
lenSq2 = x2*x2 + y2*y2 + z2*z2
angle = acos(dot/sqrt(lenSq1 * lenSq2))

Edit: Note that some comments and alternate answers advise against the use of acos for numeric reasons, in particular if the angles to be measured are small.

Plane embedded in 3D

One special case is the case where your vectors are not placed arbitrarily, but lie within a plane with a known normal vector n. Then the axis of rotation will be in direction n as well, and the orientation of n will fix an orientation for that axis. In this case, you can adapt the 2D computation above, including n into the determinant to make its size 3×3.

dot = x1*x2 + y1*y2 + z1*z2
det = x1*y2*zn + x2*yn*z1 + xn*y1*z2 - z1*y2*xn - z2*yn*x1 - zn*y1*x2
angle = atan2(det, dot)

One condition for this to work is that the normal vector n has unit length. If not, you'll have to normalize it.

As triple product

This determinant could also be expressed as the triple product, as @Excrubulent pointed out in a suggested edit.

det = n · (v1 × v2)

This might be easier to implement in some APIs, and gives a different perspective on what's going on here: The cross product is proportional to the sine of the angle, and will lie perpendicular to the plane, hence be a multiple of n. The dot product will therefore basically measure the length of that vector, but with the correct sign attached to it.

To compute angle you just need to call atan2(v1.s_cross(v2), v1.dot(v2)) for 2D case. Where s_cross is scalar analogue of cross production (signed area of parallelogram). For 2D case that would be wedge production. For 3D case you need to define clockwise rotation because from one side of plane clockwise is one direction, from other side of plane is another direction =)

Edit: this is counter clockwise angle, clockwise angle is just opposite