Euclidean distance between data sets in R using rdist from "fields" package

Question

I am using rdist to compute Euclidean distances between a matrix and itself:

> m = matrix(c(1,1,1,2,2,2,3,4,3),nrow=3, ncol=3)
> m
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    4
[3,]    1    2    3
library(fields)
> rdist(m)
      [,1]  [,2]  [,3]
[1,] 1e-10 1e+00 1e-10
[2,] 1e+00 1e-10 1e+00
[3,] 1e-10 1e+00 1e-10

What confuses me is that I think it should have 0s on the diagonal (surely the distance of a vector to itself is 0?), and for the same reason it should have 0s where it compares the first and third row. The value that I see instead (1e-10) looks way to big to be numerical noise. What's going wrong?

EDIT: rdist is from the package fields.

Answer 1

First of all 1e-10 is simply 1*10^-10 which is 0.0000000001, so numericaly very close to 0 (as it is a result of square rooting, so the actual error in the computation is of row of magnitude 1e-20). Is it "too big"? Well, library is written in fortran, and is focused on speed, so it is quite acceptable. If you analyze the exact code, you will find out how it is computed:

# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"rdist" <- function(x1, x2) {
    if (!is.matrix(x1)) 
        x1 <- as.matrix(x1)
    if (missing(x2)) 
        x2 <- x1
    if (!is.matrix(x2)) 
        x2 <- as.matrix(x2)
    d <- ncol(x1)
    n1 <- nrow(x1)
    n2 <- nrow(x2)
    par <- c(1/2, 0)
    temp <- .Fortran("radbas", nd = as.integer(d), x1 = as.double(x1), 
        n1 = as.integer(n1), x2 = as.double(x2), n2 = as.integer(n2), 
        par = as.double(par), k = as.double(rep(0, n1 * n2)))$k
    return(matrix(temp, ncol = n2, nrow = n1))
}

And the exact answer is hidden in the fortran files (in radfun.f called from radbas.f), where you can find the line

if( dtemp.lt.1e-20) dtemp =1e-20 

which treats small (even 0) values as 1e-20, which after taking square root results in 1e-10. It seems that the motivation was to speed up the process by using logarithm of the value (as a result, square rooting is simply dividing by 2), which of course is not defined for 0.