Skip to content
Prev 5518 / 398506 Next

alas, no vecnorm

Faheem Mitha <faheem at email.unc.edu> writes:
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.
.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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._