[R] alas, no vecnorm

[Douglas Bates]

Faheem Mitha <faheem at email.unc.edu> writes:

> I wanted a function that would give the euclidean distance of a vector. 
> 
> Then I was happy, because I found vecnorm listed on pg 55 of V&R (3rd
> edn) which I had just bought today.
> 
> Then I was sad, because R did not have vecnorm.
> 
> Then I was happy again, because I bethought myself that I could copy the
> function vecnorm from splus to my code.
> 
> Then I was sad again because R complained
> 
> > vecnorm(v)
> Error in .Fortran("d2norm", as.integer(length(x)), as.double(x), value =
> double(1)) : 
> 	C/Fortran function name not in load table

Check in src/appl/blas.f and you will find

*
*  dnrm2() returns the euclidean norm of a vector via the function
*  name, so that
*
*     dnrm2 := sqrt( x'*x )
*
*
*  -- This version written on 25-October-1982.
*     Modified on 14-October-1993 to inline the call to DLASSQ.
*     Sven Hammarling, Nag Ltd.
*
*
      double precision function dnrm2 ( n, x, incx )

Because this is a Fortran function and not a subroutine you can't call
it directly from R.  However, it would be an informative exercise to
take the manual "Writing R Extensions" and decide how to write a
Fortran or C wrapper function that can be called from R.

In C, for the .Call interface it could be written like 

#include <R.h>
#include <Rdefines.h>

extern double F77_NAME(dnrm2)(int *, double *, int *);

SEXP dnorm2(SEXP x)
{
    int one = 1, n;
    SEXP val;

    PROTECT(x = AS_NUMERIC(x)); n = LENGTH(x);
    PROTECT(val = NEW_NUMERIC(1));
    NUMERIC_POINTER(val)[0] =
        F77_CALL(dnrm2)(&n, NUMERIC_POINTER(x), &one);
    UNPROTECT(2);
    return(val);
}

In fact, for Intel machines and any others that use 80 bit floating
point registers, all the dancing around with scaling that is done in
dnrm2 is unnecessary.  The brute force approach of accumulating the
sum of the squares of the entries and taking the square root would
work perfectly fine.  However, since dnrm2 is available, we might as
well use it.

>From R the C routine dnorm2 would be called as
 .Call("dnorm2", x)

(Warning - I didn't test that code although I did compile it.)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._