R help - Oct 2007 - the use of the .C function

Dear All,
could someone please shed some light on the use of the .C  or .Fortran function:


I am trying load and running on R the  following function

// psi.cpp -- psi function for real arguments.

//      Algorithms and coefficient values from "Computation of Special

//      Functions", Zhang and Jin, John Wiley and Sons, 1996.

//

//  (C) 2003, C. Bond. All rights reserved.

//

//  Returns psi function for real argument 'x'.

//  NOTE: Returns 1e308 for argument 0 or a negative integer.

//

#include <iostream.h> // or <stdio.h>?

#include <math.h>

#define el 0.5772156649015329



double psi(double x)

{

    double s,ps,xa,x2;

    int n,k;

    static double a[] = {

        -0.8333333333333e-01,

         0.83333333333333333e-02,

        -0.39682539682539683e-02,

         0.41666666666666667e-02,

        -0.75757575757575758e-02,

         0.21092796092796093e-01,

        -0.83333333333333333e-01,

         0.4432598039215686};



    xa = fabs(x);

    s = 0.0;

    if ((x == (int)x) && (x <= 0.0)) {

        ps = 1e308;

        return ps;

    }

    if (xa == (int)xa) {

        n = xa;

        for (k=1;k<n;k++) {

            s += 1.0/k;

        }

        ps =  s-el;

    }

    else if ((xa+0.5) == ((int)(xa+0.5))) {

        n = xa-0.5;

        for (k=1;k<=n;k++) {

            s += 1.0/(2.0*k-1.0);

        }

        ps = 2.0*s-el-1.386294361119891;

    }

    else {

        if (xa < 10.0) {

            n = 10-(int)xa;

            for (k=0;k<n;k++) {

                s += 1.0/(xa+k);

            }

            xa += n;

        }

        x2 = 1.0/(xa*xa);

        ps = log(xa)-0.5/xa+x2*(((((((a[7]*x2+a[6])*x2+a[5])*x2+

            a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]);

        ps -= s;

    }

    if (x < 0.0)

        ps = ps - M_PI*cos(M_PI*x)/sin(M_PI*x)-1.0/x;

        return ps;

}


However, when applicated the codes
>digamma(-0.9)
[1] -9.312644  OK!

But when
> dyn.load("psi.so")
> out<-.C("psi", as.double(-0.9))
Erro en .C("psi", as.double(-0.9)) : C symbol name "psi" not
in load table
> out
Error: objeto "out" no encontrado

More information on OS:
> version
               _
platform       i486-pc-linux-gnu
arch           i486
os             linux-gnu
system         i486, linux-gnu
status
major          2
minor          4.1
year           2006
month          12
day            18
svn rev        40228
language       R
version.string R version 2.4.1 (2006-12-18)



Many thanks for any insight or comments.

Bernardo Lagos A.

Hi all,

here others doubts, when load and running the function

void gdot(double *x,double *y,int *n,double *output){
int i;
*output=0;
for (i=0;i<*n;i++){
*output+=x[i]*y[i];
}
}


as following

/examplesC$ R CMD SHLIB xby.c
gcc -std=gnu99 -I/usr/share/R/include -I/usr/share/R/include      -fpic  -g -O2
-c xby.c
-o xby.o
gcc -std=gnu99 -shared  -o xby.so xby.o   -L/usr/lib/R/lib -lR

/examplesC$ R
>x<-c(1,4,6,2)
>y<-c(3,2.4,1,9)
> dyn.load("xby.so")
>product<-.C("gdot",myx=as.double(x),myy=as.double(y),myn=as.integer(NROW(x)),myoutput=numeric())
 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: .C("gdot", myx = as.double(x), myy = as.double(y), myn =
as.integer(NROW(x)),
myoutput = numeric())

Possible actions:
1: abort (with core dump)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection:
Selection:
Selection:


Can anybody answerme that happening with this Traceback?


Thanking you in advance,
Bernardo Lagos A.

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1



The code is C++ and so is compiled
using the C++ compiler, not the C compiler.
This is because the name of the file is .cpp
(and you include iostream.h, but don't seem to make any use of it.)
As a result, the names of the routines are "mangled"
and so psi is not the actual name of the compiled routine,
but something more like _Z3psid which encodes the
types of the parameters in the name (to allow overloaded routines).

To avoid the mangling, use

#ifdef __cplusplus
extern "C"
#endif
double psi(double)

Alternatively, you can use the registration mechanism to
associate a name with a routine directly for R.  See
"Writing R Extensions" which would be a good thing to read
as trying to get this sort of thing working by trial and
error can be both very time consuming, and also very confusing
to the point that you "unlearn" things.

After that you will still need to make some changes
to use the routine with the .C() interface.
That expects inputs as pointers, i.e. the double
should be a double * and you fetch the value from
that. And a routine called by .C() should return nothing
directly, but provide the results in space provided
by the inputs

So
void R_psi(double *x)
{
   *x = psi(*x);
}

and make certain R_psi is not mangled.


Also, you might want to upgrade to the latest version of R.

Bernardo Lagos Alvarez wrote:
> Dear All,
> could someone please shed some light on the use of the .C  or .Fortran
function:
> 
> 
> I am trying load and running on R the  following function
> 
> // psi.cpp -- psi function for real arguments.
> 
> //      Algorithms and coefficient values from "Computation of Special
> 
> //      Functions", Zhang and Jin, John Wiley and Sons, 1996.
> 
> //
> 
> //  (C) 2003, C. Bond. All rights reserved.
> 
> //
> 
> //  Returns psi function for real argument 'x'.
> 
> //  NOTE: Returns 1e308 for argument 0 or a negative integer.
> 
> //
> 
> #include <iostream.h> // or <stdio.h>?
> 
> #include <math.h>
> 
> #define el 0.5772156649015329
> 
> 
> 
> double psi(double x)
> 
> {
> 
>     double s,ps,xa,x2;
> 
>     int n,k;
> 
>     static double a[] = {
> 
>         -0.8333333333333e-01,
> 
>          0.83333333333333333e-02,
> 
>         -0.39682539682539683e-02,
> 
>          0.41666666666666667e-02,
> 
>         -0.75757575757575758e-02,
> 
>          0.21092796092796093e-01,
> 
>         -0.83333333333333333e-01,
> 
>          0.4432598039215686};
> 
> 
> 
>     xa = fabs(x);
> 
>     s = 0.0;
> 
>     if ((x == (int)x) && (x <= 0.0)) {
> 
>         ps = 1e308;
> 
>         return ps;
> 
>     }
> 
>     if (xa == (int)xa) {
> 
>         n = xa;
> 
>         for (k=1;k<n;k++) {
> 
>             s += 1.0/k;
> 
>         }
> 
>         ps =  s-el;
> 
>     }
> 
>     else if ((xa+0.5) == ((int)(xa+0.5))) {
> 
>         n = xa-0.5;
> 
>         for (k=1;k<=n;k++) {
> 
>             s += 1.0/(2.0*k-1.0);
> 
>         }
> 
>         ps = 2.0*s-el-1.386294361119891;
> 
>     }
> 
>     else {
> 
>         if (xa < 10.0) {
> 
>             n = 10-(int)xa;
> 
>             for (k=0;k<n;k++) {
> 
>                 s += 1.0/(xa+k);
> 
>             }
> 
>             xa += n;
> 
>         }
> 
>         x2 = 1.0/(xa*xa);
> 
>         ps = log(xa)-0.5/xa+x2*(((((((a[7]*x2+a[6])*x2+a[5])*x2+
> 
>             a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]);
> 
>         ps -= s;
> 
>     }
> 
>     if (x < 0.0)
> 
>         ps = ps - M_PI*cos(M_PI*x)/sin(M_PI*x)-1.0/x;
> 
>         return ps;
> 
> }
> 
> 
> However, when applicated the codes
>> digamma(-0.9)
> [1] -9.312644  OK!
> 
> But when
>> dyn.load("psi.so")
>> out<-.C("psi", as.double(-0.9))
> Erro en .C("psi", as.double(-0.9)) : C symbol name
"psi" not in load table
> 
>> out
> Error: objeto "out" no encontrado
> 
> More information on OS:
>> version
>                _
> platform       i486-pc-linux-gnu
> arch           i486
> os             linux-gnu
> system         i486, linux-gnu
> status
> major          2
> minor          4.1
> year           2006
> month          12
> day            18
> svn rev        40228
> language       R
> version.string R version 2.4.1 (2006-12-18)
> 
> 
> 
> Many thanks for any insight or comments.
> 
> Bernardo Lagos A.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.7 (Darwin)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iD8DBQFHEQMO9p/Jzwa2QP4RAiupAJ40LfoWt3Jm6dPXU0aYHe3gXPAhQgCdFAig
VScTISDpaa4pzGw+eNkkxKY=Kd5o
-----END PGP SIGNATURE-----