[R] Solving a tridiagonal system

[Will Harvey]

I need to find solutions to a tridiagonal system. By
this I mean a set of linear equations Ax = d where A
is a square matrix containing elements A[i,i-1],
A[i,i] and A[i,i+1] for i in 1:nrow, and zero
elsewhere. R is probably not the ideal way to do this,
but this is part of a larger problem that requires R.

In my application it is much easier (and much faster)
to generate the diagonal and off-diagonal elements of
A as vectors, i.e. a = A[i,i-1], b = A[i,i] and c =
A[i,i+1]. So I have three vectors that define A, along
with a solution vector d. The conventional method of
solving such systems is to use the so-called "Thomas
algorithm", see e.g. 
<http://www.enseeiht.fr/hmf/travaux/CD0001/travaux/optmfn/hi/01pa/hyb74/node24.html>.
This is very easy to code, but much more difficult to
"vectorize". Is anyone aware of a library that
contains a fast implementation of this algorithm?

Another alternative is to use backsolve. I can easily
eliminate the lower diagonal a, but I'm still left
with b and c, whereas backsolve requires a matrix.
Again, I can write a function to read b and c into a
matrix, but this requires loops, and is too slow. Is
there a vectorized way of doing it? Of course, the
diag command works for b, but what about c? In Octave,
diag allows for an offset, but R apparently does not.

I would appreciate any and all assistance you experts
can offer. Thanks in advance.

Will Harvey