Okay, the thing I thought about yesterday would have been okay for cylinders whose central axes intersect. However, this is not the case in general, and the larger the offset, the greater the error would have been.
Since I'm still not confident I understand your problem, could you draw me a formulation of the problem, tube-tube case first?
Meanwhile, for line-line intersection (I'm not sure whether you have already implemented this), the following should work:
P1 and P2 are points on the corresponding lines measured from GLOBAL origin, and thus vectors P1 and P2 are straightforward to define from the global origin. Additionally, you'll need the direction vectors for the both lines, be them n1 and n2 correspondingly. Normalize these vectors to unity.
First, compute n1 x n2. If the result is zero vector, the lines are parallel to each other and can only intersect if P1 = P2, otherwise no intersection is possible. If they do intersect, they'll intersect the whole distance they are defined at.
Second step, since we now know that n1 x n2 != 0. We need to rule out the possibility that the three dimensional lines are skew, and thus could not possibly intersect. Compute (P1 - P2) * (n1 x n2). If the result is zero, intersection is possible. If the result non-zero, the lines are skew and cannot intersect.
Third step is the traditional solving free parameters r and s for P1 + n1r = P2 + n2s [1]
You could open up this to three equations for each coordinate axis, and solve for parameters s and r and substitute them to the last remaining equation to check whether it remains valid. The problem with this approach is that it gives rise for many possibilities of dividing by zero, or that some component of either normal vector is zero. Not good.
Instead, let's do the vector approach. Take a cross-product of [1] with n2. This eliminates one free parameter completely, one of the great things in vector algebra.
Result is n2 x (P1 + n1r) = n2 x (P2 + n2s), and because cross-product is distributive over addition, this can be written as
n2 x P1 - n2 x n1r = n2 x P2 + n2 x n2s, where n2 x n2 cancels. Re-arrange the remaining terms over right side, and leave r for the left, which leads to
n2 x n1 r = n2 x P2 - n2 x P1 = n2 x (P2 - P1). It sure would be nice to be able to leave r alone, but because this is a vector equation, we cannot divide the left side just with a vector. Instead, it pays off to multiply the equation with the dot-product of n2 x n1.
This results in (n2 x n1) * (n2 x n1) r = (n2 x n1) * [n2 x (P2 - P1)], and since the dot product of the vector itself is a scalar, it is now possible to divide the whole set with r's coefficient, leading to
r = [(n2 x n1) * [n2 x (P2 - P1)] / | n2 x n1 |2
Substitute r to the P1 + n1r and you'll get the intersection coordinates. I'm not sure whether denominator normalizes to unity, so it is better to leave it that way.
Now, and this is a very big NOW, you cannot use these generally to determine the intersection points of the cylinders, especially if the cylinders have large base radii. Why? Because while skew lines do not intersect, skew cylinders can and will intersect when you least expect it, and this is the complicative issue.
The reason for the order of the three steps is the following. Since there is a multitude of pieces in the geometry and many of them needed to be tested, the intersection test should terminate as soon as possible if no intersection becomes apparent, and only go to heavy calculation as a last step - basically, we are sort of assuming that the intersection does not happen and this tends to be the case.