Strictly speaking, neither r nor R^2 can be computed with those vectors (so yes, you could consider it a bug).   This is because no line can be fitted with finite values for $c0 and $c1.  And when there's no line, there are no residuals, etc.

When you change the values slightly to

2,2,2,2.0001
2,2,2,1
[download]

a line can be fitted (so you also get correct values on that site), but the fitted line parameters are already rather large:

my $c0 = 20001.9999999578;    # as computed by gsl_fit_linear()
my $c1 = -9999.9999999789;
say $c0 + $c1 * $_ for 2, 2.0001;

__END__
2
1
[download]

The smaller you make the deviation from 2 for the last value of vector x, the larger the fitted parameters become, eventually approaching delicately balancing +/- "infinities".  And when all values are exactly 2, no fit can be computed any longer...   (try 2.0000001 and 2.00000001 on the linked site, and you'll already get nonsensical values (correct values should always be r = -1, r^2 = 1) — which means numerical precision is rather low).