Since no one else has mentioned this, I guess I will. Why not use quaternions, rather than rotation matrices, to represent your rotations? Quaternions on the unit sphere and 3-d rotations are isomorphic, and quaternions don't require the redundant storage and calculation that 3x3 matrices do.

A quaternion may be thought of as an entity of the form [s,x], where s is a scalar and x is a 3-vector. Multiplication of quaternions is given by [s1,x1] * [s2,x2] = [s1 * s2 - x1 dot x2, s1*x2 + s2*x1 + x1 cross x2]. A unit quaternion is one that satisfies s*s + x dot x = 1. A unit quaternion may also be thought of as a rotation of angle 2 arccos s about the axis v. To rotate a vector v by a rotation quaternion q to get a vector w, use the formula [0,w] = inv(q) * [0,v] * q, where inv(q) * q = [1,0], and inv([s,v]) = [s,-v]. Or, if you prefer, form the equivalent rotation matrix

1 - 2 (x2*x2 + x3*x3)     2 (x1*x2 + s * x3)     2 (x1*x3 - s*x2)
2 (x1*x2 - s*x3)         1 - 2 (x1*x1 + x3*x3)   2 (x2*x3 + s*x1)
2 (x1*x3 + s*x2)          2 (x2*x3 - s*x1)      1 - 2 (x1*x1 + x2*x2)

and use that.

The basic algorithm, then, to display vectors V = [v1 ; v2 ; v3 ; ... ; vn] rotating by q every frame is

rot = [1,0]
do-forever
    R = rotation matrix associated with rot
    DISP = R * V
    display all vectors in DISP
    rot = rot * q
    norm = rot.s * rot.s + rot.x1 * rot.x1 + rot.x2 * rot.x2 + rot.x3 * rot.x3
    if (abs(norm - 1) > tolerance)
        norm = sqrt(norm)
        rot.s = rot.s/norm
        rot.x1 = rot.x1/norm
        rot.x2 = rot.x2/norm
        rot.x3 = rot.x3/norm
    endif
enddo

That's the general idea, anyway. For a less terse exposition, see Ken Shoemake, "Animating Rotation with Quaternion Curves", COMPUTER GRAPHICS Vol 19 No 3, pp. 245-254.

Incidentally, similar quaternion techniques can be used for 4-d rotations. I haven't been able to get a handle on higher dimensions, though.