Re: compiling C++ and C functions for octave

[Geraint]

"Iago Mosqueira" <address@hidden> writes:

> Just one question. For this code to accept vectors, would it need to run a
> loop? Or can the octave C++ functions take care of that?
> 
> Thanks,
> 
> iago
> 

It can accept vectors with a slight modification (you might want to improve the 
error handling, check that all the arguments are the same length, etc.). I 
don't know anything about the R mathlib, but if the R functions accept vector 
inputs you could probably eliminate the for loop in the DLD.

/Geraint.


octave:1> matern(1,2,3)
ans = 0.96965
octave:2> matern(2,3,4)
ans = 0.96396
octave:3> matern([1;2],[2;3],[3;4])
ans =

  0.96965
  0.96396


matern.cc: (compile with 'mkoctfile -lRmath matern.cc')

#define MATHLIB_STANDALONE
#include <Rmath.h>
#include <math.h>
#include <octave/oct.h>

extern "C"
{
  float matern(float x, float phi, int kappa)
  {
    /*Returns the Matern function for x, phi and kappa.*/

    /* Variables */
    float ans, cte;
    float uphi=x/phi;

    /* Matern */

    if (uphi==0) return 1;
    else{
      if (kappa==0.5) 
        ans = exp(-uphi);
      else {
        cte = R_pow(2, (-(kappa-1)))/gammafn(kappa); 
        ans = cte * R_pow(uphi, kappa) * bessel_k(uphi,kappa,1); 
      }
    }

    return ans;
  }
}

DEFUN_DLD (matern, args, ,
           "Matern variogram function")
{

  /* Number of arguments*/
  int nargs = args.length ();
        
  /* Error checking*/
  if (nargs!=3){
    usage("matern (x, phi, kappa)");
    return octave_value(0.0);
  }

  ColumnVector x = args(0).column_vector_value();
  ColumnVector phi = args(1).column_vector_value();
  ColumnVector kappa = args(2).column_vector_value();
  
  int n = x.length();
  ColumnVector retval = ColumnVector(n);
  for( int i = 0; i < n; i++ ){
    retval(i) = matern (x(i), phi(i), static_cast<int>(kappa(i)));
  }  

  return retval;
}

                                        




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------