One way to compute it, in the finite case, is by iterating over the members of the domain, and in each step compute R⁺ limited to the members considered thus far. This is known as Floyd's algorithm.

In A theorem on boolean matrices, Stephan Warshall's paper in the Journal of the ACM, 9(1):11-12, 1962, it was observed that when the incidence matrix of the relation is represented as a set of bit vectors, this algorithm is a matrix multiplication at each step. If this is a constant time operation, the total complexity of transitive closure becomes quadratic.

Warshall's algorithm can only be applied if the number of domain elements is limited to the size of your bit vectors; that is, the word size of your computer, or the vector size on specially built array processors that were popular for a while.

If the size of the domain is unknown in advance, but the relation is still dense, it's still a good idea to use the fast algorithms for matrix multiplication. When the relation is sparse, it becomes too wasteful to represent it as an incidence matrix, and one has to resort to Floyd's algorithm.